///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Pendulum.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.     //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

// A program to study the driven pendulum under damping
// via the fourth-order Runge-Kutta algorithm.

import java.lang.*;
import java.io.*;

public class Lion_RK {
    public static void main(String argv[]) {
  	double xtemp,ytemp, ztemp;
  	double dt = 0.1;
  	int niter=2000;
  	StdDraw.clear(StdDraw.LIGHT_GRAY);
    StdDraw.setXscale(0, 200);
    StdDraw.setYscale(0, 800);
    
    FileOutputStream out;
    PrintStream p;
    
    double y1[] = new double[niter+1];
    double y2[] = new double[niter+1];
    double y3[] = new double[niter+1];
    double y[] = new double[3];
    
 // Set up time step and initial values
    
    y1[0] = y[0] = 4; 		//lions
    y2[0] = y[1] = 110;     //foxes
    y3[0] = y[2] = 2200;    //rabbits
    
    
 // Perform the 4th-order Runge-Kutta integration
    for (int i=0; i<niter; ++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];
      
      xtemp = y[0];
      ytemp = y[1];
      ztemp = y[2];
      
    // Writing Data to the File "OutputLion.txt"
   if(i % 10 == 0){
    
   try
    {
     out = new FileOutputStream("OutputLion.txt",true);
     p = new PrintStream(out);
              
     p.println(y[0]+"    "+y[1]+" "+y[2]+" "+dt*i);           
     p.close();
    }
    catch(Exception e)
    {
      System.err.println("Error in printing to file");
    }     
    }    
         
      StdDraw.setPenColor(StdDraw.RED);
      StdDraw.point(t, xtemp);
   
      StdDraw.setPenColor(StdDraw.RED);
      StdDraw.point(t, ytemp);
    
      StdDraw.setPenColor(StdDraw.GREEN);
      StdDraw.point(t, ztemp);  
           
  /*  StdDraw.setPenColor(StdDraw.BLUE);
      StdDraw.point(xtemp,ztemp);
    
      StdDraw.setPenColor(StdDraw.WHITE);
      StdDraw.point(xtemp,ytemp);
      
      StdDraw.setPenColor(StdDraw.MAGENTA);
      StdDraw.point(ytemp,ztemp);
    */
      
      if (i % 10 == 0) {
      	StdDraw.show(10);
      	System.out.println("i: " + i);
      	System.out.println("xtemp: " + xtemp);
      	System.out.println("ytemp: " + ytemp);
      	System.out.println("ztemp: " + ztemp);
      	}
    }
    
  }

// 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);
    for (int i=0; i<l; ++i) c2[i] = y[i] + dt*c1[i]/2;
  	c2 = g(c2);
    for (int i=0; i<l; ++i) c3[i] = y[i] + dt*c2[i]/2;
  	c3 = g(c3);
    for (int i=0; i<l; ++i) c4[i] = y[i] + dt*c3[i];
  	c4 = g(c4);
    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 derivatives

  public static double[] g(double y[]) {
    int l = y.length;
  
    double a1 = 0.009, a2 = 0.0009, a3 = 2.7;
    double a4 = 0.0016, a5 = 2.55, a6 = 0.01;
    double a7 = 5.15, a8 = 0.03, a9 = 0.3;
 
    double v[] = new double[l];
    v[0] = y[0]*(a1*y[1]+a2*y[2]-a3);
    v[1] = y[1]*(a4*y[2]-a5-a6*y[0]);
    v[2] = y[2]*(a7-a8*y[1]-a9*y[0]);
    return v;
  }
}