#include <stdio.h>
#include <math.h>
#include <vector>

using namespace std; 

#ifndef _Meteo_
#define _Meteo_
int    PGT2LMOH(int pgt,double *lmo,double *hcl);
double sigmaV_instable(double z,double ustar,double lmo,double hcl,double lmomin,int flagRelief);
double sigmaV_neutre(double z,double ustar,double hcl,double lmomin,int flagRelief);
double sigmaV_stable(double z,double ustar,double hcl,double z0,double lmomin,int flagRelief);
double calUstar(double z2, double vv, double lmo,double h,double z0);
double Z2Vit(double lmo,double h,double z2,double z0,double ustar);
double cal_sigySTAB_ADMS(double t,double sigv,double hcl,double ustar);
double cal_sigyCBL_ADMS(double t,double sigv,double hcl,double ustar,double lmo);
float  AngleSol(float lat,float lon,int hTU, int ijul);
int    ClassePasquill(float vitesse, float rayonnement, float nebulosite);
float  VectoriserU(float direction, float vitesse);
float  VectoriserV(float direction, float vitesse);
float  DirectionUV(float u, float v);
float  VitesseUV(float u, float v);

//---------------------------------------------------------------------------
// Calcul de la vitesse du vent en fonction de z (altitude au dessus du sol)
//---------------------------------------------------------------------------
double Z2Vit(double lmo,double h,double z2,double z0,double ustar)
{
    double z,x,x0,psi,vv;
	z=z2-z0;
	if(z <= h)
	{
		if(lmo < 0)
		{
			x=16.*(z+z0)/lmo;
			x=std::pow(1-x,(double)1./4.);
			x0=16.*z0/lmo;
			x0=std::pow(1-x0,(double)1./4.);
			psi=(std::pow(1+x,2)*(1+(x*x)))/(std::pow(1+x0,2)*(1+(x0*x0)));
			psi=log(psi)-2*(atan(x)-atan(x0));
		} else
		{
			psi=0.7*(z0/lmo)+0.75*((z0/lmo)-(5/0.35))*exp(-0.35*z0/lmo)-0.7*((z+z0)/lmo);
			psi=psi-0.75*(((z+z0)/lmo)-(5/0.35))*exp(-0.35*(z+z0)/lmo);
		}
		vv=(ustar/0.4)*(log((z+z0)/z0)-psi);
	} else
	{
		if(lmo < 0)
		{
			x=16.*(h+z0)/lmo;
			x=std::pow(1-x,(double)1./4.);
			x0=16.*z0/lmo;
			x0=std::pow(1-x0,(double)1./4.);
			psi=(std::pow(1+x,2)*(1+(x*x)))/(std::pow(1+x0,2)*(1+(x0*x0)));
			psi=log(psi)-2.*(atan(x)-atan(x0));
		} else
		{
			psi=0.7*(z0/lmo)+0.75*((z0/lmo)-(5./0.35))*exp(-0.35*z0/lmo)-0.7*((h+z0)/lmo);
			psi=psi-0.75*(((h+z0)/lmo)-(5./0.35))*exp(-0.35*(h+z0)/lmo);
		}
		vv=(ustar/0.4)*(log((h+z0)/z0)-psi);
	}
	return(vv);
}




//---------------------------------------------------------------------------
// Fonction retournant la valeur de sigy
// Formule valable pour les cas convectifs
// Formulation utilisee = celle d'ADMS
//---------------------------------------------------------------------------
double cal_sigyCBL_ADMS(double t,double sigv,double hcl,double ustar,double lmo)
{
    double s1,s2,s,wstar,sigvc,sigvn;
	wstar=(hcl*ustar*ustar*ustar)/(0.4*fabs(lmo));
    wstar=std::pow(wstar,(double)1/(double)(3));
	sigvc=sqrt(0.3*wstar*wstar);
	sigvn=sqrt((sigv*sigv)-(sigvc*sigvc));
	s1=sigvc*t/sqrt(1+0.91*(t*wstar)/hcl);
	s2=sigvn*t/sqrt(1+2.5*(t*ustar)/hcl);
	s=sqrt(s1*s1+s2*s2);
	return(s);
}

//---------------------------------------------------------------------------
// Fonction retournant la valeur de sigy
// Formule valable pour les cas stables
// Formulation utilisee = celle d'ADMS
//---------------------------------------------------------------------------
double cal_sigySTAB_ADMS(double t,double sigv,double hcl,double ustar)
{
    double s;
	s=sigv*t/sqrt(1+2.5*(t*ustar)/hcl);
	return(s);
}

//---------------------------------------------------------------------------
// Calcul de Ustar à partir d'une altitude au dessus du sol, d'une vitesse de vent à cette altitude, d'une valeur de LMO, d'une hauteur de CL et de la rugosité
//---------------------------------------------------------------------------
double calUstar(double z2, double vv, double lmo,double h,double z0)
{
    double z,x,x0,psi,ustar;
	z=z2-z0;
	if(z <= h)
	{
		if(lmo < 0)
		{
			x=16.*(z+z0)/lmo;
			x=std::pow(1-x,(double)1./4.);
			x0=16.*z0/lmo;
			x0=std::pow(1-x0,(double)1./4.);
			psi=(std::pow(1+x,2)*(1+(x*x)))/(std::pow(1+x0,2)*(1+(x0*x0)));
			psi=log(psi)-2*(atan(x)-atan(x0));
		} else
		{
			psi=0.7*(z0/lmo)+0.75*((z0/lmo)-(5/0.35))*exp(-0.35*z0/lmo)-0.7*((z+z0)/lmo);
			psi=psi-0.75*(((z+z0)/lmo)-(5/0.35))*exp(-0.35*(z+z0)/lmo);
		}
		ustar=(0.4*vv)/(log((z+z0)/z0)-psi);
	} else
	{
		if(lmo < 0)
		{
			x=16.*(h+z0)/lmo;
			x=std::pow(1-x,(double)1./4.);
			x0=16.*z0/lmo;
			x0=std::pow(1-x0,(double)1./4.);
			psi=(std::pow(1+x,2)*(1+(x*x)))/(std::pow(1+x0,2)*(1+(x0*x0)));
			psi=log(psi)-2.*(atan(x)-atan(x0));
		} else
		{
			psi=0.7*(z0/lmo)+0.75*((z0/lmo)-(5./0.35))*exp(-0.35*z0/lmo)-0.7*((h+z0)/lmo);
			psi=psi-0.75*(((h+z0)/lmo)-(5./0.35))*exp(-0.35*(h+z0)/lmo);
		}
		ustar=(0.4*vv)/(log((h+z0)/z0)-psi);
	}
	return(ustar);
}


//---------------------------------------------------------------------------
// Calcul sigV dans le cas instable
// Formulation ADMS
// flagRelief : mettre a 1 si l'option terrain complexe est utilisee dans ADMS ou si on est en présence d'un environnement complexe (bati ou relief) : Fixe un minimum au sigma
//---------------------------------------------------------------------------
double sigmaV_instable(double z,double ustar,double lmo,double hcl,double lmomin,int flagRelief)
{
    double twn,sigv2,tt,rapp,wstar;
    wstar=(hcl*ustar*ustar*ustar)/(0.4*fabs(lmo));
    wstar=std::pow(wstar,(double)1/(double)(3));
    //twc=2.1*std::pow(z/hcl,1.0/3.0)*(1-(0.8*z/hcl));
    rapp=z/hcl;
    if(rapp > 1.1)   // Dans ADMS c'est 1.2 mais 1.2 peut conduire a des problemes (exposant réel de valeur négative) dans le cas stable
    {
        rapp=1.1;
	}
	twn=1.0-(0.8*rapp);
	sigv2=0.3*wstar*wstar+4.0*twn*twn*ustar*ustar;
	tt=sqrt(sigv2);

	if(lmomin > 10 && lmomin < 30)
	{
		tt=std::max(tt,0.01*(lmomin-10));
	} else
	{
		if(lmomin >= 30)
		{
			tt=std::max(tt,(double)0.2);
		}
	}
	
	if(flagRelief==1)
	{
		tt=std::max(tt,(double)0.1);
	}

	return(tt);
}

//---------------------------------------------------------------------------
// Calcul sigV dans le cas neutre
// Formulation ADMS
// flagRelief : mettre a 1 si l'option terrain complexe est utilisee dans ADMS ou si on est en présence d'un environnement complexe (bati ou relief) : Fixe un minimum au sigma
//---------------------------------------------------------------------------
double sigmaV_neutre(double z,double ustar,double hcl,double lmomin,int flagRelief)
{
	double twn,tt,rapp;
    rapp=z/hcl;
    if(rapp > 1.1)   // Dans ADMS c'est 1.2 mais 1.2 peut conduire a des problemes (exposant réel de valeur négative) dans le cas stable
    {
		rapp=1.1;
    }
	twn=1.-(0.8*rapp);
	tt=2.*ustar*twn;

	if(lmomin > 10 && lmomin < 30)
	{
		tt=std::max(tt,0.01*(lmomin-10));
	} else
	{
		if(lmomin >= 30)
		{
			tt=std::max(tt,(double)0.2);
		}
	}

	if(flagRelief==1)
	{
		tt=std::max(tt,(double)0.1);
	}

	return(tt);
}

//---------------------------------------------------------------------------
// Calcul sigV dans le cas stable 
// Formulation ADMS
// flagRelief : mettre a 1 si l'option terrain complexe est utilisee dans ADMS ou si on est en présence d'un environnement complexe (bati ou relief) : Fixe un minimum au sigma et modifie la valeur d'un parametre
//---------------------------------------------------------------------------
double sigmaV_stable(double z,double ustar,double hcl,double z0,double lmomin,int flagRelief)
{
	double alphaS,tt,rapp;
	rapp=z/hcl;
	if(rapp > 1.1)   // Dans ADMS c'est 1.2 mais 1.2 peut conduire a des problemes (exposant réel de valeur négative) dans le cas stable
	{
		rapp=1.1;
	}
	if(flagRelief==1)
	{
		alphaS=0.5;
	} else
	{
		if(z0 <= 0.01)
		{
			alphaS=0.9;
		} else
		{
			if(z0 < 0.1)
			{
				alphaS=0.9-0.4*((z0-0.01)/0.09);
			} else
			{
				alphaS=0.5;
			}
		}
	}
	tt=1-(alphaS*rapp);
	tt=2*ustar*std::pow(tt,(double)3./4.);

	if(lmomin > 10 && lmomin < 30)
	{
		tt=std::max(tt,0.01*(lmomin-10));
	} else
	{
		if(lmomin >= 30)
		{
			tt=std::max(tt,(double)0.2);
		}
	}

	if(flagRelief==1)
	{
		tt=std::max(tt,(double)0.1);
	}

	return(tt);
}


//---------------------------------------------------------------------------
// Retourne LMO et H pour une classe de Pasquill
// Issue de la doc ADMS
// pgt : 1=Classe A, 6=Classe F
int PGT2LMOH(int pgt,double *lmo,double *hcl)
{
        if(pgt==1)
        {
                *lmo=-2;
                *hcl=1300;
        } else
        {
                if (pgt==2)
                {
                        *lmo=-10;
                        *hcl=900;
                } else
                {
                        if (pgt==3)
                        {
                                *lmo=-100;
                                *hcl=850;
                        } else
                        {
                               if (pgt==4)
                               {
                                        *lmo=1e10;
                                        *hcl=800;
                               } else
                               {
                                        if (pgt==5)
                                        {
                                                *lmo=100;
                                                *hcl=400;
                                        } else
                                        {
                                                if (pgt==6)
                                                {
                                                        *lmo=20;
                                                        *hcl=100;
                                                } else
                                                {
                                                        return 1; // Classe PGT invalide
                                                }
                                        }
                               }
                        }
                }
        }
        return 0;
}

//---------------------------------------------------------------------------
// Retourne l'angle d'élévation solaire
//---------------------------------------------------------------------------
float AngleSol(float lat,float lon,int hTU, int ijul)
{
        float ds,latrad,lonrad,sinP;
        float tt1,tt2,tt3;
        ds=0.409*cos((2.*3.141592*(ijul-173.))/365.25);

        //lat2=StrToFloat(EditHeureDeb1->Text)

        latrad=lat/180.*3.141592;
        lonrad=(lon/180.)*3.141592;

        tt1=(3.141592*(hTU-12.))/12.;
        tt2=tt1+lonrad;
        tt3=cos(tt2);

        sinP=sin(latrad)*sin(ds)+cos(latrad)*cos(ds)*tt3;
        return asin (sinP);
}

//---------------------------------------------------------------------------
// Calcule la classe de Pasquill (PGT) en fonction des paramètres fournis :
// vitesse en m/s, rayonnement solaire en W/m² et nébulosité en octats
//---------------------------------------------------------------------------
int ClassePasquill(float vitesse, float rayonnement, float nebulosite)
{
int classe;

// Par défaut, on retourne "neutre", c'est à dire la classe D ou 4
classe = 4;

if(rayonnement>=20)      // Pour le jour
{
        if((vitesse>=0)&&(rayonnement>=0))       // Si les valeurs sont valides
        {
                if(rayonnement>=925)
                {
                        if(vitesse<2)
                                classe = 1;                                             // A
                        else
                        {
                                if(vitesse<3)
                                        classe = 1;                                     // A
                                else
                                {
                                        if(vitesse<5)
                                                classe = 2;                             // B
                                        else
                                        {
                                                if(vitesse<6)
                                                        classe = 3;                     // C
                                                else
                                                        classe = 3;                     // C
                                        }
                                }
                        }
                }
                else
                {
                        if(rayonnement>=675)
                        {
                                if(vitesse<2)
                                        classe = 1;                                     // A
                                else
                                {
                                        if(vitesse<3)
                                                classe = 2;                             // B
                                        else
                                        {
                                                if(vitesse<5)
                                                        classe = 2;                     // B
                                                else
                                                {
                                                        if(vitesse<6)
                                                                classe = 3;             // C
                                                        else
                                                                classe = 4;             // D
                                                }
                                        }
                                }
                        }
                        else
                        {
                                if(rayonnement>=175)
                                {
                                        if(vitesse<2)
                                                classe = 2;                             // B
                                        else
                                        {
                                                if(vitesse<3)
                                                        classe = 3;                     // C
                                                else
                                                {
                                                        if(vitesse<5)
                                                                classe = 3;             // C
                                                        else
                                                        {
                                                                if(vitesse<6)
                                                                        classe = 4;     // D
                                                                else
                                                                        classe = 4;     // D
                                                        }
                                                }
                                        }
                                }
                                else
                                {
                                        if(vitesse<2)
                                                classe = 4;                             // D
                                        else
                                        {
                                                if(vitesse<3)
                                                        classe = 4;                     // D
                                                else
                                                {
                                                        if(vitesse<5)
                                                                classe = 4;             // D
                                                        else
                                                        {
                                                                if(vitesse<6)
                                                                        classe = 4;     // D
                                                                else
                                                                        classe = 4;     // D
                                                        }
                                                }
                                        }
                                }
                        }
                }
        }
}
else                     // Pour la nuit
{
        if((vitesse>=0)&&(nebulosite>=0))       // Si les valeurs sont valides
        {
                if(nebulosite<=3)
                {
                        if(vitesse<0.8)
                                classe = 6;                                                                     // F
                        else
                        {
                                if(vitesse<1.9)
                                        classe = 6;                                                             // F
                                else
                                {
                                        if(vitesse<2.9)
                                                classe = 6;                                                     // F
                                        else
                                        {
                                                if(vitesse<3.4)
                                                        classe = 6;                                             // F
                                                else
                                                {
                                                        if(vitesse<3.9)
                                                                classe = 5;                                     // E
                                                        else
                                                        {
                                                                if(vitesse<4.9)
                                                                        classe = 5;                             // E
                                                                else
                                                                {
                                                                        if(vitesse<5.5)
                                                                                classe = 5;                     // E
                                                                        else
                                                                        {
                                                                                if(vitesse<6)
                                                                                        classe = 4;             // D
                                                                                else
                                                                                        classe = 4;             // D
                                                                        }
                                                                }
                                                        }
                                                }
                                        }
                                }
                        }
                }
                else
                {
                        if(nebulosite<=7)
                        {
                                if(vitesse<0.8)
                                        classe = 6;                                                             // F
                                else
                                {
                                        if(vitesse<1.9)
                                                classe = 6;                                                     // F
                                        else
                                        {
                                                if(vitesse<2.9)
                                                        classe = 5;                                             // E
                                                else
                                                {
                                                        if(vitesse<3.4)
                                                                classe = 5;                                     // E
                                                        else
                                                        {
                                                                if(vitesse<3.9)
                                                                        classe = 4;                             // D
                                                                else
                                                                {
                                                                        if(vitesse<4.9)
                                                                                classe = 4;                     // D
                                                                        else
                                                                        {
                                                                                if(vitesse<5.5)
                                                                                        classe = 4;             // D
                                                                                else
                                                                                {
                                                                                        if(vitesse<6)
                                                                                                classe = 4;     // D
                                                                                        else
                                                                                                classe = 4;     // D
                                                                                }
                                                                        }
                                                                }
                                                        }
                                                }
                                        }
                                }
                        }
                        else
                        {
                                if(vitesse<0.8)
                                        classe = 4;                                                             // D
                                else
                                {
                                        if(vitesse<1.9)
                                                classe = 4;                                                     // D
                                        else
                                        {
                                                if(vitesse<2.9)
                                                        classe = 4;                                             // D
                                                else
                                                {
                                                        if(vitesse<3.4)
                                                                classe = 4;                                     // D
                                                        else
                                                        {
                                                                if(vitesse<3.9)
                                                                        classe = 4;                             // D
                                                                else
                                                                {
                                                                        if(vitesse<4.9)
                                                                                classe = 4;                     // D
                                                                        else
                                                                        {
                                                                                if(vitesse<5.5)
                                                                                        classe = 4;             // D
                                                                                else
                                                                                {
                                                                                        if(vitesse<6)
                                                                                                classe = 4;     // D
                                                                                        else
                                                                                                classe = 4;     // D
                                                                                }
                                                                        }
                                                                }
                                                        }
                                                }
                                        }
                                }
                        }
                }
        }
}

return classe;
}

//---------------------------------------------------------------------------
float VectoriserU(float direction, float vitesse)
{
float u;
float deg = 3.1415926535/180.0;

if((int(direction)==-999)||(int(vitesse)==-999)||(int(direction)==-9999)||(int(vitesse)==-9999)||(isnan(direction))||(isnan(vitesse)))
        u = -999;
else
	u=-vitesse*sin(direction*deg);

return u;
}

//---------------------------------------------------------------------------
float VectoriserV(float direction, float vitesse)
{
float v;
float deg = 3.1415926535/180.0;

if((int(direction)==-999)||(int(vitesse)==-999)||(int(direction)==-9999)||(int(vitesse)==-9999)||(isnan(direction))||(isnan(vitesse)))
        v = -999;
else
        v=-vitesse*cos(direction*deg);

return v;
}

//---------------------------------------------------------------------------
float DirectionUV(float u, float v)
{
float direction;
float deg = 3.1415926535/180.0;

if((int(u)==-999)||(int(v)==-999)||(int(u)==-9999)||(int(v)==-9999)||(isnan(u))||(isnan(v)))
        direction = -999;
else
{
        if(u==0)
        {
    	        if(v>0)
        	        direction=180;
                else
                        direction=0;
	}
        else
        {
    	        if(u>0 && v>=0)
           	        direction=270-(atan(v/u)/deg);

                if(u<0 && v>=0)
                        direction=90+atan(v/-u)/deg;

		if(u<0 && v<0)
                        direction=90-atan(-v/-u)/deg;

                if(u>0 && v<0)
                        direction=270+atan(-v/u)/deg;
	}
}

return direction;
}

//---------------------------------------------------------------------------
float VitesseUV(float u, float v)
{
float vitesse;

if((int(u)==-999)||(int(v)==-999)||(int(u)==-9999)||(int(v)==-9999)||(isnan(u))||(isnan(v)))
        vitesse = -999;
else
	vitesse = sqrt(u*u + v*v);

return vitesse;
}

#endif