// ode_spiral_chk.java  initial value ODE using Runge Kutta
//                  known solution to test method
//
// x(t)   = xf(t)  = t*cos(tpi*t)
// y(t)   = yf(t)  = t*sin(tpi*t)
// x'(t)  = xp(t)  = cos(tpi*t) - t*tpi*sin(tpi*t)
// y'(t)  = yp(t)  = sin(tpi*t) + t*tpi*cos(tpi*t)
// x''(t) = xpp(t) = -2*tpi*sin(tpi*t) - t*tpi*tpi*cos(tpi*t)
// y''(t) = ypp(t) = 2*tpi*cos(tpi*t) - t*tpi*tpi*sin(tpi*t)
// cx(t) = x''(t)  = -2*tpi*sin(tpi*t) - t*tpi*tpi*cos(tpi*t)
// cy(t) = y''(t)  =  2*tpi*cos(tpi*t) - t*tpi*tpi*sin(tpi*t)
//
//     differential equations to solve, given initial conditions
//     d^2x/dt^2 = cx(t)
//     d^2y/dt^2 = cy(t)
//
// Initialize t = 0.1 Runge Kutta 4th order method
//  loop:     use yp = y for exp    fourth order method for y' = y  y=exp(x)
//    k1 = h*y                yp = f(x,y) in general
//    k2 = h*y+k1/2                f(x+h/2.0,y+k1/2.0)
//    k3 = h*y+k2/2                f(x+h/2.0,y+k2/2.0)
//    k4 = h*y+k3                  f(x+h,y+k3)
//    y = y+(k1+2.0*k2+2.0*k3+k4)/6.0
//    x = x+h                      print y at this x
// 
import java.text.*;
import java.io.*;

class ode_spiral_chk
{
  double tpi = 2.0*Math.PI;
  int i = 0;
  int n = 10; // number of points to compute  
    
  ode_spiral_chk()
  {
    double t;
    double x[] = new double[256];
    double y[] = new double[256];
    double xpv[] = new double[256];
    double ypv[] = new double[256];
    double xppv[] = new double[256];
    double yppv[] = new double[256];
    double t0 = 1.0;
    double x0 = xf(t0);
    double y0 = yf(t0);
    double xp0 = xp(t0);
    double yp0 = yp(t0);
    double dt = 0.01;
    double k1, k2, k3, k4;
    double err = 0; // used for debugging
    
    System.out.println("ode_spiral_chk_eq.java running");
    System.out.println("differential equations to solve");
    System.out.println("d^2x/dt^2 = cx(t)");
    System.out.println("d^2y/dt^2 = cy(t)");
    System.out.println("n = "+n);
    System.out.println("t0 = "+t0);
    System.out.println("dt = "+dt);
    System.out.println("x0 = "+x0);
    System.out.println("y0 = "+y0);
    System.out.println("xp0 = "+xp0);
    System.out.println("yp0 = "+yp0);
    System.out.println(" ");
    i = 0;                   // known initial conditions
    t = t0;
    x[i] = x0;
    y[i] = y0;
    xpv[i] = xp0;
    ypv[i] = yp0;
    xppv[i] = cx(t);
    yppv[i] = cy(t);
    while(true)
    {
      xppv[i] = xpp(t); // just for testing
      yppv[i] = ypp(t);
      if(i>=n) { break; }	
      k1 = cx(t);
      k2 = cx(t+dt);
      xpv[i+1] = xpv[i]+dt*(k1+k2)/2.0;
      err = xpv[i+1]-xp(t+dt);
      System.out.println("xpv="+xpv[i+1]+"  xp="+xp(t+dt));
      System.out.println("xpv["+(i+1)+"] err="+err);
      System.out.println(" ");      

      k1 = xpv[i];
      k2 = xpv[i+1];
      x[i+1] = x[i]+dt*(k1+k2)/2.0;
      err = x[i+1]-xf(t+dt);
      System.out.println("x="+x[i+1]+"  xf="+xf(t+dt));
      System.out.println("x["+(i+1)+"] err="+err);
      System.out.println(" ");      

      k1 = cy(t);
      k2 = cy(t+dt);
      ypv[i+1] = ypv[i]+dt*(k1+k2)/2.0;
      err = ypv[i+1]-yp(t+dt);
      System.out.println("ypv="+ypv[i+1]+"  yp="+yp(t+dt));
      System.out.println("ypv["+(i+1)+"] err="+err);
      System.out.println(" ");      

      k1 = ypv[i];
      k2 = ypv[i+1];
      y[i+1] = y[i]+dt*(k1+k2)/2.0;
      err = y[i+1]-yf(t+dt);
      System.out.println("y="+y[i+1]+"  yf="+yf(t+dt));
      System.out.println("ypv["+(i+1)+"] err="+err);
      System.out.println(" ");      
      i = i + 1;
      t = t + dt;

    } // end while

    t = t0;
    for(int i=0; i<n; i++)
    {
      System.out.println("i="+i+"  t="+t);
      System.out.println("x[i]= "+x[i]+" =? "+(xf(t)));
      System.out.println("xpv[i]= "+xpv[i]+" =? "+(xp(t)));
      System.out.println("xppv[i]= "+xppv[i]+" =? "+(xpp(t)));
      System.out.println("y[i]= "+y[i]+" =? "+(yf(t)));
      System.out.println("ypv[i]= "+ypv[i]+" =? "+(yp(t)));
      System.out.println("yppv[i]= "+yppv[i]+" =? "+(ypp(t)));      
      t = t+dt;
      System.out.println(" ");
    } // end i
    // datwrite
    System.out.println("finished ode_spiral_chk.java");
  } // end ode_spiral_chk


  double xf(double t)
  {
      return t*Math.cos(tpi*t);
  } // end xf

  double yf(double t)
  {
      return t*Math.sin(tpi*t);
  } // end ypp

  double xp(double t)
  {
    return Math.cos(tpi*t) - t*tpi*Math.sin(tpi*t);
  } // end xp

  double yp(double t)
  {
      return Math.sin(tpi*t) + t*tpi*Math.cos(tpi*t);
  } // end yp

  double xpp(double t)
  {
    return -2.0*tpi*Math.sin(tpi*t) - t*tpi*tpi*Math.cos(tpi*t);
  } // end xpp

  double ypp(double t)
  {
      return 2.0*tpi*Math.cos(tpi*t) - t*tpi*tpi*Math.sin(tpi*t);
  } // end ypp


  double cx(double t)
  {
    return -2.0*tpi*Math.sin(tpi*t) - t*tpi*tpi*Math.cos(tpi*t);
  } // end cx

  double cy(double t)
  {
      return 2.0*tpi*Math.cos(tpi*t) - t*tpi*tpi*Math.sin(tpi*t);
  } // end cy

  public static void main (String[] args)
  {
    new ode_spiral_chk();
  }
} // end ode_spiral_chk.java