//etude du refroidissement d'une boisson
#include <iostream>
#include <cmath>
#include "dislin.h"     //bibliotheque DISLIN

using namespace std;

int main() {
  const int STEPMAX = 100;  //nombre de pas
  double dt;                //pas d'integration
     
  double tempExt = 10.;		//temperature exterieur
  double temp[STEPMAX];		//tableau des temperatures
  double t[STEPMAX];
  double k = 0.7;       //coefficient de refroidissement degres / h
  
  //conditions intiales
  temp[0] = 20.;
  t[0] = 0.;

  //solution methode de Runge
  cout << "Temperature exterieur constante" << endl;
  dt = 0.1;
  double tempM;
  bool cond1 = false;
  for (int step=0; step<STEPMAX-1; step++) {
    //pour calculer temp dans (t_i + Delta t/2) on utilse la methode d'Euler
    tempM = temp[step] - k*(temp[step]-tempExt)*dt/2.;    
    temp[step+1] = temp[step] - k*(tempM-tempExt)*dt;
    t[step+1] = t[step] + dt;
    if (!cond1 && temp[step+1] < 12.) {
      cout << "Apres t = " << t[step+1] << 
         " heurs la temperature de la boisson a baisse a 12 degres !" << endl;
      cond1 = true;
    }
  }

  //initialisation de DISLIN
  metafl("XWIN");   //XWIN ou PDF
  disini();
  smxalf("GREEK", "{", "}", 1);
  name("temps (h)", "X");
  name("T(t)", "Y");
  titlin("Refroidissement", 1);
  titlin("dT/dt = -k*(T(t)-T_ext)",2);
  titlin("T(0) = 20, k = 0.7 / h", 3);
  graf(0., 10., 0., 1., 5., 25., 5., 1.);
  title();
 
  //dessin
  color("RED");
  incmrk(-1);   //symboles au lieu de la curbe
  marker(21);   //choix du symbole
  curve(t,temp,100);

  //temperature externe variable
  cout << "\nTemperature exterieur variable" << endl;
  dt = 0.1;
  double tempB[STEPMAX];    //temperatures
  //variation de la temperature exterieur : 
  //elle chute de 4 degres sur 8 heures
  double a = 4./8.;
  
  //conditions intiales
  tempB[0] = 20.;
  t[0] = 0.;
  
  bool cond2 = false;
  bool cond3 = false;
  for (int step=0; step<STEPMAX-1; step++) {
    //pour calculer temp dans (t_i + Delta t/2) on utilse la methode d'Euler
    tempM = tempB[step] - k*(tempB[step]-(tempExt-a*t[step]))*dt/2.;    
    tempB[step+1] = tempB[step] - k*(tempM-(tempExt-a*t[step]))*dt;
    t[step+1] = t[step] + dt;
    if (!cond2 && t[step+1]>0.333) {
      cout << t[step+1]*60 << " min apres 22h00, la temperature est de "
	       << tempB[step+1] << " degres" << endl;
      cond2 = true;
    }
    if (!cond3 && t[step+1]>8.) {
	  cout << "Apres 06h00 la temperature de la boisson a chute a " 
	       << tempB[step+1] << " degres" << endl;
	  cond3 = true;
    }
  }

  //dessin
  color("GREEN");
  incmrk(-1);   //symboles au lieu de la curbe
  marker(21);   //choix du symbole
  curve(t,tempB,100);
  
  //temperature exterieur variable - methode d'EULER
  cout << "\nTemperature exterieur variable - methode d'Euler" << endl;
  dt = 0.1;
  double tempE[STEPMAX];    //temperatures

  //conditions intiales
  tempE[0] = 20.;
  t[0] = 0.;

  cond2 = false;
  cond3 = false;
  for (int step=0; step<STEPMAX-1; step++) {
    tempE[step+1] = tempE[step] - k*(tempE[step]-(tempExt-a*t[step]))*dt;
    t[step+1] = t[step] + dt;
    if (!cond2 && t[step+1]>0.333) {
      cout << t[step+1]*60 << " min apres 22h00, la temperature est de "
           << tempE[step+1] << " degres" << endl;
      cond2 = true;
    }
    if (!cond3 && t[step+1]>8.) {
      cout << "Apres 06h00 la temperature de la boisson a chute a " 
           << tempE[step+1] << " degres" << endl;
	  cond3 = true;
    }
  }

  //dessin
  color("BLUE");
  incmrk(-1);   //symboles au lieu de la curbe
  marker(21);   //choix du symbole
  curve(t,tempE,100);

  //fin de dislin
  disfin();

  return 0;
}