/*	Migmovie.c

	Migration in the (z, x, omega)-domain at a range of velocities to
	produce a movie of the migration process.  Altered from wmig.c .
	From Claerbout (1985) p. 111.

	Usage:		migmovie >outfile

	Compile:	cc migmovie.c ctris.c colmft.c -fswitch -lm -o migmovie
									*/
/* The usual includes */
#include <stdio.h>
#include <math.h>

/* Define complex multiplication and conjugate for declared struct complex */
#define  rmul(x,y)      (x.re * y.re - x.im * y.im)
#define  imul(x,y)      (x.im * y.re + x.re * y.im)
#define rcmul(x,y)      (x.re * y.re + x.im * y.im)
#define icmul(x,y)      (x.im * y.re - x.re * y.im)

/* Define complex division -- x/y for declared struct complex */
#define  rdiv(x,y) ((x.re * y.re + x.im * y.im)/(y.re * y.re + y.im * y.im))
#define  idiv(x,y) ((x.im * y.re - x.re * y.im)/(y.re * y.re + y.im * y.im))

/* Declare the structure to hold complex numbers */
struct complex 
	{
	float re; 
	float im;
	};

main()
	{
	/* Declare all variables */
	int ix, nx, iz, nz, iw, nw, it, nt, iv, nv;
	float q[64][48], pi2, alpha, dt, dtau, dw, v, dv;
	struct complex cq[64][48], cd[48], ce[48], cf[48], aa, a, b, c, cshift;
	struct complex num, den, cq0[64][48];

	/* Initialize constants */
	nt = 64;
	nz = nt;
	nx = 48;
	pi2 = 2.0*3.141592;
	dt = 1;
	dtau = 1;
	dw = pi2/(dt*nt);
	nw = nt/2;
	/* alpha = v*v*dtau/(4*dx*dx) */
	alpha = 0.25;
	nv = 11;
	dv = 1.0/(nv-1);

	/* Zero out data volume */
	for (iz=0; iz<nz; iz++)
		for (ix=0; ix<nx; ix++)
			{
			cq0[iz][ix].re = 0;
			cq0[iz][ix].im = 0;
			}

	/* Dipping line segment source */
	for (it=nt/3; it<nt; it+=nt/4)
		{
		ix = nx/4 - 1;
		cq0[it-ix/2][ix].re = 2;
		cq0[it-ix/2+1][ix].re = 2;
		for (ix=nx/4; ix<nx/2; ix++)
			{
			cq0[it-ix/2][ix].re = 4;
			cq0[it-ix/2+1][ix].re = 4;
			}
		cq0[it-ix/2][ix].re = 2;
		cq0[it-ix/2+1][ix].re = 2;
		}
	/* Broadened impulse source
	for (it=nt/3; it<nt; it+=nt/4)
		for (ix=0; ix<4; ix++)
			{
			cq0[it][ix].re = 4 - ix;
			cq0[it+1][ix].re = 4 - ix;
			}
	*/

	/* FT over time, along the columns */
	colmft(nt, nx, cq0, 1.0, 1.0);

	/* Increment migration velocity up to the true value */
	for (iv=0; iv<nv; iv++)
		{
		/* Copy transformed array into working array */
		for (iz=0; iz<nz; iz++)
			for (ix=0; ix<nx; ix++)
				{
				cq[iz][ix].re = cq0[iz][ix].re;
				cq[iz][ix].im = cq0[iz][ix].im;
				q[iz][ix] = 0;
				}
		v = iv*dv;
		v = v*v;
		/* iz and iw loops interchangeable */
		for (iz=0; iz<nz; iz++)
			for (iw=1; iw<nw; iw++)
				{
				/* Tridiagonal matrix coefficients */
				/* aa = -v*alpha/(-i*iw*dw) */
				num.re = -v*alpha;
				num.im = 0;
				den.re = 0;
				den.im = -iw*dw;
				aa.re = rdiv(num,den);
				aa.im = idiv(num,den);
				/* Diagonals of the matrix */
				a.re = -aa.re;
				a.im = -aa.im;
				b.re = 1 + 2*aa.re;
				b.im = 2*aa.im;
				c.re = -aa.re;
				c.im = -aa.im;
	
				/* Apply finite difference star */
				for (ix=1; ix<nx-1; ix++)
					{
					/* cd[ix]=aa*cq[iw][ix+1]
					+(1-2*aa)cq[iw][ix] +aa*cq[iw][ix-1] */
					num.re = 1 - 2*aa.re;
					num.im = -2*aa.im;
					/* Can only put one complex
					multiplication macro on each line */
					den.re = rmul(aa,cq[iw][ix-1]);
					den.im = imul(aa,cq[iw][ix-1]);
					den.re = rmul(num, cq[iw][ix]) + den.re;
					den.im = imul(num, cq[iw][ix]) + den.im;
					cd[ix].re=rmul(aa,cq[iw][ix+1])+den.re;
					cd[ix].im=imul(aa,cq[iw][ix+1])+den.im;
					}
	
				/* Set side conditions */
				cd[0].re = 0;
				cd[0].im = 0;
				cd[nx-1].re = 0;
				cd[nx-1].im = 0;
	
				/* Set the left and right ends of the diagonal*/
				/* Put endl into a complex buffer */
				num.re = -a.re;
				num.im = -a.im;
				/* Put endr into a complex buffer */
				den.re = -c.re;
				den.im = -c.im;
				/* Solve the complex tridiagonal system */
				if (v != 0)
					ctris(nx, &num, &a, &b, &c, &den, cd,
					&cq[iw][0], ce, cf);
	
				/* cshift = e^(-i*iw*dw*dtau) */
				cshift.re = (float)(cos((double)(-iw*dw*dtau)));
				cshift.im = (float)(sin((double)(-iw*dw*dtau)));
	
				/* Shifting term */
				for (ix=0; ix<nx; ix++)
					{
					num.re = rmul(cq[iw][ix],cshift);
					num.im = imul(cq[iw][ix],cshift);
					cq[iw][ix].re = num.re;
					cq[iw][ix].im = num.im;
					}
				
				/* Imaging condition: q(t=0) = Sum Q(w) */
				for (ix=0; ix<nx; ix++)
					q[iz][ix] += cq[iw][ix].re;
	
				}
	
		/* Write data array to standard output in rows of z,
			one at a time */
		for (ix=0; ix<nx; ix++)
			for (iz=0; iz<nz; iz++)
				write(1, &q[iz][ix], sizeof(float));
		
		}

	/* All done */
	exit(0);
	}