// LV1
// Lotka-Volterra predator-prey model (Euler method)
//   dx/dt =  a*x - b*x*y   (x : prey)
//   dy/dt = -c*y + d*x*y   (y : predator)

float a = 1.0, b = 0.1, c = 1.5, d = 0.075;  // coefficients
float x0 = 10, y0 = 5;                       // initial populations
float dt = 0.001, tmax = 30;                 // time step, end time

// left graph : population vs time
float tx0 = 40, ty0 = 350;    // origin
float kt = 10, kn = 6;        // 1 unit of t = 10 dots, 1 animal = 6 dots
// right graph : phase plane (x - y)
float px0 = 380, py0 = 350;   // origin
float kx = 5, ky = 10;        // 1 prey = 5 dots, 1 predator = 10 dots

void setup() {
  size(1280, 800);
  scale(2);  // draw in 640x400 coordinates, twice the size
  background(0);

  // axes and ticks (every 10 units)
  stroke(0, 255, 0);
  line(tx0, ty0, tx0 + 300, ty0);   // t axis
  line(tx0, ty0, tx0, ty0 - 300);   // N axis
  for (int i = 1; i <= 3; i++) line(tx0 + 10*i*kt, ty0, tx0 + 10*i*kt, ty0 - 5);
  for (int i = 1; i <= 5; i++) line(tx0, ty0 - 10*i*kn, tx0 + 5, ty0 - 10*i*kn);
  line(px0, py0, px0 + 250, py0);   // x axis
  line(px0, py0, px0, py0 - 300);   // y axis
  for (int i = 1; i <= 5; i++) line(px0 + 10*i*kx, py0, px0 + 10*i*kx, py0 - 5);
  for (int i = 1; i <= 3; i++) line(px0, py0 - 10*i*ky, px0 + 5, py0 - 10*i*ky);

  // labels
  fill(0, 255, 255);
  text("Lotka-Volterra model", 240, 24);
  text("1 division = 10", 260, 380);
  text("0", 28, 364);   text("t", 344, 364);   text("N", 36, 44);
  text("0", 368, 364);  text("x", 634, 364);   text("y", 376, 44);
  fill(255, 255, 0);  text("prey x", 120, 70);
  fill(255, 0, 255);  text("predator y", 180, 70);

  // calculation
  float t = 0, x = x0, y = y0;
  while (t <= tmax) {
    float dx = ( a*x - b*x*y) * dt;
    float dy = (-c*y + d*x*y) * dt;
    x = x + dx;
    y = y + dy;
    t = t + dt;
    stroke(255, 255, 0);  point(tx0 + kt*t, ty0 - kn*x);   // prey vs time
    stroke(255, 0, 255);  point(tx0 + kt*t, ty0 - kn*y);   // predator vs time
    stroke(255);          point(px0 + kx*x, py0 - ky*y);   // phase plane
  }
}