///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Fourier.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 discrete Fourier transform
// with function f(x)=x(1-x) with x in [0,1].  The
// inverse transform is also performed for comparison.

import java.lang.*;
public class Fourier {
  static final int n = 128, m = 8;
  public static void main(String argv[]) {
  double x[] = new double[n];
  double fr[] = new double[n];
  double fi[] = new double[n];
  double gr[] = new double[n];
  double gi[] = new double[n];
  double h = 1.0/n;
  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;
    }
    dft(fr, fi, gr, gi);

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

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

// Method to perform the discrete Foruier transform.  Here
// fr[] and fi[] are the real and imaginary parts of the
// data with the corresponding outputs gr[] and gi[].
 
  public static void dft(double fr[], double fi[],
    double gr[], double gi[]) {
    int n = fr.length;
    double x = 2*Math.PI/n;
    for (int i=0; i<n; ++i) {
      for (int j=0; j<n; ++j) {
        double q = x*j*i;
        gr[i] += fr[j]*Math.cos(q)+fi[j]*Math.sin(q);
        gi[i] += fi[j]*Math.cos(q)-fr[j]*Math.sin(q);
      }
    }
  }
}