/*	tdif.java

	Time domain 15-degree diffraction movie.
	Star:	w=p(t,z)	y=p(t,z+1)
		u=p(t+1,z)	v=p(t+1,z+1)
	Adapted from Jon Claerbout's Imaging the Earth's Interior (1985), p. 136,
	by John N. Louie, 30 November 1998.
*/

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

/* Declare a class to contain tdif as a main method, referring to the
   class imhf.rtris() static method. */
/* When you change the code in this file, make a copy to a new file with
    a new name like "ex2tdif.java", and be sure to make the same change
    to the class name below - e.g. "class ex2tdif" */
class tdif
	{
	/* The main() method compiles into a file named "tdif.class" that
	    you run in the Java virtual machine as a stand-alone application
	    with a command such as "java tdif"; or whatever you have
	    renamed the class and file. */
	public static void main(String args[]) /* No runtime arguments needed */
		{
		System.out.println("Starting tdif...");
		/* Declare variables */
		int i, ix, nx, iz, nz, it, nt, ip;
		float alfa, beta, diag, offdi, apb, amb, dd;
		
		/* Initialize constants */
		nx = 36;
		nz = 96;
		nt = 96;
		/* v*dz*dt/(8*dx*dx) */
		alfa = 0.125F;
		/* Accurate x derivative parameter; simplest case b=0 */
		beta = 0.140F;

		/* Declare and initialize volume and vectors to zero values */
		FltVol p = new FltVol(nz/2+1, nz, nx+1);
		/* Note that tdif operates in rows of x, for the
		   convenience of the rtris() method. The number of planes
		   is linked to the first loop below (and to the velocity used);
		   and there is an extra column for a depth pointer. */
		float[] u = new float[nx];
		float[] w = new float[nx];
		float[] v = new float[nx];
		float[] y = new float[nx];
		float[] e = new float[nx];
		float[] f = new float[nx];
		float[] d = new float[nx];
		
		/* Set up initial model of 4 band-limited "point" scatterers */
		for (iz=nz/5-1; iz<nz; iz+=nz/4)
			for (it=0; it<15; it++)
				for (ix=0; ix<4 && iz+it<nz; ix++)
					p.vec[0*p.nv*p.ne + (iz+it)*p.ne + ix]
					   = (4F-ix)*(7F-it)
					   * (float)(Math.exp(-0.1*(it-7.0)*(it-7.0)));

		/* Tridiagonal coefficients */
		apb = alfa + beta;
		amb = alfa - beta;
		diag = 1F + 2F*amb;
		offdi = -amb;
		
		/* Climb up in steps of 2 z-levels with iz.
		   ip keeps track of current plane of FltVol p */
		for (iz=nz-1, ip=1; iz>0 && ip<p.planes; iz-=2, ip++)
			{
			/* Copy last plane to current one */
			for (it=0; it<nz; it++)
				for (ix=0; ix<nx; ix++)
					p.vec[ip*p.nv*p.ne + it*p.ne + ix]
						= p.vec[(ip-1)*p.nv*p.ne + it*p.ne + ix];

			/* Pointer in plot to current z-level, in extra column */
			p.vec[(ip-1)*p.nv*p.ne + iz*p.ne + nx] = 5F;

			for (ix=0; ix<nx; ix++)
				{
				u[ix] = p.vec[ip*p.nv*p.ne + (iz-1)*p.ne + ix];
				v[ix] = u[ix];
				}

			for (it=iz; it<nt; it++)
				{
				/* Update the differencing star */
				for (ix=0; ix<nx; ix++)
					{
					w[ix] = u[ix];
					y[ix] = v[ix];
					v[ix] = p.vec[ip*p.nv*p.ne + it*p.ne + ix];
					}
				dd = (1-apb)*(v[0]+w[0]) + apb*(v[1]+w[1]);
				d[0] = dd - diag*y[0] - offdi*(y[0]+y[1]);
				for (ix=1; ix<nx-1; ix++)
					{
					dd = (1-2*apb)*(v[ix]+w[ix]);
					dd = dd + apb*(v[ix-1]+w[ix-1]+v[ix+1]+w[ix+1]);
					d[ix] = dd - diag*y[ix]-offdi*(y[ix-1]+y[ix+1]);
					}
				dd = (1-apb)*(v[nx-1]+w[nx-1])+apb*(v[nx-2]+w[nx-2]);
				d[nx-1] = dd - diag*y[nx-1] - offdi*(y[nx-1]+y[nx-2]);
				/* Solve the tridiagonal system */
				imhf.rtris(nx, diag+offdi, offdi, diag, offdi, diag+offdi,
					d, u, e, f);
				for (ix=0; ix<nx; ix++)
					p.vec[ip*p.nv*p.ne + it*p.ne + ix] = u[ix];
				}
			}
	
		/* Plot data volume in new frame, after setting label info */
		float rms = p.rms();
		p.setAmpUnits("Amplitude");
		p.setUnitsofDim(0, "Distance X");
		p.setDUnitofDim(0, 1F);
		p.setUnitsofDim(1, "Time-Depth Z");
		p.setDUnitofDim(1, 1F);
		p.setUnitsofDim(2, "=t");
		p.setDUnitofDim(2, 2F);
		p.viewinFrame("Lab 5: tdif.java", false, 0.5F,
			RgCtab.vagray, false, 3F*rms);
		
		/* All Done with tdif's main() application method */
		}

	/* End of methods in class tdif */
	}