//************************************************************************ 
// Lotka-Volterra competition model with the classical RK4 method
//  x'=ax(b-x-cy)
//  y'=dy(e-y-fx), x=x(t), y=y(t), T>t>0, x(0)=x_0, y(0)=y_0 
//************************************************************************
/*
compile with	
gcc -o competitionRK4 competitionRK4.c  -L/usr/local/dislin -ldislnc -lm -lX11
*/


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

#define NX 25 
#define NY 25
#define L 100
//
int step = 0.0;
int max_steps = 100000.0; // final time
//
double x = 70.0; //initial conditions
double y = 40.0;
//
double h = 0.001; // time step
//
double a = 0.004; // constants of the problem
double b = 50.0;
double c = 0.75;
double d = 0.001;
double e = 100.0;
double f = 3.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;
//
float xv[NX][NY], yv[NX][NY], xp[NX], yp[NY];

//define functions for r.h.s.
double funct_x(double x, double y) {
  return a*x*(b-x-c*y);
}

double funct_y(double x, double y) {
  return d*y*(e-y-f*x);
}

main ()
{ int i, j, iret , iclr;

  double xstep, ystep;
// define grid for a vector plot
  xstep = L / (NX - 1);; 
  ystep = L / (NY - 1);; 
  for (i = 0;  i < NX; i++)
  { xp[i] = (float) (-1. +  i * xstep);  
    for (j = 0; j < NY; j++) 
    { yp[j] = (float) (-1. +  j *ystep); 
      xv[i][j] = (float) a*xp[i]*(b-xp[i]-c*yp[j]);
      yv[i][j] = (float) d*yp[j]*(e-yp[j]-f*xp[i]);
    }
  }
  

  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 (0., L, 0., 10., 0., L, 0., 10.); 

  height (50);
  title ();

 setrgb(0,0,1); // set blue color

 vecmat ((float *) xv, (float *) yv, NX, NY, xp, yp, -1); // plot vector field

setrgb(1,0,0); //plot the fixed points
hsymbl(50);
rlsymb(15,0,0);
rlsymb(21,50,0);
rlsymb(21,0,100);
rlsymb(15,20,40);
//-----------------------------------------------------------------
// classical RK4 method
for ( step = 0; step < max_steps; step++ ) {

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

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

    m_3x = h * funct_x(x + m_2x, y + m_2y);
    m_3y = h * funct_y(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;

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

 /* stmmod ("rk4", "integration");
  stmmod ("on", "close");
  color ("blue");
  iclr = intrgb (0.0, 0.0, 0.0);
  vecclr (iclr);
  stream ((float *) xv, (float *) yv, NX, NY, xp, yp, NULL, NULL, 0); //plot the stream
*/
  disfin ();
}