/* JRG, the Resource Geology Seismic Processing System for Java,
   and Viewmat for Java are Copyright 1998, 1999 by John N. Louie

The software and methods here are the subject of academic research,
not commercial products. I would like to know what use you make of
my methods, and have your feedback on their success or failure.
Also, by letting me know who you are, I can inform you if bugs or
errors are discovered.
 
Please send me an email message with: your name and email address;
whether you are a student or faculty member, consultant or employee;
the name of your university or company, and department; and a sentence
or two describing what use you intend to make of these methods.
Send your message to louie@seismo.unr.edu */


import java.util.*;
import java.io.*;
import java.lang.Math;

class testcomplexplane
	{
	public static void main(String args[])
		{
		int i;
		
		FltPlane fltplane = new FltPlane(210, 1000);
		fltplane.readBinaryFile("test.flt", 0);
		ComplexPlane cp = new ComplexPlane(fltplane, fltplane);
		
		cp.realtoFltPlane(fltplane);
		System.out.println("Real max=" + fltplane.max() + " rms=" + fltplane.rms()
			+ " dot=" + fltplane.dot(fltplane));
		cp.imtoFltPlane(fltplane);
		System.out.println("Imaginary max=" + fltplane.max() + " rms=" + fltplane.rms()
			+ " dot=" + fltplane.dot(fltplane));
		
		/*
		FltPlane fltplane = new FltPlane(210, 1000);
		fltplane.readBinaryFile("test.flt", 0);
		*/
		cp.powtoFltPlane(fltplane);
		float rms = fltplane.rms();
		fltplane.viewinFrame("test.flt.pow", true, (0.001F*2000F/2F)/15.24F,
			RgCtab.rwb, false, 3F*rms);
		
		/*
		FltPlane fp = new FltPlane(100, 178);
		fp.readBinaryFile("dixie-migvels", 0);	
		fp.addc(-(fp.min()));	
		fp.viewinFrame("dixie-migvels", false, 1F, RgCtab.yglobe, true, -1F);
		*/
		}
	}
	
public class ComplexPlane extends ComplexVec
	{
	public int ne;
	public int vecs;
	public String VecUnits;
	public float Vec0;
	public float DVec;

	public static void filt(ComplexPlane cdata, int ld, int lu,
		int hu, int hd, int spec)
		{
		int i;
		for (int i1=0; i1<cdata.vecs; i1++)
			{
			/*fftr(cdata, i1*cdata.ne, cdata.ne, -1F,
				(float)(Math.sqrt(1.0/(double)(cdata.ne))));*/
			fftr(cdata, i1*cdata.ne, cdata.ne, -1F, cdata.getDUnitofDim(0));
			if (spec == 1 || spec == 2)
				continue;
			for(i=0; i<ld ; i++)
				{       
				cdata.vec[i1*cdata.ne+i].re = 0F;
				cdata.vec[i1*cdata.ne+i].im = 0F;
				}
			for(i=ld; i<lu ; i++)
				{       
				cdata.vec[i1*cdata.ne+i].re =
					cdata.vec[i1*cdata.ne+i].re * (float)(i+1-ld)/(float)(lu-ld);
				cdata.vec[i1*cdata.ne+i].im =
					cdata.vec[i1*cdata.ne+i].im * (float)(i+1-ld)/(float)(lu-ld);
				}
			for(i=hu; i<hd ; i++)
				{       
				cdata.vec[i1*cdata.ne+i].re =
					cdata.vec[i1*cdata.ne+i].re * (float)(hd-i)/(float)(hd-hu);
				cdata.vec[i1*cdata.ne+i].im =
					cdata.vec[i1*cdata.ne+i].im * (float)(hd-i)/(float)(hd-hu);
				}
			for(i=hd; i<cdata.ne/2 ; i++)
				{       
				cdata.vec[i1*cdata.ne+i].re = 0F;
				cdata.vec[i1*cdata.ne+i].im = 0F;
				}
			cdata.vec[i1*cdata.ne+0].im=0.0F;
			cdata.vec[i1*cdata.ne+cdata.ne/2].re=0.0F;
			cdata.vec[i1*cdata.ne+cdata.ne/2].im=0.0F;
			for(i=1;i<cdata.ne/2;i++)
				{
				cdata.vec[i1*cdata.ne+cdata.ne-i].re
					= cdata.vec[i1*cdata.ne+i].re;
				cdata.vec[i1*cdata.ne+cdata.ne-i].im
					= (-cdata.vec[i1*cdata.ne+i].im);
				}
			/*fftr(cdata, i1*cdata.ne, cdata.ne, 1F,
				(float)(Math.sqrt(1.0/(double)(cdata.ne))));*/
			fftr(cdata, i1*cdata.ne, cdata.ne, 1F,
				1F/cdata.getDUnitofDim(0)/2F/cdata.ne/2F );
			}
		}
		
	/* 2-d FFT in-place method for class ComplexPlane data structures,
	    using Rocca's columnwise FFT scheme */
	static void ft2dc(ComplexPlane cp, float sign1, float sign2)
		{
		/* Declare local variables */
		int n1, n2;
		int i1, i2;
		float sc;
		
		/* Obtain n1=nx and n2=nt from the input ComplexPlane class
		    data structure, using RGData interface methods */
		n1 = cp.getNMembofDim(1); /* Dim=1 is vectors, or rows */
		n2 = cp.getNMembofDim(0); /* Dim=0 is fastest-varying, the elements */
		
		/* Adding String objects concatenates them */
		System.out.print("ft2dc: starting " + n2 + "x" + n1 + " 2-d FFT... ");
		
		/* Transform over the fast dimension, the rows */
		if (sign1 > 0F) /* time forward */
			sc = cp.getDUnitofDim(0);
		else /* time inverse */
			sc = 1F/cp.getDUnitofDim(0)/2F/n2/2F;
		/* Loop over n1, number of rows */
		for (i1=0; i1<n1; i1++)
			/* 1-d Fourier transform of row of plane cp, done n1 times.
			    Here use a form of the fft method that does not require
			    copying the row to a scratch vector, and back into cp */
			/*fftr(cp, i1*n2, n2, sign1, 1F);*/
			fftr(cp, i1*n2, n2, sign1, sc);
		
		/* Transform over the slow dimension, the columns; assumes
		    the complex planes are small enough to all fit in memory */
		if (sign1 < 0F) /* space forward */
			sc = cp.getDUnitofDim(1);
		else /* space inverse */
			sc = 1F/cp.getDUnitofDim(1)/2F/n1/2F;
		/*colmft(cp, sign2, 1F);*/
		colmft(cp, sign2, sc);
		
		/* Standard output terminated with a newline */
		System.out.println("Done.");
		/* Return to calling method */
		}

	/* Fabio Rocca's "row" Fourier transform, from Claerbout (1985) p. 73-74,
	    modified by J. Louie for the indexing scheme used by C (Dec. 1987)
	    and Java (Nov. 1998). */
	static void colmft(ComplexPlane cx, float sign2, float scale)
		{
		/* Declare all local variables */
		int i, i1, i2, j, m, lstep, istep, n1, n2;
		float arg;
		
		Complex cw = new Complex();
		Complex cdel = new Complex();
		/* Obtain n1=nx and n2=nt from the input ComplexPlane class
		    data structure, using RGData interface methods */
		n1 = cx.getNMembofDim(1); /* Dim=1 is vectors, or rows */
		n2 = cx.getNMembofDim(0); /* Dim=0 is fastest-varying, the elements */
		
		/* Multiply each element in matrix by scale factor */
		for (i1=0; i1<n1; i1++)
			for (i2=0; i2<n2; i2++)
				cx.vec[i1*n2+i2].scale(scale);
		
		j = 0;
		for (i=0; i<n1; i++)
			{
			if (i <= j)
				twid1(n2, cx, i*n2+0, cx, j*n2+0);
			m = n1/2;
			while (j > m-1)
				{
				j = j - m;
				m = m/2;
				if (m < 1)
					break;
				}
			j = j + m;
			}
		lstep = 1;
		do
			{
			istep = 2*lstep;
			cw.re = 1F;
			cw.im = 0F;
			arg = (float)(sign2*Math.PI/lstep);
			cdel.re = (float)(Math.cos((double)(arg)));
			cdel.im = (float)(Math.sin((double)(arg)));
			for (m=0; m<lstep; m++)
				{
				for (i=m; i<n1; i+=istep)
					twid2(n2, cw, cx, i*n2+0, cx, (i+lstep)*n2+0);
				cw.mul(cw, cdel);
				}
			lstep = istep;
			}
			while (lstep < n1);
		/* Return to calling method */
		}
	static void twid1(int n, ComplexPlane cx, int j0, ComplexPlane cy, int k0)
		{
		int i;
		Complex ct = new Complex();
		for (i=0; i<n; i++)
			{
			ct.re = cx.vec[j0+i].re;
			ct.im = cx.vec[j0+i].im;
			cx.vec[j0+i].re = cy.vec[k0+i].re;
			cx.vec[j0+i].im = cy.vec[k0+i].im;
			cy.vec[k0+i].re = ct.re;
			cy.vec[k0+i].im = ct.im;
			}
		}
	static void twid2(int n, Complex cw, ComplexPlane cx, int j0,
		ComplexPlane cy, int k0)
		{
		int i;
		Complex ctemp = new Complex();
		for (i=0; i<n; i++)
			{
			ctemp.mul(cw, cy.vec[k0+i]);
			cy.vec[k0+i].sub(cx.vec[j0+i], ctemp);
			cx.vec[j0+i].add(cx.vec[j0+i], ctemp);
			}
		}

	/* Simple in-place 1-d FFT method for a vector within a class ComplexPlane
	    data structure. The number of complex elements in each vector
	    MUST be a power of 2. j0 is the index of the first element of the
	    vector to process, and lx is the number of elements per vector. */
	static void fftr(ComplexPlane cx, int j0, int lx, float signi, float sc)
		{
		/* Declare all local variables */
		int i, j, k, m, istep;
		float arg;
		Complex cw = new Complex();
		Complex ct = new Complex();
		/* Here the lx argument supplies the vector length, since we are
		    working on one vector within the cx class ComplexPlane */
		
		j = 0;
		k = 1;
		/* Left out for compatibility with fttestc.colmft
		sc = (float)(Math.sqrt(1.0/(double)(lx)));
		*/
		for (i=0; i<lx; i++)
			{
			if (i <= j)
				{
				ct.re = sc * cx.vec[j0+j].re;
				ct.im = sc * cx.vec[j0+j].im;
				cx.vec[j0+j].re = sc * cx.vec[j0+i].re;
				cx.vec[j0+j].im = sc * cx.vec[j0+i].im;
				cx.vec[j0+i].re = ct.re;
				cx.vec[j0+i].im = ct.im;
				}
			m = lx/2;
			while (j > m-1)
				{
				j = j - m;
				m = m/2;
				if (m < 1)
					break;
				}
			j = j + m;
			}
		do
			{
			istep = 2*k;
			for (m=0; m<k; m++)
				{
				arg = (float)(Math.PI*signi*m/k);
				cw.re = (float)(Math.cos((double)(arg)));
				cw.im = (float)(Math.sin((double)(arg)));
				for (i=m; i<lx; i+=istep)
					{
					/* class Complex multiplication method ct = cw*cx[i+k] */
					ct.mul(cw, cx.vec[j0+i+k]);
					/* class Complex subtraction method cx[i+k] = cx[i] - ct */
					cx.vec[j0+i+k].sub(cx.vec[j0+i], ct);
					/* class Complex addition method cx[i] = cx[i] + ct */
					cx.vec[j0+i].add(cx.vec[j0+i], ct);
					}
				}
			k = istep;
			}
			while (k < lx);
		/* Return to calling method */
		}
	public static int getPowerof2(int N)
		{
		int pad2;
		pad2 = 1;
		while ( pad2 < N )
        		pad2 *= 2;
		return pad2;
		}

	public void setElem(int Vec, int Elem, Complex cpx)
		{
		this.vec[Vec*this.ne + Elem].re = cpx.re;
		this.vec[Vec*this.ne + Elem].im = cpx.im;
		}
	public void setElemReal(int Vec, int Elem, float Real)
		{
		this.vec[Vec*this.ne + Elem].re = Real;
		}
	public void setElemIm(int Vec, int Elem, float Im)
		{
		this.vec[Vec*this.ne + Elem].im = Im;
		}
	public Complex getElem(int Vec, int Elem)
		{
		return this.vec[Vec*this.ne + Elem];
		}
	public float getElemReal(int Vec, int Elem)
		{
		return this.vec[Vec*this.ne + Elem].re;
		}
	public float getElemIm(int Vec, int Elem)
		{
		return this.vec[Vec*this.ne + Elem].im;
		}
	public void realtoFltPlane(FltPlane fp)
		{
		int i;
		for (i=0; i<this.vec.length; i++)
			fp.vec[i] = this.vec[i].re;
		fp.NDims = this.NDims;
		fp.AmpUnits = new String("Real Part of " + this.AmpUnits);
		fp.ElemUnits = new String(this.ElemUnits);
		fp.Amp0 = this.Amp0;
		fp.AmpFac = this.AmpFac;
		fp.Elem0 = this.Elem0;
		fp.DElem = this.DElem;
		fp.ne = this.ne;
		fp.vecs = this.vecs;
		fp.VecUnits = new String(this.VecUnits);
		fp.Vec0 = this.Vec0;
		fp.DVec = this.DVec;
		}
	public void imtoFltPlane(FltPlane fp)
		{
		int i;
		for (i=0; i<this.vec.length; i++)
			fp.vec[i] = this.vec[i].im;
		fp.NDims = this.NDims;
		fp.AmpUnits = new String("Imaginary Part of " + this.AmpUnits);
		fp.ElemUnits = new String(this.ElemUnits);
		fp.Amp0 = this.Amp0;
		fp.AmpFac = this.AmpFac;
		fp.Elem0 = this.Elem0;
		fp.DElem = this.DElem;
		fp.ne = this.ne;
		fp.vecs = this.vecs;
		fp.VecUnits = new String(this.VecUnits);
		fp.Vec0 = this.Vec0;
		fp.DVec = this.DVec;
		}
	public void powtoFltPlane(FltPlane fp)
		{
		int i;
		for (i=0; i<this.vec.length; i++)
			fp.vec[i] = this.vec[i].pow();
		fp.NDims = this.NDims;
		fp.AmpUnits = new String("Power of " + this.AmpUnits);
		fp.ElemUnits = new String(this.ElemUnits);
		fp.Amp0 = this.Amp0;
		fp.AmpFac = this.AmpFac;
		fp.Elem0 = this.Elem0;
		fp.DElem = this.DElem;
		fp.ne = this.ne;
		fp.vecs = this.vecs;
		fp.VecUnits = new String(this.VecUnits);
		fp.Vec0 = this.Vec0;
		fp.DVec = this.DVec;
		}
	public void amptoFltPlane(FltPlane fp)
		{
		int i;
		for (i=0; i<this.vec.length; i++)
			fp.vec[i] = this.vec[i].amp();
		fp.NDims = this.NDims;
		fp.AmpUnits = new String("Magnitude of " + this.AmpUnits);
		fp.ElemUnits = new String(this.ElemUnits);
		fp.Amp0 = this.Amp0;
		fp.AmpFac = this.AmpFac;
		fp.Elem0 = this.Elem0;
		fp.DElem = this.DElem;
		fp.ne = this.ne;
		fp.vecs = this.vecs;
		fp.VecUnits = new String(this.VecUnits);
		fp.Vec0 = this.Vec0;
		fp.DVec = this.DVec;
		}
	public void phasetoFltPlane(FltPlane fp)
		{
		int i;
		for (i=0; i<this.vec.length; i++)
			fp.vec[i] = this.vec[i].phase();
		fp.NDims = this.NDims;
		fp.AmpUnits = new String("Phase, radians");
		fp.ElemUnits = new String(this.ElemUnits);
		fp.Amp0 = this.Amp0;
		fp.AmpFac = this.AmpFac;
		fp.Elem0 = this.Elem0;
		fp.DElem = this.DElem;
		fp.ne = this.ne;
		fp.vecs = this.vecs;
		fp.VecUnits = new String(this.VecUnits);
		fp.Vec0 = this.Vec0;
		fp.DVec = this.DVec;
		}
	public void phasedegtoFltPlane(FltPlane fp)
		{
		int i;
		for (i=0; i<this.vec.length; i++)
			fp.vec[i] = this.vec[i].phasedeg();
		fp.NDims = this.NDims;
		fp.AmpUnits = new String("Phase, degrees");
		fp.ElemUnits = new String(this.ElemUnits);
		fp.Amp0 = this.Amp0;
		fp.AmpFac = this.AmpFac;
		fp.Elem0 = this.Elem0;
		fp.DElem = this.DElem;
		fp.ne = this.ne;
		fp.vecs = this.vecs;
		fp.VecUnits = new String(this.VecUnits);
		fp.Vec0 = this.Vec0;
		fp.DVec = this.DVec;
		}
	public void scalebyFltPlane(FltPlane fp)
		{
		int i;
		for (i=0; i<this.vec.length; i++)
			this.vec[i].scale(fp.vec[i]);
		}
	public void scalebyFltVec(FltVec fv)
		{
		int i, j;
		for (j=0; j<this.vecs; j++)
			{
			for (i=0; i<this.ne; i++)
 				{
				this.vec[j*this.ne + i].scale(fv.vec[i]);
				}
			}
		}

	public ComplexPlane demult()
       		{
		int i, j;
		Notice("demult: demultiplexing " + this.vecs
			+ " vectors of " + this.ne + " complex elements to new plane with "
			+ this.ne + " vectors of " + this.vecs + " elements... ");
		ComplexPlane cp = new ComplexPlane(this.ne, this.vecs);
		cp.AmpUnits = new String(this.AmpUnits);
		cp.ElemUnits = new String(this.VecUnits);
		cp.VecUnits = new String(this.ElemUnits);
		cp.Amp0 = this.Amp0;
		cp.AmpFac = this.AmpFac;
		cp.Elem0 = this.Vec0;
		cp.DElem = this.DVec;
		cp.Vec0 = this.Elem0;
		cp.DVec = this.DElem;
		for (j=0; j<cp.vecs; j++)
			{
			for (i=0; i<cp.ne; i++)
 				{
				cp.vec[j*cp.ne + i].re = this.vec[i*this.ne + j].re;
				cp.vec[j*cp.ne + i].im = this.vec[i*this.ne + j].im;
				}
			}
		Noticeln("done.");
		return cp;
		}
	
	public void windowof(ComplexPlane cp, int v0, int dv, int e0, int de)
		{
		int i, j;
		Notice("windowof: windowing " + this.ne + " elements and "
			+ this.vecs + " vectors from complex plane, with v0="
			+ v0 + " dv=" + dv + " e0=" + e0 + " de=" + de + "... ");
		for (j=0; j<this.vecs; j++)
			for (i=0; i<this.ne; i++)
				{
				this.vec[j*this.ne+i].re =
					cp.vec[(cp.ne*(v0+(j*dv))) + (e0+(i*de))].re;
				this.vec[j*this.ne+i].im =
					cp.vec[(cp.ne*(v0+(j*dv))) + (e0+(i*de))].im;
				}
		this.AmpUnits = new String(cp.AmpUnits);
		this.Amp0 = cp.Amp0;
		this.AmpFac = cp.AmpFac;
		this.ElemUnits = new String(cp.ElemUnits);
		this.Elem0 = cp.Elem0 + e0;
		this.DElem = cp.DElem*de;
		this.VecUnits = new String(cp.VecUnits);
		this.Vec0 = cp.Vec0 + v0;
		this.DVec = cp.DVec*dv;
		Noticeln("done.");
		}

       /* Constructors */
	ComplexPlane(int NVec, int NElem)
		{
		super(NVec*NElem);
		this.NDims = 2;
		this.vecs = NVec;
		this.ne = NElem;
		this.VecUnits = new String("Vector");
		this.Vec0 = 0F;
		this.DVec = 1F;
		}
	/* Construct a copy */
	ComplexPlane(ComplexPlane cp) 
		{
		super(cp.vec.length);
		int i;
		for (i=0; i<cp.vec.length; i++)
			{
			this.vec[i].re = cp.vec[i].re;
			this.vec[i].im = cp.vec[i].im;
			}
		this.NDims = cp.NDims;
		this.vecs = cp.vecs;
		this.ne = cp.ne;
		this.AmpUnits = new String(cp.AmpUnits);
		this.Amp0 = cp.Amp0;
		this.AmpFac = cp.AmpFac;
		this.ElemUnits = new String(cp.ElemUnits);
		this.Elem0 = cp.Elem0;
		this.DElem = cp.DElem;
		this.VecUnits = new String(cp.VecUnits);
		this.Vec0 = cp.Vec0;
		this.DVec = cp.DVec;
		}
	/* Construct a copy of a real plane */
	ComplexPlane(FltPlane fp) 
		{
		super(fp.vec.length);
		int i;
		for (i=0; i<fp.vec.length; i++)
			{
			this.vec[i].re = fp.vec[i];
			this.vec[i].im = 0F;
			}
		this.NDims = fp.NDims;
		this.vecs = fp.vecs;
		this.ne = fp.ne;
		this.AmpUnits = new String(fp.AmpUnits);
		this.Amp0 = fp.Amp0;
		this.AmpFac = fp.AmpFac;
		this.ElemUnits = new String(fp.ElemUnits);
		this.Elem0 = fp.Elem0;
		this.DElem = fp.DElem;
		this.VecUnits = new String(fp.VecUnits);
		this.Vec0 = fp.Vec0;
		this.DVec = fp.DVec;
		}
	/* Construct a copy from separate real and imaginary planes */
	ComplexPlane(FltPlane rp, FltPlane ip) 
		{
		super(rp.vec.length);
		int i;
		for (i=0; i<rp.vec.length; i++)
			{
			this.vec[i].re = rp.vec[i];
			this.vec[i].im = ip.vec[i];
			}
		this.NDims = rp.NDims;
		this.vecs = rp.vecs;
		this.ne = rp.ne;
		this.AmpUnits = new String(rp.AmpUnits);
		this.Amp0 = rp.Amp0;
		this.AmpFac = rp.AmpFac;
		this.ElemUnits = new String(rp.ElemUnits);
		this.Elem0 = rp.Elem0;
		this.DElem = rp.DElem;
		this.VecUnits = new String(rp.VecUnits);
		this.Vec0 = rp.Vec0;
		this.DVec = rp.DVec;
		}
		
	/* RGData interface methods */
	public int getNMembofDim(int Dim)
		{
		if (Dim == 0)
			return(ne);
		else if (Dim == 1)
			return(vecs);
		else if (Dim >= NDims)
			return(1);
		else
			return(0);
		}
	public String getUnitsofDim(int Dim)
		{
		if (Dim == 0)
			return(ElemUnits);
		else if (Dim == 1)
			return(VecUnits);
		return("");
		}
	public void setUnitsofDim(int Dim, String Dimunits)
		{
		if (Dim == 0)
			ElemUnits = new String(Dimunits);
		else if (Dim == 1)
			VecUnits = new String(Dimunits);
		}
	public float getUnit0ofDim(int Dim)
		{
		if (Dim == 0)
			return Elem0;
		else if (Dim == 1)
			return Vec0;
		return 0F;
		}
	public void setUnit0ofDim(int Dim, float Unit0)
		{
		if (Dim == 0)
			Elem0 = Unit0;
		else if (Dim == 1)
			Vec0 = Unit0;
		}
	public float getDUnitofDim(int Dim)
		{
		if (Dim == 0)
			return DElem;
		else if (Dim == 1)
			return DVec;
		return 1F;
		}
	public void setDUnitofDim(int Dim, float DUnit)
		{
		if (Dim == 0)
			DElem = DUnit;
		else if (Dim == 1)
			DVec = DUnit;
		}
	}