/* 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 */


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

class testcomplexvec
	{
	public static void main(String args[])
		{
		int i;
		System.out.println("Starting construction and reading.");
		FltVec fltvec = new FltVec(210*1000);
		
		/* read test */
		fltvec.readBinaryFile("test.flt", 0);
		ComplexVec cvec = new ComplexVec(fltvec, fltvec);
		
		cvec.realtoFltVec(fltvec);
		System.out.println("Real max=" + fltvec.max() + " rms=" + fltvec.rms()
			+ " dot=" + fltvec.dot(fltvec));
		cvec.imtoFltVec(fltvec);
		System.out.println("Imaginary max=" + fltvec.max() + " rms=" + fltvec.rms()
			+ " dot=" + fltvec.dot(fltvec));
		
		cvec.amptoFltVec(fltvec);
		System.out.println("Amp max=" + fltvec.max() + " rms=" + fltvec.rms()
			+ " dot=" + fltvec.dot(fltvec));
		cvec.phasedegtoFltVec(fltvec);
		System.out.println("Phase max=" + fltvec.max() + " rms=" + fltvec.rms()
			+ " dot=" + fltvec.dot(fltvec));
		}
	}


public class ComplexVec implements RGData, RGNotice
	{
	public int NDims;
	public String AmpUnits;
	public float Amp0;
	public float AmpFac;
	public String ElemUnits;
	public float Elem0;
	public float DElem;
	public boolean noticeon;
	public Complex vec[];
	
	public static void filt(ComplexVec cdata, int ld, int lu,
		int hu, int hd, int spec)
		{
		int i;
		/*fft(cdata, -1F, (float)(Math.sqrt(1.0/(double)(cdata.vec.length))));*/
		fft(cdata, -1F, cdata.getDUnitofDim(0));
		if (spec == 1 || spec == 2)
			return;
		for(i=0; i<ld ; i++)
			{       
			cdata.vec[i].re = 0F;
			cdata.vec[i].im = 0F;
			}
		for(i=ld; i<lu ; i++)
			{       
			cdata.vec[i].re =
				cdata.vec[i].re * (float)(i+1-ld)/(float)(lu-ld);
			cdata.vec[i].im =
				cdata.vec[i].im * (float)(i+1-ld)/(float)(lu-ld);
			}
		for(i=hu; i<hd ; i++)
			{       
			cdata.vec[i].re =
				cdata.vec[i].re * (float)(hd-i)/(float)(hd-hu);
			cdata.vec[i].im =
				cdata.vec[i].im * (float)(hd-i)/(float)(hd-hu);
			}
		for(i=hd; i<cdata.vec.length/2 ; i++)
			{       
			cdata.vec[i].re = 0F;
			cdata.vec[i].im = 0F;
			}
		cdata.vec[0].im = 0.0F;
		cdata.vec[cdata.vec.length/2].re = 0.0F;
		cdata.vec[cdata.vec.length/2].im = 0.0F;
		for(i=1; i<cdata.vec.length/2; i++)
			{
			cdata.vec[cdata.vec.length-i].re
				= cdata.vec[i].re;
			cdata.vec[cdata.vec.length-i].im
				= (-cdata.vec[i].im);
			}
		/*fft(cdata, 1F, (float)(Math.sqrt(1.0/(double)(cdata.vec.length))));*/
		fft(cdata, 1F, 1F/cdata.getDUnitofDim(0)/2F/cdata.vec.length/2F);
		}
		
	/* Simple 1-d FFT method for class ComplexVec data structures.
	    The number of complex elements in the vector MUST be a power of 2 */
	static void fft(ComplexVec cx, float signi, float sc)
		{
		/* Declare all local variables */
		int i, j, k, m, istep, lx;
		float arg;
		Complex cw = new Complex();
		Complex ct = new Complex();
		/* Obtain Complex vector length lx from the input class ComplexVec
		    data structure, using RGData interface method */
		lx = cx.getNMembofDim(0); /* Dim=0 is fastest-varying, the elements */
		
		j = 0;
		k = 1;
		for (i=0; i<lx; i++)
			{
			if (i <= j)
				{
				ct.re = sc * cx.vec[j].re;
				ct.im = sc * cx.vec[j].im;
				cx.vec[j].re = sc * cx.vec[i].re;
				cx.vec[j].im = sc * cx.vec[i].im;
				cx.vec[i].re = ct.re;
				cx.vec[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[i+k]);
					/* class Complex subtraction method cx[i+k] = cx[i] - ct */
					cx.vec[i+k].sub(cx.vec[i], ct);
					/* class Complex addition method cx[i] = cx[i] + ct */
					cx.vec[i].add(cx.vec[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 Elem, Complex cpx)
		{
		this.vec[Elem].re = cpx.re;
		this.vec[Elem].im = cpx.im;
		}
	public void setElemReal(int Elem, float Real)
		{
		this.vec[Elem].re = Real;
		}
	public void setElemIm(int Elem, float Im)
		{
		this.vec[Elem].im = Im;
		}
	public Complex getElem(int Elem)
		{
		return this.vec[Elem];
		}
	public float getElemReal(int Elem)
		{
		return this.vec[Elem].re;
		}
	public float getElemIm(int Elem)
		{
		return this.vec[Elem].im;
		}
	public void realtoFltVec(FltVec fv)
		{
		int i;
		for (i=0; i<this.vec.length; i++)
			fv.vec[i] = this.vec[i].re;
		fv.NDims = this.NDims;
		fv.AmpUnits = new String("Real Part of " + this.AmpUnits);
		fv.ElemUnits = new String(this.ElemUnits);
		fv.Amp0 = this.Amp0;
		fv.AmpFac = this.AmpFac;
		fv.Elem0 = this.Elem0;
		fv.DElem = this.DElem;
		}
	public void imtoFltVec(FltVec fv)
		{
		int i;
		for (i=0; i<this.vec.length; i++)
			fv.vec[i] = this.vec[i].im;
		fv.NDims = this.NDims;
		fv.AmpUnits = new String("Imaginary Part of " + this.AmpUnits);
		fv.ElemUnits = new String(this.ElemUnits);
		fv.Amp0 = this.Amp0;
		fv.AmpFac = this.AmpFac;
		fv.Elem0 = this.Elem0;
		fv.DElem = this.DElem;
		}
	public void powtoFltVec(FltVec fv)
		{
		int i;
		for (i=0; i<this.vec.length; i++)
			fv.vec[i] = this.vec[i].pow();
		fv.NDims = this.NDims;
		fv.AmpUnits = new String("Power of " + this.AmpUnits);
		fv.ElemUnits = new String(this.ElemUnits);
		fv.Amp0 = this.Amp0;
		fv.AmpFac = this.AmpFac;
		fv.Elem0 = this.Elem0;
		fv.DElem = this.DElem;
		}
	public void amptoFltVec(FltVec fv)
		{
		int i;
		for (i=0; i<this.vec.length; i++)
			fv.vec[i] = this.vec[i].amp();
		fv.NDims = this.NDims;
		fv.AmpUnits = new String("Magnitude of " + this.AmpUnits);
		fv.ElemUnits = new String(this.ElemUnits);
		fv.Amp0 = this.Amp0;
		fv.AmpFac = this.AmpFac;
		fv.Elem0 = this.Elem0;
		fv.DElem = this.DElem;
		}
	public void phasetoFltVec(FltVec fv)
		{
		int i;
		for (i=0; i<this.vec.length; i++)
			fv.vec[i] = this.vec[i].phase();
		fv.NDims = this.NDims;
		fv.AmpUnits = new String("Phase, radians");
		fv.ElemUnits = new String(this.ElemUnits);
		fv.Amp0 = this.Amp0;
		fv.AmpFac = this.AmpFac;
		fv.Elem0 = this.Elem0;
		fv.DElem = this.DElem;
		}
	public void phasedegtoFltVec(FltVec fv)
		{
		int i;
		for (i=0; i<this.vec.length; i++)
			fv.vec[i] = this.vec[i].phasedeg();
		fv.NDims = this.NDims;
		fv.AmpUnits = new String("Phase, degrees");
		fv.ElemUnits = new String(this.ElemUnits);
		fv.Amp0 = this.Amp0;
		fv.AmpFac = this.AmpFac;
		fv.Elem0 = this.Elem0;
		fv.DElem = this.DElem;
		}
	public void scalebyFltVec(FltVec fv)
		{
		int i;
		for (i=0; i<this.vec.length; i++)
			this.vec[i].scale(fv.vec[i]);
		}
	public void windowof(ComplexVec cv, int e0, int de)
		{
		int i;
		Notice("windowof: windowing " + vec.length
			+ " elements from complex array of " + cv.vec.length
			+ ", e0=" + e0 + " de=" + de + "... ");
		for (i=0; i<this.vec.length; i++)
			{
			this.vec[i].re = cv.vec[e0+(i*de)].re;
			this.vec[i].im = cv.vec[e0+(i*de)].im;
			}
		this.AmpUnits = new String(cv.AmpUnits);
		this.Amp0 = cv.Amp0;
		this.AmpFac = cv.AmpFac;
		this.ElemUnits = new String(cv.ElemUnits);
		this.Elem0 = cv.Elem0 + e0;
		this.DElem = cv.DElem*de;
		Noticeln("done.");
		}
	public void windowof(ComplexPlane cp, int v0, int e0, int de)
		{
		int i;
		Notice("windowof: windowing " + vec.length
			+ " elements from vector " + v0 + ", of "
			+ cp.ne + "elements, from complex plane with e0=" + e0
			+ " de=" + de + "... ");
		for (i=0; i<this.vec.length; i++)
			{
			this.vec[i].re = cp.vec[(cp.ne*v0) + (e0+(i*de))].re;
			this.vec[i].im = cp.vec[(cp.ne*v0) + (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;
		Noticeln("done.");
		}
	public void tileto(ComplexPlane cp, int v0, int e0, int de)
		{
		int i;
		Notice("windowof: copying " + vec.length
			+ " elements into vector " + v0 + ", of "
			+ cp.ne + "elements, to complex plane with e0=" + e0
			+ " de=" + de + "... ");
		for (i=0; i<this.vec.length; i++)
			{
			cp.vec[(cp.ne*v0) + (e0+(i*de))].re = this.vec[i].re;
			cp.vec[(cp.ne*v0) + (e0+(i*de))].im = this.vec[i].im;
			}
		Noticeln("done.");
		}

	/* Constructors */
	ComplexVec(int NElem)
		{
		int i;
		vec = new Complex[NElem];
		for (i=0; i<NElem; i++)
			vec[i] = new Complex();
		this.NDims = 1;
		this.AmpUnits = new String("Complex Amplitude");
		this.ElemUnits = new String("Element");
		this.Amp0 = 0F;
		this.AmpFac = 1F;
		this.Elem0 = 0F;
		this.DElem = 1F;
		this.noticeon = true;
		}
	/* Construct a copy */
	ComplexVec(ComplexVec cv) 
		{
		int i;
		this.vec = new Complex[cv.vec.length];
		for (i=0; i<cv.vec.length; i++)
			this.vec[i] = new Complex(cv.vec[i]);
		this.NDims = cv.NDims;
		this.AmpUnits = new String(cv.AmpUnits);
		this.ElemUnits = new String(cv.ElemUnits);
		this.Amp0 = cv.Amp0;
		this.AmpFac = cv.AmpFac;
		this.Elem0 = cv.Elem0;
		this.DElem = cv.DElem;
		this.noticeon = cv.noticeon;
		}
	/* Construct a copy of a real vector */
	ComplexVec(FltVec fv) 
		{
		int i;
		this.vec = new Complex[fv.vec.length];
		for (i=0; i<fv.vec.length; i++)
			this.vec[i] = new Complex(fv.vec[i]);
		this.NDims = fv.NDims;
		this.AmpUnits = new String(fv.AmpUnits);
		this.ElemUnits = new String(fv.ElemUnits);
		this.Amp0 = fv.Amp0;
		this.AmpFac = fv.AmpFac;
		this.Elem0 = fv.Elem0;
		this.DElem = fv.DElem;
		this.noticeon = fv.noticeon;
		}
	/* Construct a copy from separate real and imaginary vectors */
	ComplexVec(FltVec rv, FltVec iv) 
		{
		int i;
		this.vec = new Complex[rv.vec.length];
		for (i=0; i<rv.vec.length; i++)
			this.vec[i] = new Complex(rv.vec[i], iv.vec[i]);
		this.NDims = rv.NDims;
		this.AmpUnits = new String(rv.AmpUnits);
		this.ElemUnits = new String(rv.ElemUnits);
		this.Amp0 = rv.Amp0;
		this.AmpFac = rv.AmpFac;
		this.Elem0 = rv.Elem0;
		this.DElem = rv.DElem;
		this.noticeon = rv.noticeon;
		}
	
	/* RGData interface methods */
	public int getNDims()
		{ return NDims; }
	public int getNMembofDim(int Dim)
		{
		if (Dim == 0)
			return vec.length;
		else if (Dim >= NDims)
			return 1;
		else
			return 0;
		}
	public String getAmpUnits()
		{
		return AmpUnits;		
		}
	public void setAmpUnits(String Ampunits)
		{
		AmpUnits = new String(Ampunits);
		}
	public float getAmp0()
		{
		return Amp0;		
		}
	public void setAmp0(float Amp0)
		{
		this.Amp0 = Amp0;
		}
	public float getAmpFac()
		{
		return AmpFac;		
		}
	public void setAmpFac(float Ampfac)
		{
		this.AmpFac = Ampfac;
		}
	public String getUnitsofDim(int Dim)
		{
		if (Dim == 0)
			return ElemUnits;
		return "";
		}
	public void setUnitsofDim(int Dim, String Dimunits)
		{
		if (Dim == 0)
			ElemUnits = new String(Dimunits);
		}
	public float getUnit0ofDim(int Dim)
		{
		if (Dim == 0)
			return Elem0;
		return 0F;
		}
	public void setUnit0ofDim(int Dim, float Unit0)
		{
		if (Dim == 0)
			Elem0 = Unit0;
		}
	public float getDUnitofDim(int Dim)
		{
		if (Dim == 0)
			return DElem;
		return 1F;
		}
	public void setDUnitofDim(int Dim, float DUnit)
		{
		if (Dim == 0)
			DElem = DUnit;
		}
	public void Errout(String Error)
		{
		System.out.println(Error);
		System.out.flush();
		}
	public void Noticeln(String Note)
		{
		if (noticeon)
			{
			System.out.println(Note);
			System.out.flush();
			}
		}
	public void Notice(String Note)
		{
		if (noticeon)
			{
			System.out.print(Note);
			System.out.flush();
			}
		}
	public void setNoticeOn(boolean on)
		{
		this.noticeon = on;
		}
	}