/* From: "A SURVEY OF COMPUTATIONAL PHYSICS" 
   by RH Landau, MJ Paez, and CC BORDEIANU 
   Copyright Princeton University Press, Princeton, 2008.
   Electronic Materials copyright: R Landau, Oregon State Univ, 2008;
   MJ Paez, Univ Antioquia, 2008; and CC BORDEIANU, Univ Bucharest, 2008.
   Support by National Science Foundation                              
   */
// java rk4 solution of predator & prey dynamics

import java.io.*;
public class PredatorPrey {
  
  public static void main(String[] argv) throws IOException,
    FileNotFoundException {                        // open  file .dat for output data
    PrintWriter w = new PrintWriter(new FileOutputStream("Pp30.dat"), true);
    PrintWriter q = new PrintWriter(new FileOutputStream("Pp31.dat"), true);
    PrintWriter l = new PrintWriter(new FileOutputStream("Pp32.dat"), true);
    double h, t, Tmin = 0.0, Tmax = 500.0;                               // endpoints
    double y[] = new double[2];  
    int Ntimes=1000;
    y[0]=2.0;  y[1] = 1.3;                                              // initialize
    h = (Tmax-Tmin)/Ntimes; 
    t = Tmin;
    for (t = Tmin;  t <= Tmax; t += h) {
      System.out.println(" t=" +t+" , x= "+y[0]+", v= "+y[1]);//printout
      w.println(""+t +" " +y[0]+" ");//output to files
      q.println(""+t +" " +y[1]+" ");
      l.println(""+y[0] +" " +y[1]+" ");
      rk4(t, y, h, 2);
    }
    System.out.println("Data stored in Pp2.dat");
  }
    
  // ====== rk4 method, *DO NOT MODIFY*, Instead, modify f method 
  public static void rk4(double t, double y[], double h, int Neqs) {
    int i; 
    double F[]  = new double[Neqs]; double ydumb[] = new double[Neqs];
    double k1[] = new double[Neqs]; double k2[] = new double[Neqs]; 
    double k3[] = new double[Neqs]; double k4[] = new double[Neqs];  
    f(t, y, F);
    for (i=0; i<Neqs; i++){ k1[i] = h*F[i]; ydumb[i] = y[i] + k1[i]/2;}
    f(t + h/2, ydumb, F);
    for (i=0; i<Neqs; i++) { k2[i] = h*F[i];  ydumb[i] = y[i] + k2[i]/2;}
    f(t + h/2, ydumb, F);
    for (i=0; i<Neqs; i++) { k3[i]=  h*F[i]; ydumb[i] = y[i] + k3[i];}
    f(t + h, ydumb, F);
    for (i=0; i<Neqs; i++) {
    k4[i] = h*F[i];  y[i] = y[i] + (k1[i] + 2*(k2[i]+k3[i]) + k4[i])/6;}
  }
  
  // FUNCTION of your choice below
  public static void f(double t, double y[], double F[]) { 
    // Lotka-Volterra model  y[0]=1.7 y[1]=1.0 t=[0,500]
    // F[0] = 0.2*y[0]-0.1*y[0]*y[1];                        // RHS of first equation
    // F[1] =-0.1*y[1]+0.1*y[0]*y[1];                       // RHS of second equation
	  // Lotka-Volterra model with prey limitation 
   //  y[0]=1.7 y[1]=1.0 t=[0,500]
   F[0] = 0.2*y[0]*(1-(y[0]/(20.0)))-0.1*y[0]*y[1];          // RHS of first equation
   F[1] =-0.1*y[1]+0.1*y[0]*y[1];                           // RHS of second equation
    //Lotka-Volterra model with functional and numerical responses 
    // y[0]=1.7 y[1]=1.0 t=[0,500] 
    // F[0] = 0.2*y[0]*(1-(y[0]/500.0))-(0.3*y[0]*y[1]/(1+0.3*y[0]*0.5));
    // F[1] = 0.1*y[1]*(1-(y[1]/(0.2*y[0])));
} }