/*  This file is part of the source code for 3D-XplorMath-J, Version 1.0 (January 2008).
 *  Copyright (c) 2008 The 3D-XplorMath Consortium (http://3d-xplormath.org).
 *  This source code is released under a BSD License, which allows redistribution   
 *  in source and binary form, with or without modification, provided copyright
 *  and license information are included, and with no warranty or guarantee of
 *  any kind.  For details, see http://3d-xplormath.org/j/source/BSDLicense.txt
 */
 
 package vmm3d.fourier;

import vmm3d.core.Complex;

public class FastFourier {
	static final int MAXFFTSIZE = 512;

	static Complex twid[]; 
	static Complex Cadd(Complex a,Complex b) {return a.plus(b);}
	static Complex Csub(Complex a,Complex b) {return a.minus(b);}
    static Complex Cmul(Complex a, Complex b) {return a.times(b);}
	static Complex Cdiv(Complex a,Complex b) {return a.dividedBy(b);}
	static Complex RCmul(double a, Complex b) {return b.times(a);}
	static Complex Cexp(Complex a) {return a.exponential();}
	static Complex Comp(double a,double b){Complex c = Complex.ZERO_C; c.assign(a,b);return c;}	
	static Complex cexp(double x,double y) { return RCmul(Math.exp(x),Comp(Math.cos(y),Math.sin(y)));}
	static Complex Cconjg(Complex c) { return c.conj();}
	static double Cabs(Complex c) {return c.r();}
	static Complex Cexpi(double y) {return Comp(Math.cos(y),Math.sin(y));}

	private static void SetupTwid() {
		if (twid == null) {
			 twid = new Complex[MAXFFTSIZE];
	         double h = 2*Math.PI/MAXFFTSIZE;
	         for (int j = 0; j< MAXFFTSIZE;j++) {
	    	         twid[j] = new Complex();
	                 twid[j].re = Math.cos(j*h);
	                 twid[j].im = Math.sin(j*h);
	      }
	   }
	}
	
	static void fft( Complex[] c, int J) {
		int i,j,bit;
		int numtrans,sizetrans,halfsize;
		int level,in,sub;
		int iminus,iplus;
		Complex twiddle;
		Complex but1,but2;
		Complex temp;
		SetupTwid();
		/* bit reversal */
		i=0;
		for (j=1; j< J; j++) { /* count with j */
		/* binary add 1 to i in mirror */
		bit=J/2;
		while ( i >= bit ) {  /* until you encounter 0, change 1 to 0 */
		i=i-bit;
		bit=bit/2;}
		i=i+bit;		/* then change 0 to 1 */
		if (i 1)  {
		numtrans=numtrans/2;     /* at each level, do half as many */
		halfsize=sizetrans;
		sizetrans=sizetrans*2;   /* subtransforms of twice the size   */
		for (in=0; in