entity Rabbits {
  property real alpha; // = 0.1;
  property real beta;  // = 0.01;
  property real pop;
  property real dpop;
 
  service real getPop() {
  	return pop;
  }
  service setPop(real newPop) {
  	pop = newPop;
  }
  service real diff_birth(real x) { 
  	return alpha*x;
  }
  service real diff_eaten(real x, real y) {
  	return -beta*x*y;
  }
}

entity Foxes {
  property real gamma; // = 0.05; 
  property real delta; // = 0.001;
  property real pop;
  property real dpop;

  service real getPop() {
  	return pop ;
  }
  service setPop(real newPop) {
  	pop = newPop;
  }
  service real diff_death(real y) {
  	return -gamma*y;
  }
  service real diff_eating(real y, real x) {
  	return delta*y*x;
  }
}

relation Predation[a, b] {

  service euler_evpop(real delta_t) {
    real d_a=0; 
    real d_b=0;
    real xPop = a.getPop();
    real yPop = b.getPop();
    d_a = d_a + a.diff_birth(xPop)*delta_t;
    d_b = d_b + b.diff_death(yPop)*delta_t;
    d_a = d_a + a.diff_eaten(xPop,yPop)*delta_t;
    d_b = d_b + b.diff_eating(yPop,xPop)*delta_t;
    a.setPop(xPop + d_a);
    b.setPop(yPop + d_b);
  }
  
  service rk4_evpop(real delta_t) {
    real d_a=0;
    real d_b=0;
    real k1x ; real k2x ; real k3x; real k4x;
    real k1y ; real k2y ; real k3y; real k4y;
    real xPop = a.getPop();
    real yPop = b.getPop();

    k1x = a.diff_birth(xPop) + a.diff_eaten(xPop,yPop);
    k1y = b.diff_death(yPop) + b.diff_eating(yPop,xPop);
    k2x = a.diff_birth(xPop+k1x*delta_t/2) + a.diff_eaten(xPop+k1x*delta_t/2, yPop+k1y*delta_t/2);
    k2y = b.diff_death(yPop+k1y*delta_t/2) + b.diff_eating(yPop+k1y*delta_t/2, xPop+k1x*delta_t/2);
    k3x = a.diff_birth(xPop+k2x*delta_t/2) + a.diff_eaten(xPop+k2x*delta_t/2, yPop+k2y*delta_t/2);
    k3y = b.diff_death(yPop+k2y*delta_t/2) + b.diff_eating(yPop+k2y*delta_t/2, xPop+k2x*delta_t/2);
    k4x = a.diff_birth(xPop+k3x*delta_t) + a.diff_eaten(xPop+k3x*delta_t, yPop+k3y*delta_t);
    k4y = b.diff_death(yPop+k3y*delta_t) + b.diff_eating(yPop+k3y*delta_t, xPop+k3x*delta_t);
    d_a = d_a + (k1x+2*k2x+2*k3x+k4x)*delta_t/6;
    d_b = d_b + (k1y+2*k2y+2*k3y+k4y)*delta_t/6;	

    a.setPop(xPop + d_a);
    b.setPop(yPop + d_b);
  }
  
}

scenario Main{

    real dt = 0.1; 

    rabbits1 = Rabbits{pop = 50; alpha = 0.1; beta  = 0.01;};
    foxes1 = Foxes{pop = 15; gamma = 0.05; delta = 0.001;};
    rp = Predation[Rabbits, Foxes];
    rp.connect(rabbits1,foxes1);    	
    for (i in 1..2000){
       //rp.euler_evpop(dt);
       rp.rk4_evpop(dt);
       print(i);
       print(rabbits1.getPop()); 
       print(foxes1.getPop());
       print("-");
    }
}