// ode_spiral_chk2.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)
//
// 
import java.text.*;
import java.io.*;

class ode_spiral_chk2
{
  double tpi = 2.0*Math.PI;
  int i = 0;
  int n = 10; // number of points to compute  
    
  ode_spiral_chk2()
  {
    double t;
    double x[] = new double[256];
    double y[] = new double[256];
    double xpv[] = new double[256];
    double ypv[] = 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_chk2_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;
    while(true)
    {
      if(i>=n) { break; }	
      k1 = cx(t);
      k2 = cx(t+dt/2.0);
      k3 = cx(t+dt);
      xpv[i+1] = xpv[i]+dt*(k1+2.0*k2+k3)/4.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/2.0);
      k3 = cy(t+dt);
      ypv[i+1] = ypv[i]+dt*(k1+2.0*k3+k2)/4.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


    // datwrite
    System.out.println("finished ode_spiral_chk2.java");
  } // end ode_spiral_chk2


  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_chk2();
  }
} // end ode_spiral_chk2.java