#include "dislin.h"
#include <iostream>
#include <stdio.h>

//Quelques constantes pour l'equation de Lotka-Volterra
double alpha = 7;
double beta  = 1;
double gamma = 80;
double delta = 1.8;

//Un facteur por modifier les constantes
double f = 0.004;

//Premiere equation de Lotka-Volterra
double dxdt(double x, double y) {
    return f*(+ alpha*x - beta*x*y);
}

//Deuxieme equation de Lotka-Volterra
double dydt(double x, double y) {
	return f*(- gamma*y + delta*x*y); 
}

int main() {

	//Tableaux pour les temps
	double t[500] = {0};

	//Tableaux pour les solutions utilisant le methode de Runge
	double x[500] = {0};	//Nombre de proies
	double y[500] = {0};	//Nombre de predateurs

	//Tableaux pour les solutions utilisant le methode d'Euler
	double x_Euler[500] = {0}; 
	double y_Euler[500] = {0};
	
	//Conditions initiales
	x[0] = 60;
	y[0] = 6;
	
	//Meme conditions initiales pour la solution utilisant le methode d'Euler
	x_Euler[0] = 60;
	y_Euler[0] = 6;
	
	//Boucle dans le temps
    for (int time=1; time<500; ++time) {

		//Definissez le temps
		t[time] = time;
		
		//Euler	
		x_Euler[time] = x_Euler[time-1] + dxdt(x_Euler[time-1], y_Euler[time-1]);
		y_Euler[time] = y_Euler[time-1] + dydt(x_Euler[time-1], y_Euler[time-1]);

		//Runge
		double xm = x[time-1] + 0.5*dxdt(x[time-1],y[time-1]);
		double ym = y[time-1] + 0.5*dydt(x[time-1],y[time-1]);
		x[time] = x[time-1] + dxdt(xm,ym);
		y[time] = y[time-1] + dydt(xm,ym);
	}

	//Ecrivez les constants aux graphiques
	char constants[40];
	sprintf(constants, "alpha=%.3f, beta=%.3f, gamma=%.3f, delta=%.3f", alpha*f, beta*f, gamma*f, delta*f);

	//Premiere graphique
	metafl("PDF");	disini();
	name("temps t (unites arbitraires)", "x");	
	name("Nombres de proies x et de predateurs y", "y");	
	graf(0.,500.,0.,100.,0.,70.,0.,10.);
	titlin("Solution to the Lotka-Volterra equation",-1);
	titlin(constants,2);	
	titlin("BLEU : nombre de proies",3);
	titlin("ROUGE : nombre de predateurs",4);
	title();	
	color("blue");	curve(t, x, 500);	
	color("red");	curve(t, y, 500);
	disfin();
	
	//Deuxieme graphique
	metafl("PDF");	disini();
	name("nombre de proies x", "x");	
	name("nombre de predateurs y", "y");	
	graf(20.,70.,20.,5.,0.,25.,0.,5.);
	titlin("Solution to the Lotka-Volterra equation",-1);
	titlin(constants,2);	
	titlin("TRAIT PLEIN : Methode de Runge", 3);
	titlin("TRAIT POINTILLE : Methode d'Euler (insuffisante!)", 4);
	title();
	curve(x, y, 500);
	dash(); curve(x_Euler, y_Euler, 500);
	disfin();

	return 0;
}