/* built with Studio Sketchpad: 
 *   https://sketchpad.cc
 * 
 * observe the evolution of this sketch: 
 *   https://studio.sketchpad.cc/sp/pad/view/ro.sVq$IfV5dIS/rev.1792
 * 
 * authors: 
 *   
 *   
 *   Gustav Delius
 *   Owen Petchey
 *   Gian Marco

 * license (unless otherwise specified): 
 *   creative commons attribution-share alike 3.0 license.
 *   https://creativecommons.org/licenses/by-sa/3.0/ 
 */ 



int N = 250; // largest possible prey number
int M = 150; // largest possible predator number (we don't need it)
int n0 = 125; // initial prey number
int m0 = 125; // initial predator number

// functional response parameters
float a = 0.08;
// why does a handling time (h) (or b) of 0 result in prey extinction?
// I (Owen) would think that h=0 gives a linear functional response
// which at least for a deterministic predator-prey model will
// lead to stable limit cycles.
// Note that we change the death rate function of the predators.
float b;
float h = 1000;
float e = 0.1;

float lambdan(int n, int m) {
  float lambda1 = 0.06; // intrinsic prey birth rate constant
  return n*lambda1*(1-n/N);
}
  
float mun(int n, int m) {
  float mu1 = 0.01; // intrinsic prey deat rate constant
  return mu1*n*(1+n/N) + a*n*m/(1+b*n);    
}
float lambdam(int n, int m) {
  float lambda2 = 0.05; // intrinsic predator birth rate constant
  return  /* m*lambda2*(1-m/M)+ */  e*a*n*m/(1+b*n);
}
  
float mum(int n, int m) {
  float mu2 = 0.035; // intrinsic predator death rate constant
  // the following is a predator density dependent predator death rate
  //return mu2*m*(1+m/M);    
  // we may also use a density independent predator death rate
  return mu2*m;    
  
  
}

int numWalkers = 100;
Walker[] w1;
//Walker[] w2;

void setup () {
  b = a*h;
  
  randomSeed(100);
  size (600, 300);
  smooth ();
  frameRate(100);
  fill(0);
  textAlign(CENTER, CENTER);
  
  // Divide canvas into two panels
  line(280, 0, 280, height);
  
  // Phase plane panel
  // Draw vertical axis
  line(25, 20, 25, height);
  // label vertical axis
  text("Predator", 25, 7);
  for (int i=50; i<300; i = i+50) {
    line(20, height-20-i, 30, height-20-i);
    text(i, 8, height-20-i);
  }
  // Draw horizontal axis
  line(0, height-20, 270, height-20);
  // label horizontal axis
  text("Prey", 255, height-8);
  for (int i=50; i<250; i = i+50) {
    line(25+i, height-25, 25+i, height-15);
    text(i, 25+i, height-8);
  }
  
  // Time evolution panel
  // Draw vertical axis
  line(25+290, 10, 25+290, height);
  // label vertical axis
  stroke(20, 200, 60);
  line(350, 7, 380, 7);
  text("Prey", 400, 7);
  stroke(255, 0, 0);
  line(460, 7, 490, 7);
  stroke(0, 0, 0);
  text("Predator", 520, 7);
  for (int i=50; i<300; i = i+50) {
    line(20+290, height-20-i, 30+290, height-20-i);
    text(i, 8+290, height-20-i);
  }
  // Draw horizontal axis
  line(0+290, height-20, width-10, height-20);
  // label horizontal axis
  text("Time", width-25, height-8);
  for (int i=50; i<250; i = i+50) {
    line(25+i+290, height-25, 25+i+290, height-15);
    text(i, 25+i+290, height-8);
  }
  
  // create many walkers and put them into the array w
  w1 = new Walker[numWalkers]; 
  //w2 = new Walker[numWalkers]; 
  for (int i=0; i<numWalkers; i++) {
    w1[i] = new Walker(0, n0, m0, 1); 
    //w2[i] = new Walker(0, 150, 2); 
  }
  
}

void draw () {
  // update all walkers
  for (int i=0; i<numWalkers; i++) {
    w1[i].update();
    //w2[i].update();
  }
}

class Walker {
  float t;  // time
  int n, m; // predator and prey numbers
  int type; // type of walker
  
  Walker (float tCon, int nCon, int mCon, int typeCon) {
    type = typeCon;
    t = tCon;
    n = nCon;
    m = mCon;
  }

  void update () {
    // do nothing if predator or prey extinct
    if ((n==0) || (m==0)) {
        return;
    }

    float bn = lambdan(n,m);
    float dn = mun(n,m);
    float bm = lambdam(n,m);
    float dm = mum(n,m);
    float h = bn + dn + bm + dm;
    
    // determine time of next process
    float tn = t - log(random(1))/h;

    // initialise population values for next timestep to default
    int nn = n;
    int mn = m;

    // randomly determine which process takes place;
    float r = random(1);
    if (r<dn/h) {
      nn = n-1;
    } else if (r<(dn+bn)/h) {
      nn = n+1;
    } else if (r < (dn + bn + dm)/h) {
      mn = m-1;
    } else {
      mn = m+1;
    }

    // update time series plot
    // draw prey in green
    stroke(20, 200, 60, 10);
    line(t+315, height-20-n, tn+315, height-20-nn);
    // draw predator in red
    stroke(255, 0, 0, 10);
    line(t+315, height-20-m, tn+315, height-20-mn);
    
    // update phase space plot
    stroke(0, 0, 255, 10);
    line(n+25, height-20-m, nn+25, height-20-mn);
    
    // update location
    t = tn;
    n = nn;
    m = mn;
  }
}