/*************************************************************************/
/* Program file name: Ex4_8.java                                         */
/*                                                                       */
/* © Tao Pang 2006                                                       */
/*                                                                       */
/* Last modified: June 15, 2006                                          */
/*                                                                       */
/* (1) This Java program is created for 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.     */
/*                                                                       */
/*************************************************************************/

// A program to study the Lorenz model via the fourth-order
// Runge-Kutta algorithm.

import java.lang.*;
public class Ex4_8 {
  static final int n = 100, nt = 10, m = 5;
  public static void main(String argv[]) {
    double y1[] = new double[n+1];
    double y2[] = new double[n+1];
    double y3[] = new double[n+1];
    double y[] = new double[3];

 // Set up time step and initial values
    double dt = 2*Math.PI/nt;
    y1[0] = y[0] = 0;
    y2[0] = y[1] = 1;
    y3[0] = y[2] = 2;

 // Perform the 4th-order Runge-Kutta integration
    for (int i=0; i<n; ++i) {
      double t = dt*i;
      y = rungeKutta(y, t, dt); 
      y1[i+1] = y[0];
      y2[i+1] = y[1];
      y3[i+1] = y[2];
    }

 // Output the result in every m time steps
    for (int i=0; i<=n; i+=m)
      System.out.println(y1[i] + " " + y2[i] + " " + y3[i]);
  }

// Method to complete one Runge-Kutta step.

  public static double[] rungeKutta(double y[],
    double t, double dt) {
    int l = y.length;
    double c1[] = new double[l];
    double c2[] = new double[l];
    double c3[] = new double[l];
    double c4[] = new double[l];
    c1 = g(y, t);
    for (int i=0; i<l; ++i) c2[i] = y[i] + dt*c1[i]/2;
    c2 = g(c2, t+dt/2);
    for (int i=0; i<l; ++i) c3[i] = y[i] + dt*c2[i]/2;
    c3 = g(c3, t+dt/2);
    for (int i=0; i<l; ++i) c4[i] = y[i] + dt*c3[i];
    c4 = g(c4, t+dt);
    for (int i=0; i<l; ++i)
      c1[i] = y[i] + dt*(c1[i]+2*(c2[i]+c3[i])+c4[i])/6;
    return c1;
  }

// Method to provide the generalized velocity vector.

  public static double[] g(double y[], double t) {
    int l = y.length;
    double a0 = 0.5, b0 = 1.0, c0 = 1.5;
    double v[] = new double[l];
    v[0] = a0*(y[1]-y[0]);
    v[1] = (b0-y[2])*y[0] - y[1];
    v[2] = y[0]*y[1] - c0*y[2];
    return v;
  }
}