//************************************************************************ 
// Lotka-Volterra predator-prey model with the classical RK4 method
//  x'=ax-bxy
//  y'=cxy-dy, x=x(t), y=y(t), T>t>0, x(0)=x_0, y(0)=y_0 
//************************************************************************
/*
compile with	
gcc -o Lotka Lotka.c  -L/usr/local/dislin -ldislnc -lm -lX11
*/


#include "/usr/local/dislin/dislin.h"
#include <math.h>
#include <stdio.h>


#define PI 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
// constants 
double a=2.;
double b=0.02;
double c = 0.0002;
double d = 0.8;
double ep=0.;

//
double h=0.005; //time step
//
int step=0.0;
int max_steps = 12000.; // final time
//
double x = 5000.0; //initial conditions
double y = 120.0;
//
//
double m_0x = 0.0;
double m_1x = 0.0;
double m_2x = 0.0;
double m_3x = 0.0;
double m_0y = 0.0;
double m_1y = 0.0;
double m_2y = 0.0;
double m_3y = 0.0;
double t=0.0;
//
//define functions for r.h.s.
double funct_x(double t, double x, double y) {
  return a*(1.0+ep*sin(3.*PI*t/2.))*x-b*x*y;
}

double funct_y(double t, double x, double y) {
  return c*x*y-d*y;
}
main ()
{ 
  winsiz (600, 600);
  page (2600, 2600);
  sclmod ("full");
  scrmod ("revers");
  metafl ("xwin");

  disini ();
  pagera ();
  hwfont ();
  name ("x", "x");
  name ("y", "y");

  axspos (350, 2300);
  axslen (2000,2000);
  titlin ("Phase diagram", 4);
  graf (1000.0, 10000.0, 1000., 2000., 40., 180.0, 40., 20.); 

  height (50);
  title ();


setrgb(1,0,0); //plot the fixed points
hsymbl(20);
rlsymb(21,d/c,a/b);

//-----------------------------------------------------------------
// classical RK4 method
for ( step = 0; step < max_steps; step++ ) {

    m_0x = h * funct_x(t, x, y);
    m_0y = h * funct_y(t, x, y);
    
    m_1x = h * funct_x(t+0.5*h, x + 0.5 * m_0x, y + 0.5 * m_0y);
    m_1y = h * funct_y(t+0.5*h, x + 0.5 * m_0x, y + 0.5 * m_0y);
    

    m_2x = h * funct_x(t+0.5*h, x + 0.5 * m_1x, y + 0.5 * m_1y);
    m_2y = h * funct_y(t+0.5*h, x + 0.5 * m_1x, y + 0.5 * m_1y);

    m_3x = h * funct_x(t+h, x + m_2x, y + m_2y);
    m_3y = h * funct_y(t+h, x + m_2x, y + m_2y);

  
    x += (m_0x + 2 * m_1x + 2 * m_2x + m_3x) / 6.0;
    y += (m_0y + 2 * m_1y + 2 * m_2y + m_3y) / 6.0;
    t=t+h;

    setrgb(0,0,1);
    hsymbl(10);
    rlsymb(21,x,y);
   }
//-----------------------------------------------------------------

  disfin ();
}