// A program to model the Lotka Volterra model using the fourth-order Runge-Kutta algorithm.

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

public class Lotka_RK {
    public static void main(String argv[]) {
  	double xtemp,ytemp;
  	double dt = 0.1;
  	int niter=6000;
  	StdDraw.clear(StdDraw.LIGHT_GRAY);
    StdDraw.setXscale(0, 50);
    StdDraw.setYscale(0, 50);
    
    FileOutputStream out;
    PrintStream p;
    
    double y1[] = new double[niter+1];
    double y2[] = new double[niter+1];
    double y[] = new double[2];
    
 // Set up time step and initial values
    
    y1[0] = y[0] = 50;				//Rabbits
    y2[0] = y[1] = 10;				//Foxes
    
      
 // 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];
      
      xtemp = y[0];
      ytemp = y[1];
      
 // Writing Data to the File "OutputLotka.txt"
    try
    {
     out = new FileOutputStream("OutputLotka.txt",true);
     p = new PrintStream(out);            
     p.println(y[0]+"    "+y[1]+" "+dt*i);           
     p.close();
    }
    catch(Exception e)
    {
      System.err.println("Error in printing to file");
    }     
            
      StdDraw.setPenColor(StdDraw.RED);
      StdDraw.point(t, ytemp);
      StdDraw.setPenColor(StdDraw.BLUE);
      StdDraw.point(t, xtemp);
            
      StdDraw.setPenColor(StdDraw.YELLOW);
      StdDraw.point(xtemp,ytemp);
      
      if (i % 5 == 0) {
      	StdDraw.show(10);
      	System.out.println("t: " + dt*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);
    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 of population

  public static double[] g(double y[]) {
    int l = y.length;
    double a = 0.04, b = 0.005, c = 0.02, e = 0.5;
    double v[] = new double[l];
    v[0] = y[0]*(a-b*y[1]);
    v[1] = (c*y[0]-e)*y[1];
    return v;
  }
}