#include "ibex.h"
#include "vibes.cpp"


using namespace ibex;


void plot_simu(const simulation* sim)
{
  std::list<solution_g>::const_iterator iterator_list;
  for(iterator_list=sim->list_solution_g.begin();iterator_list!=sim->list_solution_g.end();iterator_list++) {
    vibes::drawBox(iterator_list->time_j.lb(), iterator_list->time_j.ub(),
		   iterator_list->box_j1->operator[](0).lb(), iterator_list->box_j1->operator[](0).ub(), "red[red]");

    vibes::drawBox(iterator_list->time_j.lb(), iterator_list->time_j.ub(),
		   iterator_list->box_j1->operator[](1).lb(), iterator_list->box_j1->operator[](1).ub(), "blue[blue]");

    vibes::drawBox(iterator_list->time_j.lb(), iterator_list->time_j.ub(),
		   iterator_list->box_j1->operator[](2).lb(), iterator_list->box_j1->operator[](2).ub(), "black[black]");
  }
}


int main(){

  const int n= 3;
  Variable y(n);

  IntervalVector yinit(n);
  yinit[0] = Interval(9.98);
  yinit[1] = Interval(0.01);
  yinit[2] = Interval(0.01);
  Interval a(0.3);

  Function ydot = Function(y,Return(-y[0]*y[1]/(1+y[0]),
				    y[0]*y[1]/(1+y[0]) - a*y[1],
				    a*y[1]));

  NumConstraint csp1(y,y[0]+y[1]+y[2] -10.0 = 0);
  NumConstraint csp2(y,y[0]>=0);
  NumConstraint csp3(y,y[1]>=0);
  NumConstraint csp4(y,y[2]>=0);

  Array<NumConstraint> csp(csp1,csp2,csp3,csp4);

  ivp_ode problem = ivp_ode(ydot,0.0,yinit,csp,AUTODIF);

  simulation simu = simulation(&problem,100.0,KUTTA3,1e-6);

  simu.active_monotony();

  simu.run_simulation();

  vibes::beginDrawing ();
  vibes::newFigure("PDS");
  vibes::setFigureProperty("plot","time",10);
  vibes::setFigureProperty("plot","y",10);
  vibes::setFigureProperty("plot","width",600);
  vibes::setFigureProperty("plot","height",600);

  plot_simu(&simu);

  vibes::endDrawing();

  return 0;

}