///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Wavelet.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 continuous wavelet
// transform with function f(t)=t(1-t) for 0<t<1.

import java.lang.*;
public class Wavelet{
  static final int nt = 21, nv = 101;
  public static void main(String argv[]) {
  double t[] = new double[nt];
  double f[] = new double[nt];
  double fi[] = new double[nt];
  double v[] = new double[nv];
  double g[][] = new double[nv][nv];
  double dt = 1.0/(nt-1), dv = 4.0/(nv-1);
  double rescale = 100;

 // Assign the data, dilate, and translate
    for (int i=0; i<nt; ++i) {
      t[i] = dt*i;
      f[i] = t[i]*(1-t[i]);
    }
    for (int i=0; i<nv; ++i) v[i] = dv*(i+1);

 // Perform the transform
    for (int i=0; i<nv; ++i){
      double sa = Math.sqrt(Math.abs(v[i]));
      for (int j=0; j<nv; ++j){
        for (int k=0; k<nt; ++k){
          fi[k]  = f[k]*(theta(t[k]-v[j])
                 - 2*theta(t[k]-v[j]-v[i]/2)
                 + theta(t[k]-v[j]-v[i]))/sa;
        }
        g[i][j] = simpson(fi,dt);
      }
    }

 // Output the coefficients obtained
    for (int i=0; i<nv; ++i){
      for (int j=0; j<nv; ++j){
        double g2 = rescale*g[i][j]*g[i][j];
        System.out.println(v[i] + " " + v[j] + " " + g2);
      }
    }
  }

// Method to create a step funciton.

  public static double theta(double x) {
    if (x>0) return 1.0;
    else return 0.0; 
  }

// Method to achieve the evenly spaced Simpson rule.

  public static double simpson(double y[], double h) {
    int n = y.length-1;
    double s0 = 0, s1 = 0, s2 = 0;
    for (int i=1; i<n; i+=2) {
      s0 += y[i];
      s1 += y[i-1];
      s2 += y[i+1];
    }
    double s = (s1+4*s0+s2)/3;

 // Add the last slice separately for an even n+1
    if ((n+1)%2 == 0)
      return h*(s+(5*y[n]+8*y[n-1]-y[n-2])/12);
    else
      return h*s;
  }
}