// test_fourier.java 
//                                                             
// y_approx(xj) = a0/2 + sum_i=j,n( ai*cos(i*2*Pi*xj/T)+       
//                                  bi*sin(i*2*Pi*xj/T) )      
//            ai = 2/T * integral_0,T( f(x)*cos(i*2*Pi*xj/T)dx 
//            bi = 2/T * integral_0,T( f(x)*sin(i*2*Pi*xj/T)dx 
//                                                             
// for numerical computation, x is scaled to be in range [0,1] 
// xs = (x-xmin)/(xmax-xmin)  T=1                              
// then y_approx is computed  y_approx((x-xmin)/(xmax-xmin))   
//                                                             
// a "smoother" fit can be obtained by a Fejer approximation.  
// many coefficients may be near zero.                         

public class test_fourier
{
  double f(double x) // simple step function 0 to 2pi
                     // poor for odd number of points
  {
    if(x<Math.PI) return -0.5;
    return 0.5;
  }

  public test_fourier()
  {  
    boolean debug = true;
    double Pi = 3.1415926535897932384626433832795028841971;
    int n = 101;
    double x[] = new double[101]; // 0.0 to 2Pi
    double y[] = new double[101]; // 0.0 to pi
    double yf[] = new double[101]; // fit
    final double tpi = 2.0*Math.PI;
    final double dx1 = tpi/100.0;
    int nf = 10; // could be nt  number of terms
    double a[] = new double[nf];
    double b[] = new double[nf];
    double xmin, xmax, ymin, ymax, yavg;
    double xs, xs1, xs2, ti, dx;
    double err, avgerr=0.0, rmserr=0.0, maxerr=0.0;
    double sum, sumsq=0.0;
    double max_err = 0.0;
    int nvert, npoly, ipoly;
    double vert[][] = new double[n+6][3]; // for datwrite
    int poly[] = new int[3*n+3*5];        // for datwrite
    
    System.out.println("test_fourier.java running");
    System.out.println("data to fit x,y  find yf=f(x)");
    for(int i=0; i<n; i++)
    {
      x[i] = (double)i*dx1;
      y[i]=f(x[i]);
      if(debug) System.out.println("x["+i+"]="+x[i]+", y["+i+"]="+y[i]);    
    } // end i
    System.out.println(" ");
    xmin = x[0]; // generally not known
    xmax = x[0];
    ymin = y[0];
    ymax = y[0];
    for(int j=1; j<n; j++)
    {
      xmin = Math.min(x[j],xmin);
      xmax = Math.max(x[j],xmax);
      ymin = Math.min(y[j],ymin);
      ymax = Math.max(y[j],ymax);
    }
    yavg = (ymin+ymax)/2.0;
    System.out.println("xmin="+xmin+", xmax="+xmax+
                       ", ymin="+ymin+", ymax="+ymax+
		       ", yavg="+yavg);

    // integration, simple trapezoidal rule 
    System.out.println("coefficients");
    for(int i=0; i<nf; i++)
    {
      a[i] = 0.0;
      b[i] = 0.0;
      ti = (double)i;
      for(int j=0; j<n-1; j++)
      {
        xs1 = (x[j]-xmin)/(xmax-xmin);
        xs2 = (x[j+1]-xmin)/(xmax-xmin);
        dx = xs2-xs1;
        a[i] = a[i] + dx*(y[j]*Math.cos(ti*tpi*xs1)+y[j+1]*
                               Math.cos(ti*tpi*xs2));
        b[i] = b[i] + dx*(y[j]*Math.sin(ti*tpi*xs1)+y[j+1]*
                               Math.sin(ti*tpi*xs2));
      }
      System.out.println("a["+i+"]="+a[i]+", b["+i+"]="+b[i]);
    }
    System.out.println(" ");
    // check fit 

    for(int j=0; j<n; j++)
    {
      sum = 0.0;
      sumsq = 0.0;
      xs = (x[j]-xmin)/(xmax-xmin);
      for(int i=0; i<nf; i++)
      {
        ti = (double)i;
        sum = sum + a[i]*Math.cos(ti*tpi*xs) + b[i]*Math.sin(ti*tpi*xs);
      }
      yf[j] = sum-yavg;
      err = y[j]-yf[j];
      sumsq = sumsq + err*err;
      avgerr = avgerr + Math.abs(err);
      maxerr = Math.max(Math.abs(err),maxerr);
      if(debug) System.out.println("y["+j+"]="+y[j]+", yf["+j+"]="+yf[j]+
                                   ", err"+err);
    }
    rmserr = Math.sqrt(sumsq/(double)n);
    avgerr = avgerr/(double)n;
    System.out.println("avgerr="+avgerr+", rmserr="+rmserr+", maxerr="+maxerr);
    System.out.println("test_fourier.java finished ");
    
    nvert = 0;
    npoly = 0;
    ipoly = 0;
    vert[nvert][0] = 0.0;   // start to point
    vert[nvert][1] = 0.0;
    vert[nvert][2] = 0.0;
    nvert = nvert+1;
    vert[nvert][0] = 0.0;   // start to point
    vert[nvert][1] = -0.5;
    vert[nvert][2] = 0.0;
    nvert = nvert+1;
    vert[nvert][0] = Math.PI;
    vert[nvert][1] = -0.5;
    vert[nvert][2] = 0.0;
    nvert = nvert+1;
    vert[nvert][0] = Math.PI; // point to end 
    vert[nvert][1] = 0.5;
    vert[nvert][2] = 0.0;
    nvert = nvert+1;
    vert[nvert][0] = 2.0*Math.PI;
    vert[nvert][1] = 0.5;
    vert[nvert][2] = 0.0;
    nvert = nvert+1;
    vert[nvert][0] = 2.0*Math.PI;
    vert[nvert][1] = 0.0;
    vert[nvert][2] = 0.0;
    nvert = nvert+1;
    poly[ipoly] = 2; // start to point
    ipoly = ipoly+1;
    poly[ipoly] = 1;
    ipoly = ipoly+1;
    poly[ipoly] = 2;
    ipoly = ipoly+1;
    npoly = npoly+1;
    poly[ipoly] = 2;
    ipoly = ipoly+1;
    poly[ipoly] = 2;
    ipoly = ipoly+1;
    poly[ipoly] = 3;
    ipoly = ipoly+1;
    npoly = npoly+1;
    poly[ipoly] = 2;
    ipoly = ipoly+1;
    poly[ipoly] = 3;
    ipoly = ipoly+1;
    poly[ipoly] = 4;
    ipoly = ipoly+1;
    npoly = npoly+1;
    poly[ipoly] = 2;
    ipoly = ipoly+1;
    poly[ipoly] = 4;
    ipoly = ipoly+1;
    poly[ipoly] = 5;
    ipoly = ipoly+1;
    npoly = npoly+1;
    poly[ipoly] = 2;
    ipoly = ipoly+1;
    poly[ipoly] = 5;
    ipoly = ipoly+1;
    poly[ipoly] = 6;
    ipoly = ipoly+1;
    npoly = npoly+1;
    for(int i=0; i<n; i++)   // Fourier points
    {
      vert[nvert][0] = x[i];
      vert[nvert][1] = yf[i];
      vert[nvert][2] = 0.0;
      nvert = nvert+1;
    } // end i
    for(int i=0; i<n-1; i++)
    {
      poly[ipoly] = 2; // connect Fourier yf points
      ipoly = ipoly+1;
      poly[ipoly] = i+7; // first 6 vert above
      ipoly = ipoly+1;
      poly[ipoly] = i+8;
      ipoly = ipoly+1;
      npoly = npoly+1;
    } // end i
    System.out.println("nvert="+nvert+"  npoly="+npoly+"  ipoly="+ipoly);
    
    new datwrite("test_fourier.dat", vert, nvert, poly, npoly);
    System.out.println("wrote test_fourier.dat");
  } // end test_fourier
  
  public static void main(String args[])
  {
    new test_fourier();
  }
} // end test_fourier.java