///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Fourier2.java                                      //
//                                                                       //
// © Tao Pang 2006                                                       //
//                                                                       //
// Last modified: January 18, 2006                                       //
//                                                                       //
// (1) This Java program is part of the book, "An Introduction to        //
//     Computational Physics, 2nd Edition," written by Tao Pang and      //
//     published by Cambridge University Press on January 19, 2006.      //
//                                                                       //
// (2) No warranties, express or implied, are made for this program.     //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

// An example of performing the fast Fourier transform
// with function f(x) = x(1-x) in [0,1].  The inverse
// transform is also performed for comparison.

import java.lang.*;
public class Fourier2 {
  static final int n = 128, l = 8;
  public static void main(String argv[]) {
    double x[] = new double[n];
    double fr[] = new double[n];
    double fi[] = new double[n];
    int m = 7;
    double h = 1.0/(n-1);
    double f0 = 1/Math.sqrt(n);

 // Assign the data and perform the transform
    for (int i=0; i<n; ++i) {
      x[i] = h*i;
      fr[i] = x[i]*(1-x[i]);
      fi[i] = 0;
    }
    fft(fr, fi, m);

 // Perform the inverse Fourier transform
    for (int i=0; i<n; ++i) {
      fr[i] =  f0*fr[i];
      fi[i] = -f0*fi[i];
    }
    fft(fr, fi, m);
    for (int i=0; i<n; ++i) {
      fr[i] =  f0*fr[i];
      fi[i] = -f0*fi[i];
    }

 // Output the result in every l data steps
    for (int i=0; i<n; i+=l)
      System.out.println(x[i] + " " + fr[i] + " " + fi[i]);
  }

// Method to perform the fast Foruier transform.  Here
// fr[] and fi[] are the real and imaginary parts of 
// both the input and output data. 

  public static void fft(double fr[], double fi[],
    int m) {
    int n = fr.length;
    int nh = n/2;
    int np = (int) Math.pow(2, m);

 // Stop the program if the indices do not match
    if (np != n) {
      System.out.println("Index mismtch detected");
      System.exit(1);
    }

 // Rearrange the data to the bit-reversed order
    int k = 1;
    for (int i=0; i<n-1; ++i) {
      if (i < k-1) {
        double f1 = fr[k-1];
        double f2 = fi[k-1];
        fr[k-1] = fr[i];
        fi[k-1] = fi[i];
        fr[i] = f1;
        fi[i] = f2;
      }
      int j = nh;
      while (j < k) {
        k -= j;
        j /= 2;
      }
      k += j;
    }

 // Sum up the reordered data at all levels
    k = 1;
    for (int i=0; i<m; ++i) {
      double w = 0;
      int j = k;
      k = 2*j;
      for (int p=0; p<j; ++p) {
        double u =  Math.cos(w);
        double v = -Math.sin(w);
        w += Math.PI/j;
        for (int q=p; q<n; q+=k) {
          int r = q+j;
          double f1 = fr[r]*u-fi[r]*v;
          double f2 = fr[r]*v+fi[r]*u;
          fr[r] = fr[q]-f1;
          fr[q] += f1;
          fi[r] = fi[q]-f2;
          fi[q] += f2;
        }
      }
    }
  }
}