/* 	Stolt migration methods without cosine weight.
	Adapted from Jon Claerbout's Imaging the Earth's Interior (1985), p. 71-72,
       by John N. Louie, 25 November 1998.
*/

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

/* Declare a class to contain the Stolt method and a main to test it */
/* When you change the code in this file, make a copy to a new file with
    a new name like "ex10stolttest.java", and be sure to make the same change
    to the class name below - e.g. "class ex10stolttest" */
class stolttest
	{
	/* The main() method compiles into a file named "stolttest.class" that
	    you run in the Java virtual machine as a stand-alone application
	    with a command such as "java stolttest"; or however you have
	    renamed it. */
	public static void main(String args[]) /* No runtime arguments needed */
		{
		System.out.println("Starting stolttest...");
		/* Declare variables */
		int it, nt, ix, nx;
		float vdtodx;
		/* Assign sizes */
		nx = 64;
		nt = 256;
		vdtodx = 1F/4F; /* vdtodx = v dt / dx */

		/* Declare and initialize planes to zero values */
		FltPlane p = new FltPlane(nx, nt);
		ComplexPlane cp = new ComplexPlane(nx, nt);
		
		/* Insert spikes into the real part of the complex array.
		    The class ComplexPlane has the method setElemReal(ix,it,val)
		    as a shortcut for the statement: cp.vec[ix*cp.ne + it].re = val; */
		cp.setElemReal(8,31,1F); cp.setElemReal(16,63,1F);
		cp.setElemReal(32,127,1F);

		/* Migrate */
		stolt(cp, vdtodx);

		/* Copy the real part of the transformed complex plane to a real plane
		   to allow plotting */
		cp.realtoFltPlane(p);
		
		/* Plot the real plane in a new window titled "stolttest", with
		    the fastest-varying direction pointing down, 0.25 vertical
		    exaggeration, using white color for negative amplitudes and
		    black color for positive amplitudes with near-zero amplitudes
		    gray (like looking at a WTVA wiggle-trace-variable-area plot from
		    a distance), scaling from the 3*rms clip down to minus the clip. */
		float rms = p.rms();
		p.viewinFrame("stolttest", true, 0.25F, RgCtab.vagray, false, 5F*rms);
		
		/* All Done with fttest's main() application method */
		}

	/* 	Stolt migration method without cosine weight.
	       From Claerbout (1985) p. 72. */
	static void stolt(ComplexPlane cp, float vdtodx)
		{
		/* Declare all local variables */
		int nt, nx;
		int ikx, nth, iktau, iom;
		float om, vkx, wl, wh, aktau, pionth, pi;
		int ix, it;
		
		nt = cp.getNMembofDim(0); /* rows of time */
		nx = cp.getNMembofDim(1);
		ComplexVec cbf = new ComplexVec(4*nt+1);
		/* if needed for plotting intermediate results */
		/* FltPlane p1 = new FltPlane(nx, nt); */

		/* Set some variables */
		pi = (float)(Math.PI);
		nth = nt/2;
		pionth = pi/nth;

		/* Transform data into the frequency-wavenumber domain
		    using your method compiled from fttestc.java */
		fttestc.ft2dc(cp, 1F, -1F);

		/* Interpolate the frequency axis onto the x-wavenumber axis */
		for (ikx=0; ikx<nx; ikx++)
			{
			vkx = ikx*2*pi*vdtodx/nx;
			if (ikx >= nx/2)
				/* Negative x-wavenumber */
				vkx = 2*pi*vdtodx - vkx;
			/* Set up scratch buffer cbf */
			cbf.vec[0].re = 0F; cbf.vec[0].im = 0F;
			cbf.vec[nt].re = 0F; cbf.vec[nt].im = 0F;
			for (iom=0; iom<nt; iom++)
				{
				/* Omit weighting */
				cbf.vec[iom].re = cp.vec[ikx*nt+iom].re;
				cbf.vec[iom].im = cp.vec[ikx*nt+iom].im;
				}
			/* Ignore zero freq*/
			cp.vec[ikx*nt+0].re = 0F; cp.vec[ikx*nt+0].im = 0F;
	
			/* Stretch */
			for (iktau=1; iktau<=nth; iktau++)
				{
				aktau = ((float)(iktau) - 0.01F)*pionth;
				/* From the paraxial wave equation -
				    argument to Math.sqrt must be positive */
				om = (float)(Math.sqrt((double)(aktau*aktau+vkx*vkx)));
				iom = (int)(om/pionth);
				if (iom < nth-1)
					{
					/* Linear interpolation ratios */
					wl = 1 + (float)(iom) - om/pionth;
					wh = 1 - wl;
					/* Mapping with interpolation between points */
					cp.vec[ikx*nt+iktau].re = cbf.vec[iom].re * wl
								+ cbf.vec[iom+1].re * wh;
					cp.vec[ikx*nt+iktau].im = cbf.vec[iom].im * wl
								+ cbf.vec[iom+1].im * wh;
					cp.vec[ikx*nt+nt-iktau].re = cbf.vec[nt-iom].re * wl
								+ cbf.vec[nt-iom-1].re * wh;
					cp.vec[ikx*nt+nt-iktau].im = cbf.vec[nt-iom].im * wl
								+ cbf.vec[nt-iom-1].im * wh;
					}
				else
					{
					cp.vec[ikx*nt+iktau].re = 0;
					cp.vec[ikx*nt+iktau].im = 0;
					}
				}
			}

		/* Transform back into the z-x domain */
		fttestc.ft2dc(cp, -1F, 1F);

		/* Return to calling method */
		}

	/* End of methods in class stolttest */
	}