// Attention, pour compiler utiliser :
// g++ -std=c++98 -O2 odoreverse.cpp -o odoreverse -lpng
// S'il y a une erreur concernant lpng, il faut installer la librairie libpng-devel, sous Ubuntu : apt-get install libpng-dev

#include <ctime>
#include <cmath>
#include <algorithm>
#include <iostream>
#include "TDateTime.h"
#include "TStringList.h"
#include "GestionFichiers.h"
#include "Meteo.h"
#include "WritePng.h"

#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
#define __int64 long long int

//---------------------------------------------------------------------------
double RRRsigvc,RRRsigvn,RRRsigy,RRRsigv,RRRustar,RRRhcl,RRRwstar;
double resol_sigyCBL_ADMS(double t);
double resol_sigySTAB_ADMS(double t);
double zbrent(double (*func)(double), double x1, double x2, double tol, int *erreur);
// Modification 24 mars 2020
int calLmoUstar(double z0,double alb,double ppt,double vit,double zvit, double temp,double ray,double neb,double *lmo,double *ustar);
double coeffA(int pgt);

void Log(String nom_fichier_log, String message);
void Log(String nom_fichier_log, TStringList *lignes);

// -------------------------------------------------------------------------
int main(int argc, char **argv, char **env)
{
String txt, str;
float latitude, longitude;

//Modification 24 mars 2020
float zvent;
double vmes;
vector<double> ustarb;
vector<int> pgtb;

float seuil_vent_calme;
float grille_xmin, grille_ymin, grille_xmax, grille_ymax, grille_resolution, ff;
double x1,x2,x3,x4,xray,xneb,xvit,xdir,vtp,sigv,ustar,teq,xeq,sigy,sigy0,sigtheta,sigtheta0,vitmod,rcalme,xminr,xmaxr,yminr,ymaxr;
double xb1,xb2,yb1,yb2,ex,rayonf,rayon1,rayon2,dist,rapp1,uu,vv,denom,num,em,ecart,ecart1,ecart2,z0;
int f, nn,i,c,ijul,j,ip,nmeteo,iok,idp,i1,i2,indt,retour,ipp,ix,iy,fdebug;
vector<double> vit;      // Vitesse du vent a 10 m (m/s) dans le fichier meteo
vector<double> dir;      // Direction du vent a 10 m (°) dans le fichier meteo
vector<double> ray;      // Rayonnement (W/m2) dans le fichier meteo
vector<double> neb;      // Nebulosite (octa) dans le fichier meteo
vector<__int64> temps;  // Nombre de secondes depuis 1/1/1990 correspondant aux differentes donnees meteo
vector<String> Pid;  // Identifiant de chaque plainte
vector<double> Px,Py;   // Coordonnées de chaque plainte
vector<__int64> Ptemps;    // Nombre de secondes depuis 1/1/1990 correspondant au début de chaque plainte
vector<double> Pdur;    // Duree en secondes de chaque plainte
vector<String> Sid;  // Identifiant de chaque source potentielle
vector<double> Sx,Sy;   // Coordonnées de chaque source potentielle
vector<String> vpars;
vector<int> Pts,PtsOK;
vector<double> PtsX,PtsY;
TDateTime date0,date1;
FILE* fichier = NULL;
char ligne[10000]="";
int nx,ny,idpmax,indinf,indsup;
vector<double> xb,yb,zb,rb,vitb,lmob,hb;
vector<String> vars;
vector<__int64> tb;
__int64 tmin,tmax,i64;
double **meteo = NULL;   // Valeur des parametres meteo interpolles sur chaque minutes necessaires au traitement de toutes les plaintes : meteo[itemps][iparam] : iparam : 0:Vitesse, 1:direction, 2:Nebul, 3:Rayonnement
int **DurInterp = NULL ;   // Duree en secondes separant la valeur du parametre meteo de la valeur non interpolee le plus proche
int *MaxDurInterp = NULL ;  // Max de DurInterp pour chacune des plaintes
int **retro_grille = NULL ;    // retro_grille[i][j] : Nbre d'impact d'une des retrodiffusions associées a une plainte sur le point de (i,j) de la grille. (1 plainte peut etre associee a plusieurs retro : 1 par quart d'heure de la duree de la plainte)
float **retro_grille_all = NULL ;    // retro_grille_all[i][j] : Nbre d'impact d'une des retrodiffusions associées a TOUTES les plaintes sur le point de (i,j) de la grille
int **retro_source = NULL ;  // retro_source[i][j] : Nbre d'impact de l'une des retrodiffusions associées a la plainte j (Note : premiere plainte=plainte 0) sur la source i (Note : premiere source=source 0)
int *flagSource=NULL;
int *flagProbMTOPlainte=NULL; // Vaut 1 si la plainte n'a pas pu etre traitée du fait de manque de donnees météo
int annee, mois, jour, heure, minute, seconde;
TStringList *strs;

TStringList *temp_sortie = new TStringList();

// Resolution temporelle des plaintes (en secondes) : Par exemple si dtplainte=15 minutes, pour une plainte ayant duréé 1 heure, 4 rétrotrajectoire seront calculées
double dtplainte = 15 * 60;
// Resolution temporelle en seconde utilisee pour le calcul des retrodiffusion
double deltat = 60;
// Temps d'integration des mesures meteo en heure
double tinteg = 0.25;
// Altitude des bouffées en m
// Fixe notamment l'altitude du vent qui sera utilisée pour la retrodiffusion.
// 2 m peut être choisi car c'est l'altitude de perception de la gene
// Mais ce n'est pas necessairement l'altitude des sources
// Au final on choisit 10 m pour etre un peu majorant en terme de distance parcourue
double zbouf = 10;
// Durée max de rétrotrajectoire
double duree_max = 7200;

bool erreur = false;
bool erreur2 = false;

String RepertoireSite, CodeSite, NomFichierLog, IdRequete, RepSortie;
TDateTime DateDebut, DateDebut0H, DateFin;

// Modification 24 mars 2020
double fcor;


DateDebut = EncodeDateTime(1990,1,1,0,0,0,0);
DateFin = EncodeDateTime(1990,1,1,0,0,0,0);

if(argc>3)
{
	CodeSite = ToString(argv[1]);
	date0 = EncodeDateTime(1990,1,1,0,0,0,0);
	IdRequete = ToString(argv[2]);
	for(i=3;i<argc;++i)
	{
		explode(&vars, ToString(argv[i]), ':');
		if(vars.size()==5)
		{
			float px0, py0, pdur0;
			if((!valueOf(SubString(vars[1],1,4), annee))||(!valueOf(SubString(vars[1],5,2), mois))||(!valueOf(SubString(vars[1],7,2), jour))||(!valueOf(SubString(vars[1],9,2), heure))||(!valueOf(SubString(vars[1],11,2), minute))||(!valueOf(vars[2], px0))||(!valueOf(vars[3], py0))||(!valueOf(vars[4], pdur0)))
			{
				cout << "ERREUR : arguments attendus : [code_site] [id_unique_identifiant_la_requête] [plainte n°1, avec id_plainte:date_au_format_yyyymmddhhnn:X:Y:duree_en_minutes] [plainte n°2 , avec id_plainte:date_au_format_yyyymmddhhnn:X:Y:duree_en_minutes] ..." << endl;
				return 0;
			}
			else
			{
				date1 = EncodeDateTime(annee, mois, jour, heure, minute, 0, 1);
				if(DateDebut.FormatString("yyyy")=="1990")
					DateDebut = date1;
				else if(ComparerDates(date1, DateDebut)<0)
					DateDebut = date1;
				if(ComparerDates(date1, DateFin)>0)
					DateFin = date1;
				Ptemps.push_back(SecondsBetween(date0, date1));
				Pid.push_back(ToString(vars[0]));
				Px.push_back(px0);
				Py.push_back(py0);
				Pdur.push_back(pdur0/60.);
			}
		}
		else
		{
			cout << "ERREUR : arguments attendus : [code_site] [id_unique_identifiant_la_requête] [plainte n°1, avec id_plainte:date_au_format_yyyymmddhhnn:X:Y:duree_en_minutes] [plainte n°2 , avec id_plainte:date_au_format_yyyymmddhhnn:X:Y:duree_en_minutes] ..." << endl;
			return 0;
		}
	}
}
else
{
	cout << "ERREUR : arguments attendus : [code_site] [id_unique_identifiant_la_requête] [plainte n°1, avec id_plainte:date_au_format_yyyymmddhhnn:X:Y:duree_en_minutes] [plainte n°2 , avec id_plainte:date_au_format_yyyymmddhhnn:X:Y:duree_en_minutes] ..." << endl;
	return 0;
}

// Création du répertoire de sortie
RepertoireSite = "";
explode(&vpars, ToString(argv[0]), '/');
for(c=0;c<vpars.size()-1;++c)
	RepertoireSite = RepertoireSite + "/" + vpars[c];
RepertoireSite = RepertoireSite + "/site/"+CodeSite+"/";
if(!erreur)
{
	RepSortie = RepertoireSite+"data/plaintes/";
	ForceDirectories(RepSortie);
	if(!DirectoryExists(RepSortie))
	{
		Log(NomFichierLog, "Impossible de créer le répertoire de sortie : "+RepSortie);
		return 0;
	}
	else
		NomFichierLog = RepSortie+IdRequete+".log"; 
}

Log(NomFichierLog, Now().FormatString("dd/mm/yyyy hh:nn:ss")+ " DEBUT CALCUL ODOREVERSE DU "+DateDebut.FormatString("dd/mm/yyyy hh:nn")+" AU "+DateFin.FormatString("dd/mm/yyyy hh:nn")+" POUR LE SITE \""+CodeSite+"\"");

Log(NomFichierLog, "Liste des plaintes a traiter :");
for(ip=0;ip<Pid.size();ip++)
	Log(NomFichierLog, "  "+Pid[ip]+" : ptemps="+ToString(Ptemps[ip])+" x="+ToString(Px[ip])+" y="+ToString(Py[ip])+" duree="+ToString(Pdur[ip]));

// Lecture du .ini du site
latitude = -999;
longitude = -999;
seuil_vent_calme = -999;
grille_xmin = -999;
grille_ymin = -999;
grille_xmax = -999;
grille_ymax = -999;
z0 = -999;
if(FileExists(RepertoireSite+"data/web_site.ini"))
{
	TStringList *fichier = new TStringList();
	fichier->LoadFromFile(RepertoireSite+"data/web_site.ini");
	for(f=0;f<fichier->Count;++f)
	{
		//Log(NomFichierLog, SubString(fichier->Strings[f], 1, 26) );
		if(SubString(fichier->Strings[f], 1, 17)=="Station_latitude=")
			valueOf(SubString(fichier->Strings[f], 18, 20), latitude);
		else if(SubString(fichier->Strings[f], 1, 18)=="Station_longitude=")
			valueOf(SubString(fichier->Strings[f], 19, 20), longitude);
		else if(SubString(fichier->Strings[f], 1, 9)=="Rugosite=")
			valueOf(SubString(fichier->Strings[f], 10, 20), z0);
		else if(SubString(fichier->Strings[f], 1, 9)=="Altitude=")
			valueOf(SubString(fichier->Strings[f], 10, 20), zvent);
        else if(SubString(fichier->Strings[f], 1, 18)=="VitesseMinCapteur=")
            valueOf(SubString(fichier->Strings[f], 19, 20), seuil_vent_calme); 
		else if(SubString(fichier->Strings[f], 1, 16)=="Odoreverse_Xmin=")
			valueOf(SubString(fichier->Strings[f], 17, 20), grille_xmin);
		else if(SubString(fichier->Strings[f], 1, 16)=="Odoreverse_Xmax=")
			valueOf(SubString(fichier->Strings[f], 17, 20), grille_xmax);
		else if(SubString(fichier->Strings[f], 1, 16)=="Odoreverse_Ymin=")
			valueOf(SubString(fichier->Strings[f], 17, 20), grille_ymin);
		else if(SubString(fichier->Strings[f], 1, 16)=="Odoreverse_Ymax=")
			valueOf(SubString(fichier->Strings[f], 17, 20), grille_ymax);
	}
	delete fichier;		
	//printf("hauteurStation=%f\n",zvent);
	// Modification 24 mars 2020
	// Normalement il ne faudrait pas faire ça mais on le fait pour etre 
	//   conforme a ce qui est fait dans ADMS
	zvent=zvent+z0;
	// Fin modif
}
else
{
	Log(NomFichierLog, "ERREUR, fichier "+RepertoireSite+"data/web_site.ini inexistant.");
	return 0;
}

if((latitude<-998)||(longitude<-998))
{
	Log(NomFichierLog, "ERREUR : [Meteo]Station_latitude et/ou [Meteo]Station_longitude non défini.");
    erreur = true;
}

// Modification 24 mars 2020
// Calcul du facteur de coriolis
fcor=((4.*3.141592)/86400.)*sin((2.*3.141592*latitude)/360.);
// Fin modification 24 mars 2020

if(z0<-998)
{        
	Log(NomFichierLog, "ERREUR : [Site]Rugosite non défini.");
    erreur = true;
}

if(seuil_vent_calme<0.1)
{        
	Log(NomFichierLog, "ERREUR : [Vents faibles]SeuilVitesse ou [Vents faibles]SeuilVentFaible non défini, et le seuil minimal des vents calmes doit être >=0.1.");
    erreur = true;
}

if(((int)grille_xmin==-999)||((int)grille_ymin==-999)||((int)grille_xmax==-999)||((int)grille_ymax==-999))
{
	Log(NomFichierLog, "ERREUR : [Domaine]Odoreverse_Xmin ou [Domaine]Odoreverse_Xmax ou [Domaine]Odoreverse_Ymin ou [Domaine]Odoreverse_Ymax non défini.");
    erreur = true;
}

if(!erreur)
{
	if(grille_xmin>grille_xmax)
	{
		Log(NomFichierLog, "ERREUR : la coordonnée XMin est supérieure à la coordonnée XMax.");
		erreur = true;
	}
	else if(grille_ymin>grille_ymax)
	{
		Log(NomFichierLog, "ERREUR : la coordonnée YMin est supérieure à la coordonnée YMax.");
		erreur = true;
	}
	else
	{
		// On calcule la résolution de la grille de façon à ce qu'il y ait environ 400 * 400 points
		grille_resolution = (grille_xmax-grille_xmin)/400;
		Log(NomFichierLog, "Résolution de la grille fixée à "+ToString(grille_resolution)+" mètres.");
	}
}

// Lecture du fichier sources
if(!erreur)
{
	if(FileExists(RepertoireSite+"data/sources.ini"))
	{
		TStringList *fichier = new TStringList();
		fichier->LoadFromFile(RepertoireSite+"data/sources.ini");
		int col_x, col_y, col_nom, nbcols;
		col_x = -1;
		col_y = -1;
		col_nom = -1;
		if(fichier->Count>1)
		{
			explode(&vpars, fichier->Strings[0], ';');
			nbcols = vpars.size();
			for(c=0;c<nbcols;++c)
			{
				if(vpars[c]=="x")
					col_x = c;
				else if(vpars[c]=="y")
					col_y = c;
				else if(vpars[c]=="nom")
					col_nom = c;
			}
		}
		if((col_x!=-1)&&(col_y!=-1)&&(col_nom!=-1))
		{
			float sx0, sy0;
			for(int l=1;l<fichier->Count;++l)
			{
				explode(&vpars, fichier->Strings[l], ';');
				if(vpars.size()==nbcols)
				{
					if((valueOf(vpars[col_x],sx0))&&(valueOf(vpars[col_y],sy0)))
					{
						Sid.push_back(vpars[col_nom]);
                        Sx.push_back(sx0);
                        Sy.push_back(sy0);
						if((sx0<grille_xmin)||(sx0>grille_xmax)||(sy0<grille_ymin)||(sy0>grille_ymax))
                        {
							Log(NomFichierLog, "ERREUR : la source \""+vpars[col_nom]+"\" est en dehors du domaine, impossible de continuer.");
                            erreur=true;
                        }
					}
				}
			}                                             
		}
		delete fichier;
		
		if(Sid.size()==0)
		{
			Log(NomFichierLog, "ERREUR : aucune source n'a été identifiée dans "+RepertoireSite+"data/sources.ini");
            erreur=true;
		}
	}
	else
	{
		Log(NomFichierLog, "ERREUR, fichier "+RepertoireSite+"data/sources.ini inexistant.");
		return 0;
	}
}

if(!erreur)
{
	TStringList *fichier_meteo = NULL;
	TDateTime date_meteo_debut, date_meteo_fin, date_meteo_temp;
	__int64 date_meteo_fin64, date_temp64;
	int col_heure, col_minute, col_dir, col_vit, col_ray, col_neb, nbcols, l, j, h, m;
	String nom_fichier_meteo;
	float vitesse, direction, nebulosite, rayonnement;
	vector<String> vars;
	int nb_jours = DaysBetween(DateDebut, DateFin);
	date0 = EncodeDateTime(1990,1,1,0,0,0,0);
        DateDebut0H=EncodeDateTime(DateDebut.year(),DateDebut.month(),DateDebut.day(),0,0,0,0);
	for(j=-2;j<=nb_jours;j++)
	{
		//date_meteo_temp = IncDay(DateDebut, j);
		date_meteo_temp = IncDay(DateDebut0H, j);
		i64 = SecondsBetween(date0, date_meteo_temp);
		for(h=0;h<=23;++h)
			for(m=0;m<=59;++m)
			{
				temps.push_back(i64+m*60+h*60*60);
				vit.push_back (-999);
				dir.push_back (-999);
				neb.push_back (-999);
				ray.push_back (-999);
			}
	}
	
	for(j=-2;j<=nb_jours;j++)
	{
		date_meteo_temp = IncDay(DateDebut, j);
		Log(NomFichierLog, "Lecture de la météo du "+date_meteo_temp.FormatString("dd/mm/yyyy")+"...");
		nom_fichier_meteo = RepertoireSite+"data/"+date_meteo_temp.FormatString("yyyy")+"/"+date_meteo_temp.FormatString("yyyy_mm_dd")+"/meteo.txt";
		if(FileExists(nom_fichier_meteo))
		{
			Log(NomFichierLog, "   Fichier lu : "+nom_fichier_meteo);
			col_heure = -1;
			col_minute = -1;
			col_dir = -1;
			col_vit = -1;
			col_ray = -1;
			col_neb = -1;
			fichier_meteo = new TStringList();
			fichier_meteo->LoadFromFile(nom_fichier_meteo);
			if(fichier_meteo->Count>1)
			{
				explode(&vars, fichier_meteo->Strings[0], ';');
				nbcols = vars.size();
				for(c=0;c<vars.size();++c)
				{
					if(vars[c]=="heure")
						col_heure = c;
					else if(vars[c]=="minute")
						col_minute = c;
					else if(vars[c]=="direction")
						col_dir = c;
					else if(vars[c]=="vitesse")
						col_vit = c;
					else if(vars[c]=="nebulosite")
						col_neb = c;
					else if(vars[c]=="rayonnement")
						col_ray = c;
				}
			}
			if((col_heure!=-1)&&(col_minute!=-1))
			{
				for(l=1;l<fichier_meteo->Count;++l)
				{
					explode(&vars, fichier_meteo->Strings[l], ';');
					if((vars.size()==nbcols)&&(valueOf(vars[col_heure], heure))&&(valueOf(vars[col_minute], minute)))
					{
						i = (j+2)*24*60+heure*60+minute;
						vitesse = -999;								
						if(col_vit!=-1)
							if(!valueOf(vars[col_vit], vitesse))
								vitesse = -999;
						direction = -999;
						if(col_dir!=-1)
							if(!valueOf(vars[col_dir], direction))
								direction = -999;
						nebulosite = -999;
						if(col_neb!=-1)
							if(!valueOf(vars[col_neb], nebulosite))
								nebulosite = -999;
						rayonnement = -999;
						if(col_ray!=-1)
							if(!valueOf(vars[col_ray], rayonnement))
								rayonnement = -999;
						vit[i] = vitesse;

						// Modification 24 mars 2020
						//dir[i] = direction;
						if(vitesse >= (seuil_vent_calme-0.01)) 
						    dir[i] = direction;
						else
						    dir[i] = -999;
						//printf("i=%d,dir=%e\n",i,dir[i]);

						neb[i] = nebulosite;
						ray[i] = rayonnement;
						//if(i==2880) {
						//   date1 = IncSecond(date0 , temps[i]);
						//   printf("TEST dir[i]=%e,vit[i]=%e,temps[i]=%d,heure=%d,minute=%d,neb=%f\n",dir[i],vit[i],temps[i],heure,minute,neb[i]);
						//   Log(NomFichierLog, "TEST " + date1.FormatString("dd/mm/yyyy hh:nn:ss")+ " date0= "+date0.FormatString("dd/mm/yyyy hh:nn:ss"));
						//}
					}
				}
			}
			delete fichier_meteo;
			fichier_meteo = NULL;
		}
		nom_fichier_meteo = RepertoireSite+"data/"+date_meteo_temp.FormatString("yyyy")+"/"+date_meteo_temp.FormatString("yyyy_mm_dd")+"/meteomin.txt";		
		if(FileExists(nom_fichier_meteo))
		{
			Log(NomFichierLog, "   Fichier lu : "+nom_fichier_meteo);
			col_heure = -1;
			col_minute = -1;
			col_dir = -1;
			col_vit = -1;
			col_ray = -1;
			col_neb = -1;
			fichier_meteo = new TStringList();
			fichier_meteo->LoadFromFile(nom_fichier_meteo);
			if(fichier_meteo->Count>1)
			{
				explode(&vars, fichier_meteo->Strings[0], ';');
				nbcols = vars.size();
				for(c=0;c<vars.size();++c)
				{
					if(vars[c]=="heure")
						col_heure = c;
					else if(vars[c]=="minute")
						col_minute = c;
					else if(vars[c]=="direction")
						col_dir = c;
					else if(vars[c]=="vitesse")
						col_vit = c;
					else if(vars[c]=="nebulosite")
						col_neb = c;
					else if(vars[c]=="rayonnement")
						col_ray = c;
				}
			}
			if((col_heure!=-1)&&(col_minute!=-1))
			{
				for(l=1;l<fichier_meteo->Count;++l)
				{
					explode(&vars, fichier_meteo->Strings[l], ';');
					if((vars.size()==nbcols)&&(valueOf(vars[col_heure], heure))&&(valueOf(vars[col_minute], minute)))
					{
						i = (j+2)*24*60+heure*60+minute;
						vitesse = -999;								
						if(col_vit!=-1)
							if(!valueOf(vars[col_vit], vitesse))
								vitesse = -999;
						direction = -999;
						if(col_dir!=-1)
							if(!valueOf(vars[col_dir], direction))
								direction = -999;
						nebulosite = -999;
						if(col_neb!=-1)
							if(!valueOf(vars[col_neb], nebulosite))
								nebulosite = -999;
						rayonnement = -999;
						if(col_ray!=-1)
							if(!valueOf(vars[col_ray], rayonnement))
								rayonnement = -999;
						vit[i] = vitesse;

						// Modification 24 mars 2020
						//dir[i] = direction;
						if(vitesse >= (seuil_vent_calme-0.01)) 
						    dir[i] = direction;
						else
						    dir[i] = -999;

						neb[i] = nebulosite;
						ray[i] = rayonnement;
					}
				}
			}
			delete fichier_meteo;
			fichier_meteo = NULL;
		}		
	}

	 // Calcul nebulosite à partir du rayonnement pour les cas ou la nebul n'est pas disponible
	 Log(NomFichierLog, "Calcul de la nébulosité à partir du rayonnement si nécessaire...");
	 for(i=0;i<neb.size();i++)
	 {
		if(neb[i] < -998)
		{
			date1 = IncSecond(date0 , temps[i]);
			ijul = DayOfTheYear(date1);
			// Il faut bien prendre en compte les secondes memes si les données ne sont qu'a la minute près du fait de l'imprecision numerique introduite par les formats date
			valueOf(date1.FormatString("hh"), heure);
			valueOf(date1.FormatString("nn"), minute);
			valueOf(date1.FormatString("ss"), seconde);
			
			x1 = heure+minute/60.+seconde/3600.;
			x2=(double)AngleSol(latitude,longitude,x1,ijul);
			if(x2 <= 0)
				neb[i]=5;  // Hypothese car la nebul ne peut pas etre calculee de nuit a partir du rayonnement
			else if(ray[i] > -998)
			{
				x3=ray[i];
				x4=((double)990 * sin(x2))-(double)30;
				if(x4 > 2e-6)
				{
					if(ray[i] > x4 - 1e-6)
						x3 = x4 - 1e-6;
					x4=((double)1-(x3/x4))/0.75;
					x4=pow(x4,((double)1/3.4));
					neb[i]=(double)8*x4;
					if(neb[i] > 8)
						neb[i]=(double)8;
					if(neb[i] < 0)
						neb[i]=(double)0;  // Normalement ce cas ne peut pas se produire mais bon ...
				}
				else	// Rayonnement très faible ou negatif => nuit => neb=5
					neb[i]=5;
			}
		}
	 }

	 // Calcul rayonnement à partir de la nebul pour les cas ou le rayonnement n'est pas dispo
	 Log(NomFichierLog, "Calcul du rayonnement à partir de la nébulosité si nécessaire...");
	 for(i=0;i<neb.size();i++)
	 {
		if(ray[i] < -998)
		{
			date1 = IncSecond(date0 , temps[i]);
			ijul = DayOfTheYear(date1);
			valueOf(date1.FormatString("hh"), heure);
			valueOf(date1.FormatString("nn"), minute);
			valueOf(date1.FormatString("ss"), seconde);
			// Il faut bien prendre en compte les secondes memes si les données ne sont qu'a la minute près du fait de l'imprecision numerique introduite par les formats date
			x1 = heure+minute/60.+seconde/3600.;
			x2=(double)AngleSol(latitude,longitude,x1,ijul);
			if(x2 <= 0)
				ray[i]=0;     // C'est la nuit : Rayonnement nul quelle que soit la nebulosite
			else if(neb[i] > -998)
			{
				x3=(neb[i]/(double)8);
				x3=pow(x3,3.4);
				ray[i]=((double)1-0.75*x3)*((990 * sin(x2))-30);
				if(ray[i] < 0)
					ray[i]=(double)0;
			}
		}
	}
}

if(!erreur)
{
	// Allocation matrice
	Log(NomFichierLog, "Initialisation de la matrice...");
	nx=1+round((grille_xmax-grille_xmin)/grille_resolution);
	ny=1+round((grille_ymax-grille_ymin)/grille_resolution);

	retro_grille = new int *[nx];
	for (i=0;i<nx;i++)
		retro_grille[i]=new int[ny];

	retro_grille_all = new float *[nx];
	for (i=0;i<nx;i++)
		retro_grille_all[i]=new float[ny];

	retro_source = new int *[Sid.size()];
	for (i=0;i<Sid.size();i++)
		retro_source[i]=new int[Pid.size()];

	flagSource = new int[Sx.size()];
	MaxDurInterp = new int[Pid.size()];
	flagProbMTOPlainte = new int[Pid.size()];

	// Initialisation matrice
	for (i=0;i<nx;i++)
		for (j=0;j<ny;j++)
				retro_grille_all[i][j]=0;

	for (i=0;i<Sid.size();i++)
		for (j=0;j<Pid.size();j++)
				retro_source[i][j]=0;
		   
	for (i=0;i < Pid.size();i++)
		flagProbMTOPlainte[i] = 0;
	
	// Calcul "retrodiffusion" pour chaque plainte		
	for(ip=0;ip<Ptemps.size();ip++)
	{
		Log(NomFichierLog, "Calcul de la rétrodiffusion pour la plainte " + ToString(ip+1)+ "/"+ ToString(Ptemps.size())+ "...");
	
		erreur2 = false;
		// Initialisation matrice
		for (i=0;i<nx;i++)
			for (j=0;j<ny;j++)
				retro_grille[i][j]=0;

		Log(NomFichierLog, "Météo avant interpolation :");

		strs = new TStringList();
		strs->Add("");
		strs->Add("i;temps;vit;dir;neb;ray");
		for(i=0;i<temps.size();++i)
			strs->Add(ToString(i)+";"+ToString(temps[i])+";"+ToString(vit[i])+";"+ToString(dir[i])+";"+ToString(neb[i])+";"+ToString(ray[i]));
		strs->Add("");
		Log(NomFichierLog, strs);
		delete strs;
		
		Log(NomFichierLog, "Interpolation des données météo...");

		// Interpolation des données météo pour chaque minute de la periode necessaire a la plainte
		// Recherche de la periode necessaire a la plainte
		idpmax = 1+floor(((Pdur[ip] - 1e-6)*(double)3600) / dtplainte); // Nbre de retrodiffusion a realiser pour cette plainte
		tmax=Ptemps[ip]+(idpmax-1)*dtplainte;
		tmin=Ptemps[ip]-(__int64)round(duree_max);

		if(meteo!=NULL)
		{
			for(i=0;i<nmeteo;i++)
				delete[] meteo[i];
			delete[] meteo;
		}
		if(DurInterp!=NULL)
		{
			for(i=0;i<nmeteo;i++)
				delete[] DurInterp[i];
			delete[] DurInterp;
		}

		x1=tmax-tmin;
		nn=round((x1/(double)60)); // Nombre de minute dans la periode determinee
		nmeteo=nn+1; // Taille necessaire pour les matrices
		meteo = new double *[nmeteo];
		DurInterp = new int *[nmeteo];
		for (i=0;i<nmeteo;i++)
		{
			meteo[i]=new double[4];  // 4 => Vit, direction, Neb et Rayonnement
			DurInterp[i]=new int[4];  // 4 => Vit, direction, Neb et Rayonnement
			for(j=0;j<4;j++)
			{
				meteo[i][j]=(double)-999;
				DurInterp[i][j]=-999;
			}
		}

		// Interpolation vitesse
		i1=-1;
		for(j=0;j<nmeteo;j++)
		{
			i64=tmin + ((__int64)60 * (__int64)j);
			i=0;
			while( (temps[i] <= i64) & (i < temps.size()) )
			{
				if(vit[i] > -998)
					i1 = i;
				i++;
			}
			if(i1 >= 0)
			{
				if(i64==temps[i1])
				{
					meteo[j][0]=vit[i1];
					DurInterp[j][0]=0;
				} 
				else
				{
					i2=i1+1;
					if(i2 < temps.size())
					{
						while( (i64 >= temps[i2]) || ((vit[i2] < -998) & (i2 < temps.size())) )
						//while( (vit[i2] < -998) & (i2 < temps.size()) )
							i2++;
						if(i2 < temps.size())
						{
							ff=(double)(i64-temps[i1])/(double)(temps[i2]-temps[i1]);
							meteo[j][0]=(1-ff)*vit[i1]+ff*vit[i2];
							DurInterp[j][0]=min(i64-temps[i1],temps[i2]-i64);
						}
					}
				}
			}
		}


		// Interpolation nebulosite
		i1=-1;
		for(j=0;j<nmeteo;j++)
		{
			i64=tmin + ((__int64)60 * (__int64)j);
			i=0;
			while( (temps[i] <= i64) & (i < temps.size()) )
			{
				if(neb[i] > -998) 
					i1 = i;
				i++;
			}
			if(i1 >= 0)
			{
				if(i64==temps[i1])
				{
					meteo[j][2]=neb[i1];
					DurInterp[j][2]=0;
				}
				else
				{
					i2=i1+1;
					if(i2 < temps.size())
					{
						while( (i64 >= temps[i2]) || ((neb[i2] < -998) & (i2 < temps.size())) )
						//while( (neb[i2] < -998) & (i2 < temps.size()) )
							i2++;
						if(i2 < temps.size())
						{
							ff=(double)(i64-temps[i1])/(double)(temps[i2]-temps[i1]);
							meteo[j][2]=(1-ff)*neb[i1]+ff*neb[i2];
							DurInterp[j][2]=min(i64-temps[i1],temps[i2]-i64);
						}
					}
				}
			}
		}

		// Interpolation rayonnement
		i1=-1;
		for(j=0;j<nmeteo;j++)
		{
			i64=tmin + ((__int64)60 * (__int64)j);
			i=0;
			while( (temps[i] <= i64) & (i < temps.size()) )
			{
				if(ray[i] > -998)
					i1 = i;
				i++;
			}
			if(i1 >= 0)
			{
				if(i64==temps[i1])
				{
					meteo[j][3]=ray[i1];
					DurInterp[j][3]=0;
				}
				else
				{
					i2 = i1+1;
					if(i2 < temps.size())
					{
						while( (i64 >= temps[i2]) || ((ray[i2] < -998) & (i2 < temps.size())) )
						//while( (ray[i2] < -998) & (i2 < temps.size()) )
							i2++;
						if(i2 < temps.size())
						{
							ff=(double)(i64-temps[i1])/(double)(temps[i2]-temps[i1]);
							meteo[j][3]=(1-ff)*ray[i1]+ff*ray[i2];
							DurInterp[j][3]=min(i64-temps[i1],temps[i2]-i64);
						}
					}
				}
			}
		}


		// Interpolation direction
		// Modification 24 mars 2020
		//i1=-1;
		for(j=0;j<nmeteo;j++)
		{
			i64=tmin + ((__int64)60 * (__int64)j);
			i=0;
			// Modification 24 mars 2020
			i1=-1; 
			// Premiere recherche en ignorant les cas de vents faibles
			// De cette facon, en cas de vent faible, on utilise la derniere direction correcte
			while( (temps[i] <= i64) & (i < temps.size()) )
			{
				if(dir[i] > -998)
					i1 = i;
				i++;
			}

			// Modification 24 mars 2020
			if(i1 < 0)
			{
                                i=0;
			        // Par contre il se peut qu'on ait que des vents calmes de disponibles 
			        // Dans ce cas, on considere le dernier vent calme et on procedera comme suit :
			        // Si i1 et i2 correspondent a des vents calme, on fixe la direction 
			        //     a 360 (Elle ne sera pas utilisée de toute façon) 
			        // Si seul i1 correspond a un vent calme, on fixe la direction de i1 a celle de i2
				while( (temps[i] <= i64) & (i < temps.size()) )
				{
					if(vit[i] < (seuil_vent_calme-0.01) && vit[i] > -998)
						i1 = i;
					i++;
				}
			}

			if(i1 >= 0)
			{
				if(i64==temps[i1])
				{
					// Modification 24 mars 2020
					//meteo[j][1]=dir[i1];
					if(vit[i1] < (seuil_vent_calme-0.01))
					    meteo[j][1]=-360;  // Ne sera pas utilisé de toute facon
					else
					    meteo[j][1]=dir[i1];
					DurInterp[j][1]=0;
					//printf("j=%d,i1=%d\n",j,i1);
				}
				else
				{
					i2=i1+1;
					if(i2 < temps.size())
					{
						// Premiere tentative de recherche de i2 en ignorant les vents calmes
						// De cette facon, en cas de vent faible, on utilise la premiere direction correcte
						while( (i64 >= temps[i2]) || ((dir[i2] < -998) & (i2 < temps.size())) )
							i2++;
						// Modification 24 mars 2020
						if(i2 >= temps.size())
						{
                                                        i2=i1+1;
							// Par contre il se peut qu'on ait que des vents calmes de disponibles 
							// Dans ce cas, on considere le premier vent calme et on procedera comme suit :
							// Si i1 et i2 correspondent a des vents calme, on fixe la direction 
							//      a 360 (Elle ne sera pas utilisée de toute façon) 
							// Si seul i1 correspond a un vent calme, on fixe la direction de i2 a celle de i1
							while( i2 < temps.size() )
							{
                                                                //printf("dir[i2]=%e,v=%e\n",dir[i2],vit[i2]);
                                                                if(i64 <= temps[i2])
                                                                {
								   if(dir[i2] < -998 && (vit[i2] < (seuil_vent_calme-0.01) && vit[i2] > -998)) break; 
                                                                }
								i2++;
							}
						}
                                                //printf("j=%d,i1=%d,ray[i1]=%e,i2=%d,temps.size=%d\n",j,i1,ray[i1],i2,temps.size());
						if(i2 < temps.size())
						{
							// Modification 24 mars 2020
							if(vit[i1] < (seuil_vent_calme-0.01) && vit[i2] < (seuil_vent_calme-0.01))
							{
								meteo[j][1]=360;  // Ne sera pas utilisé de toute facon
							} 
							else if(vit[i1] < (seuil_vent_calme-0.01))
							{
								meteo[j][1]=dir[i2];
							}
							else if(vit[i2] < (seuil_vent_calme-0.01))
							{
								meteo[j][1]=dir[i1];
							}
							else
							{
								// Bascule dans le sens induisant le minimum de variabilité de la direction
								ff=(double)(i64-temps[i1])/(double)(temps[i2]-temps[i1]);
								x1=dir[i2]-dir[i1];
								if(x1 > 180)
									x1=x1-(double)360;
								if(x1 < -180)
									x1=(double)360+x1;
								meteo[j][1]=dir[i1]+ff*x1;
								if(meteo[j][1] > 360)
									meteo[j][1]=meteo[j][1]-(double)360;
								if(meteo[j][1] < 0)
									meteo[j][1]=meteo[j][1]+(double)360;
							   DurInterp[j][1]=min(i64-temps[i1],temps[i2]-i64);
							}
						}
					}
				}
			}
		}

		Log(NomFichierLog, "Météo après interpolation :");
		strs = new TStringList();
		strs->Add("");
		strs->Add("i;vit;dir;neb;ray");
		for(i=0;i<nmeteo;i++)
			strs->Add(ToString(i)+";"+ToString(meteo[i][0])+";"+ToString(meteo[i][1])+";"+ToString(meteo[i][2])+";"+ToString(meteo[i][3]));
		strs->Add("");
		Log(NomFichierLog, strs);
		delete strs;
		
		Log(NomFichierLog, "     Verification des données disponibles...");

		// Verification que les données necessaires sont disponibles
		// A ce stade la premiere et la derniere valeur de chaque donnees du tableau meteo doit etre > a -998
		iok=1;
		for(i=0;i<4;i++)
		{
			if(meteo[0][i] < -998)
				iok = 0;
			if(meteo[nmeteo-1][i] < -998)
				iok = 0;
		}

		if(iok==0)
		{
			Log(NomFichierLog, "     WARNING : les données météorologiques fournies ne couvrent pas la période nécessaire au traitement de la plainte d'identifiant " + ToString(Pid[ip])+ ". Cette plainte ne sera pas traitée.");
			erreur2=true;
			flagProbMTOPlainte[ip]=1;
		}

		// Vérification que les plaintes sont dans le domaine
		if((Px[ip]<grille_xmin)||(Px[ip]>grille_xmax)||(Py[ip]<grille_ymin)||(Py[ip]>grille_ymax))
		{
			Log(NomFichierLog, "     WARNING : la plainte d'identifiant " + ToString(Pid[ip])+ " est en dehors du domaine, elle ne sera pas traitée.");
			erreur2=true;
			flagProbMTOPlainte[ip]=1;
		}

		// Arrondi de DurInterp pour faire disparaitre l'erreur numerique introduite lors des calculs de date
		for(i=0;i<nmeteo;i++)
			for(j=0;j<4;j++)
				DurInterp[i][j]=60*round(DurInterp[i][j]/60);

		// Calcul du Max de DurInterp
		if(flagProbMTOPlainte[ip]==0)
		{
			MaxDurInterp[ip]=0;
			for(i=0;i<nmeteo;i++)
				for(j=0;j<4;j++)
					if(DurInterp[i][j] > MaxDurInterp[ip])
						MaxDurInterp[ip]=DurInterp[i][j];
		}
		else
			MaxDurInterp[ip]=-999;

		if(erreur2==false)
		{
			Log(NomFichierLog, "     Boucle sur les rétrodiffusions...");
			
			//Decomposition de la plainte en retrotraj de duree egale a dtplainte (dtplainte=15 minutes si rien n'a change)
			idpmax = 1+floor(((Pdur[ip] - 1e-6)*(double)3600) / dtplainte);
			for(idp=0;idp<idpmax;idp++)
			{
				// Init bouffée
				tb.clear();
				xb.clear();
				yb.clear();
				zb.clear();
				rb.clear();
				vitb.clear();
				lmob.clear();
				hb.clear();
				// Modification 24 mars 2020
				pgtb.clear();
				ustarb.clear();

				xb.push_back(Px[ip]);
				yb.push_back(Py[ip]);
				zb.push_back(zbouf);
				rb.push_back(0); // Rayon de la bouffeee
				i64=Ptemps[ip]+(__int64)(idp*dtplainte);
				tb.push_back(i64);
				// Determination indice des param meteo a utiliser
				i64=i64-tmin;
				indt=round((double)i64/(double)60);

				// Modification 24 mars 2020
				xneb=meteo[indt][2];
				xray=meteo[indt][3];
				vmes=meteo[indt][0];
				// Le seuil de vent calme est censé etre le seuil de detection du capteur
				// Quand on est en dessous de ce seuil, la valeur reelle est en fait
				//    comprise entre 0 et ce seuil
				// On utilise donc au final comme vitesse, la valeur de ce seuil / 2.
				if((vmes+0.01) < seuil_vent_calme) vmes=seuil_vent_calme/(double)2.;
				// Calcul de la LMO
				// On utilise 0.23 pour l'albedo, 15°C pour la temperture et 1 pour le parametre de Presley Taylor
				calLmoUstar(z0,0.23,1,vmes,zvent,(double)15,xray,xneb,&x1,&ustar);
				lmob.push_back(x1);
				ustarb.push_back(ustar);
				// Calcul de la hauteur de couche limite
				if(x1 > 0)
				{
					 x2=(4*1.9*0.3*ustar)/(x1*fabs(fcor)); 
					 x2=sqrt(1+x2);
					 x2=(x2-1)*(x1/3.8);
				} else
				{
					 x2=0.3*ustar/fabs(fcor);
				}
				hb.push_back(x2);
				//printf("Hcl=%e\n",x2);
				// Calcul de la vitesse a 10 m 
				xvit=Z2Vit(lmob[lmob.size()-1],hb[hb.size()-1],10.,z0,ustar);  
				// Calcul de la classe de Pasquill
				i1=round(ClassePasquill(xvit, xray, xneb));
				pgtb.push_back(i1);
				// Calcul de la vitesse a l'altitude de la bouffee
				xvit=Z2Vit(lmob[lmob.size()-1],hb[hb.size()-1],zb[zb.size()-1],z0,ustar);  
				vitb.push_back(xvit);
				//vitb.push_back(meteo[indt][0]);
				//i1=round(ClassePasquill(meteo[indt][0], meteo[indt][3], meteo[indt][2]));
				//i2=PGT2LMOH(i1,&x1,&x2);
				//lmob.push_back(x1);
				//hb.push_back(x2);
				// Fin modification 24 mars 2020

				iok=1;
				vtp=deltat;
				while(iok==1)
				{
					i64=Ptemps[ip]+(__int64)(idp*dtplainte)-(__int64)vtp;
					tb.push_back(i64);

					// Determination des parametres meteo au temps considere
					i64=i64-tmin;
					x1=(double)i64/(double)60;
					indt=floor(x1);
					x1=x1-indt;
					if(x1 < 1e-6)
					{
						xvit=meteo[indt][0];
						xdir=meteo[indt][1];
						xneb=meteo[indt][2];
						xray=meteo[indt][3];
					}
					else
					{
						xvit=(1-x1)*meteo[indt][0]+x1*meteo[indt+1][0];
						xneb=(1-x1)*meteo[indt][2]+x1*meteo[indt+1][2];
						xray=(1-x1)*meteo[indt][3]+x1*meteo[indt+1][3];
						// Modification 24 mars 2020
						// Pour la direction, il peut rester des -999 en cas de vents calme
						if(meteo[indt][1] < -998)
						{
							 if(meteo[indt+1][1] < -998)
							 {
								 xdir=360; // Dans ce cas, on est sur que la direction ne sera pas utilisé 
							 } else
							 {
								 xdir=meteo[indt+1][1];
							 }
						} else
						{
							 if(meteo[indt+1][1] < -998)
							 {
								 xdir=meteo[indt][1];
							 } else
							 {
								x2=meteo[indt+1][1]-meteo[indt][1];
								if(x2 > 180)
								   x2=x2-(double)360;
								if(x2 < -180)
								   x2=(double)360+x2;
								xdir=meteo[indt][1]+x1*x2;
								if(xdir > 360)
								   xdir=xdir-(double)360;
								if(xdir < 0)
								   xdir=xdir+(double)360;
							 }
						}
						//x2=meteo[indt+1][1]-meteo[indt][1];
						//if(x2 > 180)
						//	x2=x2-(double)360;
						//if(x2 < -180)
						//	x2=(double)360+x2;
						//xdir=meteo[indt][1]+x1*x2;
						//if(xdir > 360)
						//	xdir=xdir-(double)360;
						//if(xdir < 0)
						//	xdir=xdir+(double)360;
					}

					// Modification 24 mars 2020
					// A ce stade xvit est la vitesse à la hauteur de mesure
					vmes=xvit; // vmes=vitesse a l'altitude de la mesure
					// Le seuil de vent calme est censé etre le seuil de detection du capteur
					// Quand on est en dessous de ce seuil, la valeur reelle est en fait
					//    comprise entre 0 et ce seuil
					// On utilise donc au final comme vitesse, la valeur de ce seuil / 2.
					if((vmes+0.01) < seuil_vent_calme) vmes=seuil_vent_calme/(double)2.;
						// Calcul de la LMO
						// On utilise 0.23 pour l'albedo, 15°C pour la temperture et 1 pour le parametre de Presley Taylor
						calLmoUstar(z0,0.23,1,vmes,zvent,(double)15,xray,xneb,&x1,&ustar);
						//printf("z0,xvit,zvent,xray,xneb,lmo,ustar=%e,%e,%e,%e,%e,%e,%e\n",z0,xvit,zvent,xray,xneb,x1,ustar);
					lmob.push_back(x1);
					ustarb.push_back(ustar);
					//printf("lmo=%e\n",x1);
					//printf("fcor=%e\n",fcor);
					// Calcul de la hauteur de couche limite
					if(x1 > 0)
					{
						x2=(4*1.9*0.3*ustar)/(x1*fabs(fcor)); 
						x2=sqrt(1+x2);
						x2=(x2-1)*(x1/3.8);
					} else
					{
						x2=0.3*ustar/fabs(fcor);
					}
					hb.push_back(x2);
					// Calcul de la vitesse a 10 m 
					xvit=Z2Vit(lmob[lmob.size()-1],hb[hb.size()-1],10.,z0,ustar);  
					// Calcul de la classe de Pasquill
					i1=round(ClassePasquill(xvit, xray, xneb));
					pgtb.push_back(i1);
					//printf("ustar=%e\n",calUstar(zvent+0.5, xvit, x1,x2,0.5));
					//printf("PGT=%d,lmo=%f,H=%f\n",i1,x1,x2);
					//lmob.push_back(x1);
					//hb.push_back(x2);
					// Fin modification 24 mars 2020

					// Calcul de la vitesse a l'altitude de la bouffee
					zb.push_back(zbouf);
					//ustar=calUstar((double)10, xvit, lmob[lmob.size()-1],hb[hb.size()-1],z0);
					xvit=Z2Vit(lmob[lmob.size()-1],hb[hb.size()-1],zb[zb.size()-1],z0,ustar);  // A partir de là, xvit est la vitesse a la hauteur de la bouffee
					vitb.push_back(xvit);
					//printf("Vit=%e,Hcl=%e,lmo=%e\n",xvit,x2,x1);

					// Calcul retrotraj
					if(vmes >= seuil_vent_calme)
					{
						x1=xb[xb.size()-1]-deltat*(double)VectoriserU(xdir,xvit);
						x2=yb[yb.size()-1]-deltat*(double)VectoriserV(xdir,xvit);
						xb.push_back(x1);
						yb.push_back(x2);
					}
					else
					{
						x1=xb[xb.size()-1];
						x2=yb[yb.size()-1];
						xb.push_back(x1);
						yb.push_back(x2);
					}


					// Calcul rayon bouffee
					rcalme=0;
					for(i=(int)round(vtp/deltat);i>0;i--)
					{
							// Modification 24 mars 2020
							// Calcul de la vitesse à la hauteur de mesure
					        vmes=Z2Vit(lmob[i],hb[i],zvent,z0,ustarb[i]);  

                                                //printf("vb=%e,vmes=%e\n",vitb[i],vmes);
                                                 
						// Ajustement de la vitesse du vent a utiliser
                                                // Plus necessaire car vmes a ete recalculee a partir et ustar et lmo
                                                //   qui prennent deja en compte le seuil de vent calme
						//if(vitb[i] >= seuil_vent_calme)
						//if(vmes >= seuil_vent_calme)
						//	vitmod=vitb[i];
						//else
                                                //{ 
						//	vitmod=(seuil_vent_calme/(double)2.);
                                                //}

						// Calcul de la distance parcourue equivalente
						if(i==(int)round(vtp/deltat))
                                                {
							xeq = 0;
                                                } else
                                                { 
                                                    x1=coeffA(pgtb[i]);
                                                    x2=x1*pow((z0/0.03),0.2)*pow(tinteg,0.2);
                                                    xeq=pow(sigy0/x2,(double)1./0.9);
                                                }    
                                                //printf("xeq=%e\n",xeq);
                                                 
                                                // Calcul du sigmay
                                                x1=coeffA(pgtb[i]);
						x2=xeq+vitb[i]*deltat;
                                                sigy=x1*pow(x2,0.9);
                                                //printf("x1=%e,x2=%e\n",x1,x2);

                                                // Correction pour le pas d'integration meteo
                                                sigy=sigy*pow(tinteg,0.2);

                                                // Correction pour la rugosité
                                                sigy=sigy*pow((z0/0.03),0.2);

						// Calcul sigmaV
						//ustar=calUstar(zb[i], vitmod, lmob[i],hb[i],z0);
						//if(hb[i]/lmob[i] <= -0.3)
						//	sigv=sigmaV_instable(zb[i],ustar,lmob[i],hb[i],0,1);
						//else if(hb[i]/lmob[i] <= 1)
						//	sigv=sigmaV_neutre(zb[i],ustar,hb[i],0,1);
						//else
						//	sigv=sigmaV_stable(zb[i],ustar,hb[i],z0,0,1);
                                                //printf("zb[i]=%f,ustar=%f,z0=%f,vitb[i]=%e,lmob[i]=%e,hb[i]=%e,sigv=%f\n",zb[i],ustar,z0,vitb[i],lmob[i],hb[i],sigv);

						// Calcul du temps equivalent
						//if(i==(int)round(vtp/deltat))
						//	teq = 0;
						//else if(hb[i]/lmob[i] >= -0.3)
						//{
						//	RRRsigv=sigv;
						//	RRRhcl=hb[i];
						//	RRRustar=ustar;
						//	RRRsigy=sigy0;
						//	teq=zbrent(resol_sigySTAB_ADMS, 0., 1.e6, 0.0001,&retour);
						//}
						//else
						//{
						//	x1=(hb[i]*ustar*ustar*ustar)/(0.4*fabs(lmob[i]));
						//	RRRwstar=pow(x1,(double)1/(double)(3));
						//	x1=0.3*RRRwstar*RRRwstar;
						//	RRRsigvc=sqrt(x1);
						//	RRRsigvn=sqrt((sigv*sigv)-x1);
						//	RRRhcl=hb[i];
						//	RRRustar=ustar;
						//	RRRsigy=sigy0;
					        //	teq=zbrent(resol_sigyCBL_ADMS, 0., 1.e6, 0.0001,&retour);
						//}

						// Calcul du sigy
						//if(hb[i]/lmob[i] >= -0.3)
						//{
						//	x1=teq+deltat;
						//	sigy=cal_sigySTAB_ADMS(x1,sigv,hb[i],ustar);
						//}
						//else
						//{
						//	x1=teq+deltat;
						//	sigy=cal_sigyCBL_ADMS(x1,sigv,hb[i],ustar,lmob[i]);
                                                //        printf("t=%e,sigv=%f,hb[i]=%f,ustar=%f,lmob[i]=%f,sigy=%f\n",x1,sigv,hb[i],ustar,lmob[i],sigy);
						//}
						sigy0=sigy; // Stockage du sigy pour l'iteration suivante

						// Calcul du x equivalent pour sigmatheta (battement de vent)
						//if(i==(int)round(vtp/deltat))
						//	xeq = 0;
						//else
						//{
						//	x1=0.065*sqrt((7*tinteg)/vitmod);
						//	xeq=sigtheta0/x1;
						//}
						// Calcul du sigmatheta
						//x1=xeq+vitmod*deltat;
						//sigtheta=0.065*x1*sqrt((7*tinteg)/vitmod);
						//sigtheta0=sigtheta;  // Stockage du sigtheta pour l'iteration suivante

						// Augmentation de rcalme si on est dans une periode de vent calme
						//if(vitb[i] < seuil_vent_calme)
						if(vmes < seuil_vent_calme)
                                                {
							//rcalme=rcalme+(seuil_vent_calme*deltat);
							rcalme=rcalme+(vitb[i]*deltat);
                                                }
                                                // Fin correction 24 mars 2020
					}
                                        // Modification 24 mars 2020
					// Calcul du rayon correspondant au sigma de dispersion
					// En fait ici on doit multiplier par 4 le sigma pour que la bouffee contienne 95% de la masse de la bouffée
					// Mais on doit aussi diviser par 2 pour avoir le rayon et non pas le diametre
					//x1=sqrt(sigy*sigy+sigtheta*sigtheta);
					x1=2*sigy;
                                        //printf("sigy=%f,sigtheta=%f\n",sigy,sigtheta);

					// Calcul final du rayon
					rb.push_back(x1+rcalme);
                                        //printf("x1=%f,rcalme=%f\n",x1,rcalme);

					// On verifie si la bouffee influence toujours la grille
					if(xb[i]+rb[i] < grille_xmin || xb[i]-rb[i] > grille_xmax)
						if(yb[i]+rb[i] < grille_ymin || yb[i]-rb[i] > grille_ymax)	// Bouffee en dehors de la grille : Arret de la retrotrajectoire
							iok = 0;

					// Incrementation du temps
					vtp=vtp+deltat;
					if(vtp > round(duree_max))
						iok = 0;
				}

				// Recherche des points de grilles affectes

				// Recherche de la zone reduite : recherche de xminr,xmaxr,yminr,ymaxr
				xminr = 1.e30;
				yminr = 1.e30;
				xmaxr = -1.e30;
				ymaxr = -1.e30;
				for(i=0;i<xb.size();i++)
				{
					if(xb[i]+rb[i] > xmaxr)
						xmaxr = xb[i]+rb[i];
					if(xb[i]-rb[i] < xminr)
						xminr = xb[i]-rb[i];
					if(yb[i]+rb[i] > ymaxr)
						ymaxr = yb[i]+rb[i];
					if(yb[i]-rb[i] < yminr)
						yminr = yb[i]-rb[i];
				}

				// Preparation de la selection de point
				// On met tous les points eventuellement recouvert par la retrodiffusion dans une liste
				// Dans la suite, si un point se trouve recouvert par une bouffée, on le sort de la liste ce qui réduit le temps de calcul pour les bouffées suivantes
				Pts.clear();
				PtsX.clear();
				PtsY.clear();
				PtsOK.clear();

				for(i=0;i<nx;i++)
				{
					x1=grille_xmin+(i-1)*grille_resolution;
					for(j=0;j<ny;j++)
					{
						x2 = grille_ymin+(j-1)*grille_resolution;
						if((x1 >= xminr)&&(x2 >= yminr)&&(x1 <= xmaxr)&&(x2 <= ymaxr))
						{
							Pts.push_back(i+j*nx);
							PtsX.push_back(x1);
							PtsY.push_back(x2);
						}
					}
				}

				// Mise a zero du compteur pour les sources
				for (i=0;i<Sid.size();i++)
					flagSource[i] = 0;                               


				// Controle des points recouverts par le panache
				// On commence par les dernières bouffées : Ce sont elles qui risquent de recouvrir le plus de point
				//   donc ce sont aussi celles qui risquent d'éliminer un maximum de points pour les traitements suivants
				for(ipp=xb.size()-2;ipp>=0;ipp--)
				{
					xb1 = xb[ipp];
					yb1 = yb[ipp];
					xb2 = xb[ipp+1];
					yb2 = yb[ipp+1];

					uu = (xb2 - xb1) / deltat;
					vv = (yb2 - yb1) / deltat;
					rayon1 = rb[ipp];
					rayon2 = rb[ipp+1];
					rapp1 = (rayon2 - rayon1) / deltat;
					denom = (2 * uu * uu) + (2 * vv * vv) - (2 * rapp1 * rapp1);

					if(denom > 0)
					{
						//Il est important de parcourir les points en commençant par la fin car des élements vont etre retirer aux vecteurs lors du parcours de la boucle
						for(i=Pts.size()-1;i>=0;i--)
						{
							num = 2 * uu * (PtsX[i] - xb1) + 2 * vv * (PtsY[i] - yb1) + 2 * rayon1 * rapp1;
							// (Temps-temps debut pas de temps) correspondant au minimum dist2 - R2 (pas forcement sur l'intervalle)
							ex = num / denom;
							// Recherche de l'écart (Temps-temps debut pas de temps) correspondant au minimum dist2-R2 sur le pas de temps
							if(ex <= 0)
								em = 0;
							else if(ex <= deltat)
								em = ex;
							else
								em = deltat;

							x1=PtsX[i]-(xb1+(uu*em));
							x2=PtsY[i]-(yb1+(vv*em));
							dist=sqrt(x1*x1+x2*x2);
							rayonf = rayon1+em*rapp1;

							if(dist <= rayonf)
							{
								PtsOK.push_back(Pts[i]);
								Pts.erase(Pts.begin()+i);
								PtsX.erase(PtsX.begin()+i);
								PtsY.erase(PtsY.begin()+i);
							}
						}

						// recherche des sources se trouvant dans la bouffée
						for(i=0;i<Sx.size();i++)
						{
							num = 2 * uu * (Sx[i] - xb1) + 2 * vv * (Sy[i] - yb1) + 2 * rayon1 * rapp1;

							// (Temps-temps debut pas de temps) correspondant au minimum dist2 - R2 (pas forcement sur l'intervalle)
							ex = num / denom;
							// Recherche de l'écart (Temps-temps debut pas de temps) correspondant au minimum dist2-R2 sur le pas de temps
							if(ex <= 0)
								em = 0;
							else if(ex <= deltat)
								em = ex;
							else
								em = deltat;

							x1=Sx[i]-(xb1+(uu*em));
							x2=Sy[i]-(yb1+(vv*em));
							dist=sqrt(x1*x1+x2*x2);
							rayonf=rayon1+em*rapp1;

							if(dist <= rayonf)
								flagSource[i] = 1;
						}
					}
					else
					{
						if(denom < 0)
						{
							//Il est important de parcourir les points en commençant par la fin car des élements vont etre retirer aux vecteurs lors du parcours de la boucle
							for(i=Pts.size()-1;i>=0;i--)
							{
								num = 2 * uu * (PtsX[i] - xb1) + 2 * vv * (PtsY[i] - yb1) + 2 * rayon1 * rapp1;
								// (Temps-temps debut pas de temps) correspondant au minimum dist2 - R2 (pas forcement sur l'intervalle)
								ex = num / denom;
								// Recherche de l'écart (Temps-temps debut pas de temps) correspondant au minimum dist2-R2 sur le pas de temps
								if(ex <= 0)
								{
									em = deltat;
									x1=PtsX[i]-(xb1+(uu*em));
									x2=PtsY[i]-(yb1+(vv*em));
									dist=sqrt(x1*x1+x2*x2);
									rayonf=rayon1+em*rapp1;
									ecart = dist - rayonf;
								}
								else if(ex <= deltat)
								{
									em = deltat;
									x1=PtsX[i]-(xb1+(uu*em));
									x2=PtsY[i]-(yb1+(vv*em));
									dist=sqrt(x1*x1+x2*x2);
									rayonf=rayon1+em*rapp1;
									ecart1 = dist - rayonf;

									em = 0;
									x1=PtsX[i]-(xb1+(uu*em));
									x2=PtsY[i]-(yb1+(vv*em));
									dist=sqrt(x1*x1+x2*x2);
									rayonf=rayon1+em*rapp1;
									ecart2 = dist - rayonf;

									ecart=min(ecart1,ecart2);
								}
								else
								{
									em=0;
									x1=PtsX[i]-(xb1+(uu*em));
									x2=PtsY[i]-(yb1+(vv*em));
									dist=sqrt(x1*x1+x2*x2);
									rayonf=rayon1+em*rapp1;
									ecart = dist - rayonf;
								}

								if(ecart <= 0)
								{
									PtsOK.push_back(Pts[i]);
									Pts.erase(Pts.begin()+i);
									PtsX.erase(PtsX.begin()+i);
									PtsY.erase(PtsY.begin()+i);
								}
							}

							for(i=0;i<Sx.size();i++)
							{
								num = 2 * uu * (Sx[i] - xb1) + 2 * vv * (Sy[i] - yb1) + 2 * rayon1 * rapp1;
								// (Temps-temps debut pas de temps) correspondant au minimum dist2 - R2 (pas forcement sur l'intervalle)
								ex = num / denom;
								// Recherche de l'écart (Temps-temps debut pas de temps) correspondant au minimum dist2-R2 sur le pas de temps
								if(ex <= 0)
								{
									em = deltat;
									x1=Sx[i]-(xb1+(uu*em));
									x2=Sy[i]-(yb1+(vv*em));
									dist=sqrt(x1*x1+x2*x2);
									rayonf=rayon1+em*rapp1;
									ecart = dist - rayonf;
								}
								else if(ex <= deltat)
								{
									em = deltat;
									x1=Sx[i]-(xb1+(uu*em));
									x2=Sy[i]-(yb1+(vv*em));
									dist=sqrt(x1*x1+x2*x2);
									rayonf=rayon1+em*rapp1;
									ecart1 = dist - rayonf;

									em = 0;
									x1=Sx[i]-(xb1+(uu*em));
									x2=Sy[i]-(yb1+(vv*em));
									dist=sqrt(x1*x1+x2*x2);
									rayonf=rayon1+em*rapp1;
									ecart2 = dist - rayonf;

									ecart=min(ecart1,ecart2);
								}
								else
								{
									em=0;
									x1=Sx[i]-(xb1+(uu*em));
									x2=Sy[i]-(yb1+(vv*em));
									dist=sqrt(x1*x1+x2*x2);
									rayonf=rayon1+em*rapp1;
									ecart = dist - rayonf;
								}

								if(ecart <= 0)
									flagSource[i] = 1;
							}
						} 
						else
						{
							//Il est important de parcourir les points en commençant par la fin car des élements vont etre retirer aux vecteurs lors du parcours de la boucle
							for(i=Pts.size()-1;i>=0;i--)
							{
								num = 2 * uu * (PtsX[i] - xb1) + 2 * vv * (PtsY[i] - yb1) + 2 * rayon1 * rapp1;

								// Recherche de l'écart (Temps-temps debut pas de temps) correspondant au minimum dist2-R2 sur le pas de temps
								if(num >= 0)
									em = 0;
								else 
									em = deltat;

								x1=PtsX[i]-(xb1+(uu*em));
								x2=PtsY[i]-(yb1+(vv*em));
								dist=sqrt(x1*x1+x2*x2);
								rayonf=rayon1+em*rapp1;

								if(dist <= rayonf)
								{
									PtsOK.push_back(Pts[i]);
									Pts.erase(Pts.begin()+i);
									PtsX.erase(PtsX.begin()+i);
									PtsY.erase(PtsY.begin()+i);
								}
							}

							// recherche des sources se trouvant dans la bouffée
							for(i=0;i<Sx.size();i++)
							{
								num = 2 * uu * (Sx[i] - xb1) + 2 * vv * (Sy[i] - yb1) + 2 * rayon1 * rapp1;

								// Recherche de l'écart (Temps-temps debut pas de temps) correspondant au minimum dist2-R2 sur le pas de temps
								if(num >= 0)
									em = 0;
								else
									em = deltat;

								x1=Sx[i]-(xb1+(uu*em));
								x2=Sy[i]-(yb1+(vv*em));
								dist=sqrt(x1*x1+x2*x2);
								rayonf=rayon1+em*rapp1;

								if(dist <= rayonf)
									flagSource[i] = 1;
							}
						}
					}
				}

				// Modification des matrices de comptage pour la grille
				for (i=0;i<PtsOK.size();i++)
				{
					iy=PtsOK[i]/nx;
					ix=PtsOK[i]-(iy*nx);
					retro_grille[ix][iy]=retro_grille[ix][iy]+1;
					retro_grille_all[ix][iy]=retro_grille_all[ix][iy]+1;
				}

				// Modification des matrices de comptage pour les sources
				for(i=0;i<Sx.size();i++)
					if(flagSource[i]==1)
						retro_source[i][ip]=retro_source[i][ip]+1;
			}  // Boucle sur les "retrodiffusions" (Il peut y en avoir plusieurs pour 1 plainte

			// Ecriture du fichier grille pour la plainte consideree
			/*temp_sortie->Clear();
			temp_sortie->Add("X;Y;nb retro");
			for(i=0;i<nx;i++)
			{
				for(j=0;j<ny;j++)
				{
					x1=grille_xmin+i*grille_resolution;
					x2=grille_ymin+j*grille_resolution;
					temp_sortie->Add(ToString(x1) + ";" + ToString(x2) + ";" + ToString(retro_grille[i][j]));
				}
			}
			temp_sortie->SaveToFile(RepSortie + "resultat-grille-" + Pid[ip] + ".txt");*/
		} // test de detection d'erreur if(!erreur2)
			
	} // Boucle sur les plaintes


	// Ecriture du fichier grille global
	/*Log(NomFichierLog, "Ecriture du fichier resultat-grille-ToutesPlaintes...");
	temp_sortie->Clear();
	temp_sortie->Add("X;Y;nb retro");
	for(i=0;i<nx;i++)
	{
		for(j=0;j<ny;j++)
		{
			x1=grille_xmin+i*grille_resolution;
			x2=grille_ymin+j*grille_resolution;
			temp_sortie->Add(ToString(x1) + ";" + ToString(x2) + ";" + ToString((int)retro_grille_all[i][j]));
		}
	}
	temp_sortie->SaveToFile(RepSortie + "resultat-grille-ToutesPlaintes.txt");

	// Ecriture du fichier des sources
	Log(NomFichierLog, "Ecriture du fichier resultat-Source...");
	temp_sortie->Clear();
	txt = "ID Plainte;Position X (m);Position Y (m);Date de la plainte (TU);Duree de la plainte (H);Ecart max aux donnees meteo (s)";
	for(i=0;i<Sid.size();i++)
			txt = txt + ";" + Sid[i];
	temp_sortie->Add(txt);
	for(i=0;i<Pid.size();i++)
	{
		txt = Pid[i] + ";" + ToString(Px[i]) + ";" + ToString(Py[i]);
		date1 = IncSecond(date0 , Ptemps[i]);
		txt = txt + ";" + date1.FormatString("dd/mm/yyyy hh:nn") + ";" + ToString(Pdur[i]) + ";" + ToString(MaxDurInterp[i]);
		for(j=0;j<Sid.size();j++)
		{
			if(flagProbMTOPlainte[i]==0)
				txt = txt + ";" + ToString(retro_source[j][i]);
			else
				txt = txt + ";-999";
		}
		temp_sortie->Add(txt);
	}
	temp_sortie->SaveToFile(RepSortie + "resultat-Source.txt");*/
}

// Ecriture en PNG
if(retro_grille_all!=NULL)
{
	vector<Level> echelle_couleurs;
	Level level;
	level.Value = 0;       level.R = 255;  level.G = 255;   level.B = 255;	 level.A = 0;  	echelle_couleurs.push_back(level);
	level.Value = 1;       level.R = 0;    level.G = 173;   level.B = 221;   level.A = 255;	echelle_couleurs.push_back(level);
	level.Value = 2;       level.R = 95;   level.G = 192;   level.B = 140;   level.A = 255;	echelle_couleurs.push_back(level);
	level.Value = 3;       level.R = 235;  level.G = 232;   level.B = 64;    level.A = 255;	echelle_couleurs.push_back(level);
	level.Value = 4;       level.R = 250;  level.G = 167;   level.B = 58;    level.A = 255;	echelle_couleurs.push_back(level);
	level.Value = 5;       level.R = 238;  level.G = 54;    level.B = 36;    level.A = 255;	echelle_couleurs.push_back(level);
	level.Value = 6;       level.R = 182;  level.G = 31;    level.B = 57;    level.A = 255;	echelle_couleurs.push_back(level);
	level.Value = 7;       level.R = 130;  level.G = 26;    level.B = 85;    level.A = 255;	echelle_couleurs.push_back(level);
	level.Value = 100;     level.R = 130;  level.G = 26;    level.B = 85;    level.A = 255;	echelle_couleurs.push_back(level);
		
	Log(NomFichierLog, "Génération du fichier "+RepSortie+IdRequete+".png...");
	if(!writePNG(RepSortie+IdRequete+".png", nx, ny, retro_grille_all, echelle_couleurs))
		Log(NomFichierLog, "ERREUR : impossible de générer "+RepSortie+IdRequete+".png.");
}

// Desallocations
delete temp_sortie;

if(retro_grille!=NULL)
{
	for(i=0;i<nx;i++)
		delete[] retro_grille[i];
	delete[] retro_grille;
}
if(retro_grille_all!=NULL)
{
    for(i=0;i<nx;i++)
        delete[] retro_grille_all[i];
    delete[] retro_grille_all;
}
if(retro_source!=NULL)
{
    for(i=0;i<Sid.size();i++)
		delete[] retro_source[i];
    delete[] retro_source;
}
if(meteo!=NULL)
{
    for(i=0;i<nmeteo;i++)
        delete[] meteo[i];
    delete[] meteo;
}
if(DurInterp!=NULL)
{
    for(i=0;i<nmeteo;i++)
        delete[] DurInterp[i];
    delete[] DurInterp;
}
if(flagSource!=NULL)
	delete[] flagSource;

if(MaxDurInterp!=NULL)
	delete[] MaxDurInterp;

if(flagProbMTOPlainte!=NULL)
	delete[] flagProbMTOPlainte;

Log(NomFichierLog, "CALCUL TERMINE.");	
}

//---------------------------------------------------------------------------
// Fonction utilisee pour trouver la valeur de t racine de l'equation sigy(t)-sigy0=0
// Formule valable pour les cas convectifs
// Formulation utilisee = celle d'ADMS
//---------------------------------------------------------------------------
double resol_sigyCBL_ADMS(double t)
{
    double s1,s2,s;
	s1=RRRsigvc*t/sqrt(1+0.91*(t*RRRwstar)/RRRhcl);
	s2=RRRsigvn*t/sqrt(1+2.5*(t*RRRustar)/RRRhcl);
	s=sqrt(s1*s1+s2*s2)-RRRsigy;
	return(s);
}

//---------------------------------------------------------------------------
// Fonction utilisee pour trouver la valeur de t racine de l'equation sigy(t)-sigy0=0
// Formule valable pour les cas stables
// Formulation utilisee = celle d'ADMS
//---------------------------------------------------------------------------
double resol_sigySTAB_ADMS(double t)
{
    double s;
	s=(RRRsigv*t/sqrt(1+2.5*(t*RRRustar)/RRRhcl))-RRRsigy;
	return(s);
}

//---------------------------------------------------------------------------
#define ITMAX 100
#define EPS 3.0e-8
// Recherche d'une racine d'une fonction entre x1 et x2 avec la tolerance tol
// erreur :
//   0 : Tout s'est bien passe
//   1 : La racine de la fonction ne semble pas etre entre x1 et x2
//   2 : Nbre d'iteration maximum atteint
//---------------------------------------------------------------------------
double zbrent(double (*func)(double), double x1, double x2, double tol, int *erreur)
{
	int iter;
	double a=x1,b=x2,c=x2,d,e,min1,min2;
	double fa=(*func)(a),fb=(*func)(b),fc,p,q,r,s,tol1,xm;

    *erreur=0;

        //if (true) // TESTTTTTTTTTTTTTTTTTTTTTT
	if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0))
	{
		//ErreurType = "Root must be bracketed in zbrent";
		//ErreurCalcul = true;
		*erreur=1;
		return -999.0;
	}
	fc=fb;
	for (iter=1;iter<=ITMAX;iter++)
	{
		if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0))
		{
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (fabs(fc) < fabs(fb))
		{
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*fabs(b)+0.5*tol;
		xm=0.5*(c-b);
		if (fabs(xm) <= tol1 || fb == 0.0)
			return b;
		if (fabs(e) >= tol1 && fabs(fa) > fabs(fb))
		{
			s=fb/fa;
			if (a == c)
			{
				p=2.0*xm*s;
				q=1.0-s;
			}
			else
			{
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0)
				q = -q;
			p=fabs(p);
			min1=3.0*xm*q-fabs(tol1*q);
			min2=fabs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2))
			{
				e=d;
				d=p/q;
			}
			else
			{
				d=xm;
				e=d;
			}
		}
		else
		{
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (fabs(d) > tol1)
			b += d;
		else
			b += SIGN(tol1,xm);
		fb=(*func)(b);
	}
    //ErreurType = "Maximum number of iterations exceeded in zbrent";
    *erreur=2;
    return -999;
}
#undef ITMAX
#undef EPS

//---------------------------------------------------------------------------
void Log(String nom_fichier_log, String message)
{
	TStringList *log = new TStringList();
	if(FileExists(nom_fichier_log))
		log->LoadFromFile(nom_fichier_log);
	log->Add(Now().FormatString("dd/mm/yyyy hh:nn")+" "+message);
	log->SaveToFile(nom_fichier_log);
	delete log;
	cout << Now().FormatString("dd/mm/yyyy hh:nn:ss").c_str() << " " << message.c_str() << endl;
}

//---------------------------------------------------------------------------
void Log(String nom_fichier_log, TStringList *lignes)
{
	TStringList *log = new TStringList();
	if(FileExists(nom_fichier_log))
		log->LoadFromFile(nom_fichier_log);
	for(int l=0;l<lignes->Count;++l)
		log->Add(lignes->Strings[l]);
	log->SaveToFile(nom_fichier_log);
	delete log;
}

// Ajout 24 mars 2020
// Calcul de la LMO et de ustar a partir des variables suivantes (Methode type ADMS) :
// - RAY : Rayonnement en W/m2. 
// - NEB : Nebulosite en octas. 
// - ZVIT : Hauteur de mesure de la vitesse du vent en m
// - VIT : Vitesse du vent mesurée a la hauteur ZVIT en m/s
//         Les vitesses < 0.01 m/s sont ramenées a 0.01 m/s 
// - TEMP : Temperature en deg Celcius
// -  Z0 : Hauteur de rugosité (m)
// - ALB : Albedo
// - PPT : Parametre de Priestley Taylor
// 
// Sorties : 
// - LMO : LMO en m
// - USTAR : u* en m/s
//
int calLmoUstar(double z0,double alb,double ppt,double vit,double zvit, double temp,double ray,double neb,double *lmo,double *ustar)
{
   double UstarOld,LmoOld,tstar,hf,xx,xs,tempk;
   int iok,iok2,niter;
   double vk=0.4; // Constante von karman
   double gg=9.807; // acceleration pesanteur (m/s2)
   double cp=1012; // Chaleur specifique massique a pression constante (J/kg/°C)
   double dens=1.225;  // Densite de l'air (kg/m3)

   if(vit < 0.01) vit=0.01;

   tempk=temp+273.15;

   if(ray <= 20) // nuit
   {
       // Calcul thetastar
       tstar=0.09*(1-(0.5*(neb/8)*(neb/8)));

       // Calcul USTAR et LMO
       iok2=0;
       while(iok2==0)
       {
           LmoOld=1e30;
           UstarOld=calUstar(zvit, vit, LmoOld,3000,z0);
           iok=0;
           niter=0;
           while(iok==0)
           {
                *lmo=((1*tempk*UstarOld*UstarOld))/(vk*gg*tstar);
                // Dans les cas non convergent, lmo peut prendre des valeurs qui font planter le calcul de ustar 
                //  Note : LMO ne peut etre que positif car initialement theta* est positif et la facon dont il est modifie garantit qu'il le reste
                if(*lmo < 0.1) 
                {
                   iok=1;
                } else
                {
                   *ustar=calUstar(zvit, vit, *lmo,3000,z0);
                   if(fabs((1 / *lmo)-(1/LmoOld)) < 1e-10)
                   {
                       iok=1;
                       iok2=1;
                   }
                }

                LmoOld=*lmo;
                UstarOld=*ustar;

                niter=niter+1;
                if(niter > 100) iok=1;
          }

          // On verifie qu'une solution a bien ete trouvee
          // En effet, pour certaines valeurs de flux de chaleur sensible trop grande (en valeur absolue), aucune solution ne peut etre trouve
          // Dans ce cas on modifie le flux de chaleur sensible (par l'intermédiaire de theta* jusqu'à obtenir une solution
          if(iok2==0)
          {
                if(tstar > 0.001)
                {
                     tstar=tstar-0.001;
                } else tstar=tstar/2;  // Le cas neutre (t*=0) admet converge forcement
				
          }
       }

   } else  // Cas du jour
   {
       // Calcul du flux de chaleur sensible
       xx=((1-alb)*ray)+(5.31e-13*pow(tempk,6))-(5.67e-8*pow(tempk,4))+(60*(neb/(double)8));
       xx=0.9*(xx/1.12);
       xs=exp(0.055*(tempk-279));
       xs=((1-ppt)*xs+1)/(xs+1);
       hf=(xx*xs)-(20*ppt);
		
       // Calcul USTAR et LMO
       iok2=0;
       while(iok2==0)
       {
          LmoOld=1e30;
          UstarOld=calUstar(zvit, vit, LmoOld,3000,z0);
          iok=0;
          niter=0;
          while(iok==0)
          {
             *lmo=((-1*tempk*UstarOld*UstarOld*UstarOld))/(vk*gg*(hf/(dens*cp)));
             *ustar=calUstar(zvit, vit, *lmo,3000,z0);
             if(fabs((1 / *lmo)-(1/LmoOld)) < 1e-6)
             {
                iok=1;
                iok2=1;
             }

             niter=niter+1;
             if(niter > 100) iok=1;
             LmoOld=*lmo;
             UstarOld=*ustar;
          }

          // On verifie qu'une solution a bien ete trouvee
          // En effet, pour certaines valeurs de flux de chaleur sensible trop extreme, aucune solution ne peut etre trouve
          // Dans ce cas on modifie le flux de chaleur sensible jusqu'à obtenir une solution
          if(iok2==0)
          {
             if(hf < 0)
             {
                hf=hf+1;
             } else hf=hf-1;
          }
       }
   }


   // On verifie que LMO n'est pas trop petit
   // Dans le cas contraire, on pose LMO=LMOmin et on recalcule les autres termes
   if((1 / *lmo) > 1)
   {
      *lmo=1;
      *ustar=calUstar(zvit, vit, *lmo,3000,z0);
      //hf=(-1.* dens * cp * tempk * *ustar * *ustar * *ustar)/(*lmo * vk * gg);
   }

   return 1;
}

// Valeur du coefficient "a" dans la formulation des sigma-Y de Mesopuff-2. 
// Cf doc CALMET 5 p 2-34
// pgt : Classe de pasquill sous la forme d'un entier entre 1 et 6.
//       1=Classe A, 2=Classe B, ...
//
// Retourne -999 si la pgt n'est pas valide
double coeffA(int pgt)
{
    double a;

    if(pgt==1) a=0.36;
    else if(pgt==2) a=0.25;
    else if(pgt==3) a=0.19;
    else if(pgt==4) a=0.13;
    else if(pgt==5) a=0.096;
    else if(pgt==6) a=0.063;
    else a=-999;
    return a;
}