/*	wmig.java

	Migration in the (z, x, omega)-domain.
	Adapted from Jon Claerbout's Imaging the Earth's Interior (1985), p. 111,
	by John N. Louie, 29 November 1998.
*/

/* Include the usual classes */
import java.util.*;
import java.io.*;
import java.lang.Math;

/* Declare a class to contain wmig as a main method, referring to the
   class fttestc.colmft() static method, and re-defining the extrap.ctris()
   static method. */
/* When you change the code in this file, make a copy to a new file with
    a new name like "ex8wmig.java", and be sure to make the same change
    to the class name below - e.g. "class ex8wmig" */
class wmig
	{
	/* The main() method compiles into a file named "wmig.class" that
	    you run in the Java virtual machine as a stand-alone application
	    with a command such as "java wmig"; or whatever you have
	    renamed the class and file. */
	public static void main(String args[]) /* No runtime arguments needed */
		{
		System.out.println("Starting wmig...");
		/* Declare variables */
		int ix, nx, iz, nz, iw, nw, it, nt;
		float pi2, alpha, dt, dtau, dw, s, cang, w, v, dx;
		Complex aa = new Complex();
		Complex a = new Complex();
		Complex b = new Complex();
		Complex c = new Complex();
		Complex cshift = new Complex();
		Complex num = new Complex();
		Complex den = new Complex();
		Complex num2 = new Complex();
		Complex bl = new Complex();
		Complex br = new Complex();
		Complex qm1 = new Complex();
		Complex qnx = new Complex();
		
		/* Initialize constants */
		nt = 64;
		nz = nt;
		nx = 48;
		pi2 = (float)(2.0*Math.PI);
		dt = 1F;
		dtau = 1F;
		dw = pi2/(dt*nt);
		nw = nt/2;
		v = 1F;
		dx = 1F;
		/* alpha = 0.25; */
		alpha = v*v*dtau/(4F*dx*dx);
		
		/* Declare and initialize planes to zero values */
		FltPlane q = new FltPlane(nt, nx);
		ComplexPlane cq = new ComplexPlane(nt, nx);
		/* Note that wmig operates in rows of x, for the
		   convenience of the ctrisr() method. */
		ComplexVec cd = new ComplexVec(nx);
		ComplexVec ce = new ComplexVec(nx);
		ComplexVec cf = new ComplexVec(nx);
		
		/* Broadened impulse source */
		for (it=nt/3; it<nt; it+=nt/4)
			/* Tricky sources on sides */
			for (ix=0; ix<4; ix++)
				{
				cq.vec[it*cq.ne + ix].re = 4 - ix;
				cq.vec[(it+1)*cq.ne + ix].re = 4 - ix;
				}

		/* FT over time, along the columns */
		fttestc.colmft(cq, 1F, 1F);
		
		/* iz and iw loops interchangeable */
		for (iz=0; iz<nz; iz++)
			for (iw=1; iw<nw; iw++)
				{
				w = iw*dw;
				/* Tridiagonal matrix coefficients */
				/* aa = -alpha/(-i*iw*dw) */
				num.re = -alpha;
				num.im = 0F;
				den.re = 0F;
				den.im = -iw*dw;
				aa.div(num,den);
				/* Diagonals of the matrix */
				a.re = -aa.re;
				a.im = -aa.im;
				b.re = 1F + 2F*aa.re;
				b.im = 2F*aa.im;
				c.re = -aa.re;
				c.im = -aa.im;
				/* Left and right boundary conditions */
				bl.re = 1F;
				bl.im = 0F;
				br.re = 1F;
				br.im = 0F;
				
				/* 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 = 1F - 2F*aa.re;
					num.im = -2F*aa.im;
					/* Can only put one complex multiplication
					   method on each line */
					den.mul(aa,cq.vec[iw*cq.ne + ix-1]);
					num2.mul(num, cq.vec[iw*cq.ne + ix]);
					den.add(num2,den);
					cd.vec[ix].mul(aa,cq.vec[iw*cq.ne + ix+1]);
					cd.vec[ix].add(cd.vec[ix],den);
					}

				/* Set side conditions */
				qm1.mul(bl,cq.vec[iw*cq.ne + 0]);
				qnx.mul(br,cq.vec[iw*cq.ne + nx-1]);
				/* cd[0]=aa*cq[iw][1]+(1-2aa)cq[iw][0]+aa*cq[iw][-1] */
				num.re = 1 - 2*aa.re;
				num.im = -2*aa.im;
				/* Can only put one complex multiplication
				   method on each line */
				num2.mul(aa,qm1);
				cd.vec[0].mul(num, cq.vec[iw*cq.ne + 0]);
				cd.vec[0].add(cd.vec[0],num2);
				num2.mul(aa,cq.vec[iw*cq.ne + 1]);
				cd.vec[0].add(cd.vec[0],num2);
				/* cd[nx-1]=aa*cq[iw][nx]+(1-2aa)cq[iw][nx-1]+aa*cq[iw][nx-2] */
				num2.mul(aa,cq.vec[iw*cq.ne + nx-2]);
				cd.vec[nx-1].mul(num, cq.vec[iw*cq.ne + nx-1]);
				cd.vec[nx-1].add(cd.vec[nx-1],num2);
				num2.mul(aa,qnx);
				cd.vec[nx-1].add(cd.vec[nx-1],num2);
				
				/* Set the left and right ends of the diagonal */
				/* Put endl into a complex buffer */
				num.mul(c,bl);
				num.add(num, b);
				/* Put endr into a complex buffer */
				num2.mul(a,br);
				num2.add(num2, b);
				/* Solve the complex tridiagonal system, for a
				   row of the class ComplexPlane at cq[iw*ne+0] */
				ctrisr(nx,num,a,b,c,num2,cd,cq,iw*cq.ne,ce,cf);
				
				/* cshift = e^(-i*iw*dw*dtau) */
				cshift.re = (float)(Math.cos((double)(-iw*dw*dtau)));
				cshift.im = (float)(Math.sin((double)(-iw*dw*dtau)));

				/* Shifting term */
				for (ix=0; ix<nx; ix++)
					cq.vec[iw*cq.ne + ix].mul(cq.vec[iw*cq.ne + ix],cshift);
					
				/* Imaging condition: q(t=0) = Sum Q(w) */
				for (ix=0; ix<nx; ix++)
					q.vec[iz*q.ne + ix] += cq.vec[iw*cq.ne + ix].re;
				}

		/* Plot data plane in new frame, after setting label info */
		float rms = q.rms();
		q.setAmpUnits("Reflectivity");
		q.setUnitsofDim(0, "Distance X");
		q.setDUnitofDim(0, dx);
		q.setUnitsofDim(1, "Depth Z");
		q.setDUnitofDim(1, 1F);
		q.viewinFrame("Lab 4: Freq.-domain Migration", false, 1F,
			RgCtab.vagray, false, 2F*rms);
		
		/* All Done with wmig's main() application method */
		}

/*	ctrisr.java

		Complex tridiagonal equation solver for the system:

		|l c 0 . . . . . 0| | q[0] |   | d[0] |
		|a b c 0 . . . . .| |  .   |   |  .   |
		|0 a b c 0 . . . .| |  .   |   |  .   |
		|. .` ` `  . . . .| |  .   |   |  .   |
		|. . .` ` `  . . .| |  .   | = |  .   |
		|. . . .` ` `  . .| |  .   |   |  .   |
		|. . . . 0 a b c 0| |  .   |   |  .   | 
		|. . . . . 0 a b c| |  .   |   |  .   |
		|0 . . . . . 0 a r| |q[n-1]|   |d[n-1]|

		Adapted from Claerbout (1985) p. 98-99,
		by John N. Louie, 28 November 1998.

		Arguments:	n	dimension of arrays
				endl	corner complex value l (above)
				a,b,c	diagonals of the complex
					coefficient matrix
				endr	corner complex value r (above)
				d	ComplexVec data vector of n complex elements
				q	ComplexPlane solution plane, with vectors of n complex elements
				w0	index into ComplexPlane q of vector for solution
				e,f	ComplexVec temporary iteration vectors of
					n complex elements
										*/
	static void ctrisr(int n, Complex endl, Complex a, Complex b, Complex c,
		Complex endr, ComplexVec d, ComplexPlane q, int w0,
		ComplexVec e, ComplexVec f)
		{
		int i;
		Complex den = new Complex();
		Complex num = new Complex();
		
		/* Initialize unknowns e and f */
		e.vec[0].div(a,endl);
		e.vec[0].scale(-1F);
		f.vec[0].div(d.vec[0],endl);
		
		/* Recursions of equations 19a and 19b */
		for (i=1; i<n-1; i++)
			{
			den.mul(c,e.vec[i-1]);
			den.add(b,den);
			e.vec[i].div(a,den);
			e.vec[i].scale(-1);
			num.mul(c,f.vec[i-1]);
			num.sub(d.vec[i],num);
			f.vec[i].div(num,den);
			}

		/* Right-hand boundary condition, equation 20 */
		num.mul(c,f.vec[n-2]);
		num.sub(d.vec[n-1],num);
		den.mul(c,e.vec[n-2]);
		den.add(endr,den);
		q.vec[w0+n-1].div(num,den);
		
		/* Recursion for equation 21 */
		for (i=n-2; i>=0; i--)
			{
			num.mul(e.vec[i],q.vec[w0+i+1]);
			num.add(num,f.vec[i]);
			q.vec[w0+i].re = num.re;
			q.vec[w0+i].im = num.im;
			}
		/* End of ctrisr static method of class wmig */
		}

	/* End of methods in class wmig */
	}