/* *********************** FROGTR.C *************************

This program calculates a FROG trace for an ultrashort laser pulse with user-specified input parameters and draws it on the screen.  You may download it and use it in any manner you like.

This program must be compiled and run on an ASCII C compliant machine
with X Windows capabilities.  On most UNIX machines, the line
		cc -o frogtr -O2 frogtr.c -lm -lX11
should work fine.  Those of you who know X Windows programming will
be shocked by the blatant flaunting of conventional standards.  But
hey, it works.

Written by Ken DeLong, July 26, 1995.  Email to Ken can be directed to
71022.2774@compuserve.com, or kwdelon@bullfrog.ca.sandia.gov.  See also
our Web Page at http://www.ca.sandia.gov/ultrafrog.

The program will ask for several input parameters.  Ranges are suggested.
There are only two parameters that the user should change: the number
of points N, and the FROG geometry.  Both of these are described below.
*/

/* Header files and function prototypes. */
#include 
#include 
#include 
enum {PG, SD, SHG};
void four1(double dum[], int nn, int isign);
void DLoadFFT(double *re, double *im, double *x, int n);
void DUnloadFFT(double *re, double *im, double *x, int n);
double *set_vec(int n);
double **set_mat(int n);
void initgraph(void);
void initwindow(int win, int x, int y, int width, int height, char *title);
void closegraph(void);
void DrawGraph(int N, double y[], int win);
void ClearGraph(int win);
void Color(unsigned int color);
void PrintGraph(void);
void DrawImage(double *data, int N, int win);
void DrawText2(char *str, int win);
void OpenWindows(int m);
void DrawFrog(double **d, int win);
char Getch(void);

int N;

void main(void) {

	double *x, *Er, *Ei, **Esigr, **Esigi, **frogI;
	int flag, w, f, tau, t, tprime, Geometry, i;
	double amp, delw, A, CS, S4, arg;

	/* This is the size of the array.  Feel free to change this,
	but it must be a power of 2!! */
	N = 64;

	/* The geometry can be Polarization Gate (PG), Second-Harmonic
	Generation (SHG), or Self-Diffraction (SD).  See the paper
	DeLong, et al, JOSA B v. 11, p. 1595 (1994). */
	Geometry = PG;

	/* Setup for windows. */
	initgraph();
	OpenWindows(1);
	printf("Hit return. . .\n");
	Getch();

	x = set_vec(2*N);
	Er = set_vec(N);
	Ei = set_vec(N);
	Esigr = set_mat(N);
	Esigi = set_mat(N);
	frogI = set_mat(N);

	if(Geometry == PG) printf("Polarization Gate geometry.\n\n");
	if(Geometry == SHG) printf("Second-Harmonic Generation geometry.\n\n");
	if(Geometry == SD) printf("Self-Diffraction geometry.\n\n");

	printf("Spectral width (2 to 10)? "); scanf("%lg", &delw);
	printf("Quadratic Phase (-6 to 6)? "); scanf("%lg", &A);
	printf("Cubic Phase (-4 to 4)? "); scanf("%lg", &CS);
	printf("Quartic Phase (-2 to 2)? "); scanf("%lg", &S4);

	/* Generate the field with the 4 parameters */
	flag = 1;
	for(w = 0; w < N; w++) {
		f = w - N/2;
		amp = exp(-2.0*log(2.0)*f*f/(delw*delw));
		arg = A*pow((f/delw), 2.0) + CS*pow((f/delw), 3.0)
						 + S4*pow((f/delw),4.0);
		x[2*w] = amp*cos(arg)*flag;
		x[2*w+1] = amp*sin(arg)*flag;
		flag *= -1;
	}
	four1(x, N, -1);
	DUnloadFFT(Er, Ei, x, N);

	/* Generate a FROG trace */
	/* Make signal field */
	for(tau = 0; tau < N; tau++) {
		double inten;
		for(t = 0; t < N; t++) {
			tprime = t - (tau - N/2);  /* N/2 since we want the center of the array to be tau = 0 */
			if(tprime < 0 || tprime > N-1) {
				Esigr[tau][t] = Esigi[tau][t] = 0.0;
				continue;
			}
			if (Geometry == PG ) {
				inten = Er[tprime]*Er[tprime]
					+ Ei[tprime]*Ei[tprime];
				Esigr[tau][t] = Er[t]*inten;
				Esigi[tau][t] = Ei[t]*inten;
			}
			if (Geometry == SD) {
				double u, v;
				u = Er[t]*Er[t] - Ei[t]*Ei[t];
				v = 2*Er[t]*Ei[t];
				Esigr[tau][t] = u*Er[tprime] + v*Ei[tprime];
				Esigi[tau][t] = v*Er[tprime] - u*Ei[tprime];
			}
			if (Geometry == SHG) {
				Esigr[tau][t] = Er[t]*Er[tprime] - Ei[t]*Ei[tprime];
				Esigi[tau][t] = Er[t]*Ei[tprime] + Ei[t]*Er[tprime];
			}
		}
	}

	/* Now make the FROG signal by Fourier transform along t */
	for(tau = 0; tau < N; tau++) {
		flag = 1;
		for(t = 0; t < N; t++) {
			x[2*t] = Esigr[tau][t]*flag;
			x[2*t+1] = Esigi[tau][t]*flag;
			flag *= -1;
		}
		four1(x, N, 1);
		flag = 1;
		for(w = 0; w < N; w++) {
			i = 2*w;
			frogI[tau][w] = x[i]*x[i] + x[i+1]*x[i+1];
			flag *= -1;
		}
	}
	ClearGraph(1);
	DrawFrog(frogI, 1);
	PrintGraph();
	printf("Hit return. . .\n");
	Getch();


	free(frogI[0]); free(frogI);
	free(Esigi[0]); free(Esigi);
	free(Esigr[0]); free(Esigr);
	free(Ei);
	free(Er);
	free(x);

	closegraph();
}

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void four1(double dum[], int nn, int isign) {

	int n, mmax, m, j, istep, i;
	double wtemp, wr, wpr, wpi, wi, theta;
	double tempr, tempi, *data;

	data = dum - 1;
	n = nn << 1;
	j = 1;
	for (i = 1; i < n; i += 2) {
		if (j > i) {
			SWAP(data[j], data[i]);
			SWAP(data[j+1], data[i+1]);
		}
		m = n >> 1;
		while (m >= 2 && j > m) {
			j -= m;
			m >>= 1;
		}
		j += m;
	}
	mmax=2;
	while (n > mmax) {
		istep=2*mmax;
		theta=6.28318530717959/(isign*mmax);
		wtemp=sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi=sin(theta);
		wr=1.0;
		wi=0.0;
		for (m=1;m= 0; w--)
		for(tau = 0; tau < N; tau++)
			 data[i++] = sqrt(d[tau][w]);
	DrawImage(data, N, win);
	free(data);
}

#include 
#include 
#include 
#include 
#include 
#include 

#define PI 3.141592654    /* set up PI defintion */

#define MAXW 5

Font font;
Display *mydisplay;
Window mywindow[MAXW];
GC mygc;
unsigned long myforeground, mybackground;
int count = 0;

/* Initializes everything and sets up all the pointers */
void initgraph(void)
{
	XEvent myevent;
	KeySym mykey;
	int myscreen;

	/* set up display and foreground/background */

	mydisplay = XOpenDisplay("");

	myscreen = DefaultScreen (mydisplay);
	mybackground = WhitePixel (mydisplay, myscreen);
	myforeground = BlackPixel (mydisplay, myscreen);
}

/* Creates a window number win, at (x,y), with title */
void initwindow(int win, int x, int y, int width, int height, char *title) {

	int argc;
	char *argv[1];
	XSizeHints myhint;

	if(win > MAXW) {
		printf("\nToo many windows!!  Increase MAXW in xgraph.c\n");
		exit(0);
	}

	/* create a window with a width of WIDTH x HEIGHT */
	argc = 1; argv[0] = "frog";

	myhint.x = x;
	myhint.y = y;
	myhint.width = width;
	myhint.height = height;
	myhint.flags = PPosition | PSize;
	mywindow[win - 1] = XCreateSimpleWindow (mydisplay,
	    DefaultRootWindow (mydisplay),
	    myhint.x, myhint.y, myhint.width, myhint.height,
	    5, myforeground, mybackground);

	XSetStandardProperties (mydisplay, mywindow[win - 1], title, title,
	    None, argv, argc, &myhint);

	if(win == 1) {
		mygc = XCreateGC (mydisplay, mywindow[win-1], 0, 0);

		/* Load a font called 8x16 */

		font = XLoadFont(mydisplay,"8x16");
		XSetFont(mydisplay,mygc,font);

		XSetBackground (mydisplay, mygc, mybackground);
		XSetForeground (mydisplay, mygc, myforeground);
	}

	XSelectInput (mydisplay, mywindow[win - 1],
	    ButtonPressMask | KeyPressMask | ExposureMask);

	/* map window to the screen */

	XMapRaised (mydisplay, mywindow[win - 1]);
	XFlush(mydisplay);
	count++; /* count keeps track of the number of windows */
}

/* Turns off all the X Windows stuff */
void closegraph(void) {
	int i;

	XFreeGC (mydisplay, mygc);
	for(i = 0; i < count; i++)
		XDestroyWindow (mydisplay, mywindow[i]);
	XCloseDisplay (mydisplay);
}
void ClearGraph(int win) {

	XClearWindow(mydisplay, mywindow[win - 1]);
}
void PrintGraph(void) {
	XFlush(mydisplay);
}
unsigned long pixels[101], table = 0;
void DrawImage(double *datar, int N, int win) {

	XImage *image;
	char *idata;
	int i, err, fact, ii, j, k, mag, row, col;
	XWindowAttributes *xwa;
	double max;
	void CreateColorTable(XWindowAttributes *xwa);

	/* Allocate memory, check display depth */
	xwa = (XWindowAttributes *) malloc(sizeof(XWindowAttributes));
	XGetWindowAttributes(mydisplay, mywindow[win - 1], xwa);
	if(xwa->depth != 8) printf("Trouble: display depth %d\n", xwa->depth);

	/* Create the color table */
	if(table == 0) {
		CreateColorTable(xwa);
		/*unsigned long pm[2];
		unsigned int np, npix;
		Bool contig;
		XStandardColormap map;
		err=XGetStandardColormap(mydisplay, RootWindow(mydisplay, 0),
									&map, XA_RGB_BEST_MAP);
		printf("Return of XGetStandardColormap is %d\n", err);
		printf("map.colormap is %d\n", map.colormap);
		np = 0; npix = 100; contig = False;
		err=XAllocColorCells(mydisplay, map.colormap, contig, pm, np, pixels, npix);
		printf("Return of XAllocColorCells is %d\n", err);*/
		table = 1;
	}

	/* Find maximum of array */
	for(i = 0, max = 0; i < N*N; i++)
			if(datar[i] > max) max = datar[i];

	/* Try to fill the window with the image.  mag is magnification */
	fact = (int) pow(2, 8*sizeof(char) );   /* Scaling factor for depth*/
	fact = 100;
	mag = (int) (xwa->width < xwa->height ? xwa->width : xwa->height)/N;
	if(mag == 0) mag = 1;
	idata = (char *) malloc(mag*mag*N*N*sizeof(char) );
	for(row = 0; row < N; row++) {
		for (col = 0; col < N; col++) {
			for(j = 0; j < mag; j++){
				for(k = 0; k < mag; k++) {
					i = N*row + col;
					ii = row*N*mag*mag + mag*col + N*j*mag +k;
					idata[ii] = (char) pixels[(int) (fact*datar[i]/max)];
				}
			}
		}
	}

	image = XCreateImage( mydisplay, xwa->visual, 8, ZPixmap, 0,
		idata, N*mag, N*mag, 8*sizeof(char), 0 );

	err = XPutImage( mydisplay, mywindow[win-1], mygc, image, 0, 0, 0, 0,
		 mag*N, mag*N);

	free(idata);
	free(xwa);
}

/* Makes the custom FROG color table. */
void CreateColorTable(XWindowAttributes *xwa) {

	double crud;
	int i, k, j;
	char cname[12], junk[4];
	XColor *color;

	color = (XColor *) malloc(sizeof(XColor));

	/* Find the pixels for a Gray Scale display */
	for(i = 0; i <= 100; i++) {
		crud = i; crud /= 1000;
		/* This figures out how many leading zeros */
		for(j = 0, k = 0; j < 3; j++) {
			crud *= 10;
			if( (int) crud == 0 && k == 0) continue;
			junk[k++] = '0' + (int) crud;
			crud -= (int) crud;
		}
		junk[k] = '\0';
		strcpy(cname, "Gray");
		strcat(cname, junk);
		if(k == 0) strcpy(cname, "Gray0");
		if(i==0 || i==100) XParseColor(mydisplay, xwa->colormap, cname, color);
		else {
			/*color->red = 65535*exp(-0.001109*i*i/2);
			color->green = 65535*exp(-0.01109*(i-50)*(i-50)/2);
			color->blue = 65535*exp(-0.001109*(i-100)*(i-100)/2);*/
			double Step1 = 15, Step2 = 50, Step3 = 95;
			double red, green, blue;
			/* A smaller exponent means a faster rise, slower fall */
			if(i <= Step1) {
				red = 1.0;
				green = (1.0*i)/Step1;
				green = pow(green, 0.7);
				blue = 0.0;
			}
			if(i > Step1 && i <= Step2) {
				red = (Step2 - i)/(Step2 - Step1);
				red = pow(red, 0.7);
				green = 1.0;
				blue = 0.0;
			}
			if(i > Step2 && i <= Step3) {
				red = 0.0;
				green = (Step3 - i)/(Step3 - Step2);
				blue = (i - Step2)/(Step3 - Step2);
				green = pow(green, 0.8);
				blue = pow(blue, 0.8);
			}
			if(i > Step3) {
				red = green = (i - Step3)/(100 - Step3);
				red = pow(red, 1.0); green = pow(green, 1.0);
				blue = 1.0;
			}
			color->red = red*65535;
			color->green = green*65535;
			color->blue = blue*65535;
		}
		XAllocColor(mydisplay, xwa->colormap, color);
		pixels[i] = color->pixel;
	}
	free(color);
}
void OpenWindows(int m) {

#define WIDTH 290
#define HEIGHT 290

	int i, x, y;
	char title[35];

	/* Allocate windows, numbered from 1 to m */
	for(i = 1; i <= m; i++) {
		switch(i) {
			case 1: x = 0;
				y = 0;
				strcpy(title, "The FROG Trace");
				break;
			case 2: x = 0;
				y = 310;
				strcpy(title, "Retrieved FROG Trace");
				break;
			case 3: x = 400;
				y = 0;
				strcpy(title, "Spectrum");
				break;
			case 4: x = 400;
				y = 310;
				strcpy(title, "Experimental FROG Trace");
				break;
			default: x = 200;
				y = 200;
				strcpy(title, "Unknown quantity");
				break;
		}
		initwindow(i, x, y, WIDTH, HEIGHT, title);
	}
}
char Getch(void) {

	char c;

	fflush(stdin);
	c = getchar();
	return c;
}