1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
| double TGeomGeo::LongRayon(double La, double Ma, double Lb, double Mb )
{
// ce sont les paramètres de WGS84
double R0=6378137.00;
double f=1.0/298.257222101;
// Cette formule n'est peut-être pas très bonne ?? Si voir Notions de géodésie page 13
double s=sin((Ma+Mb)*M_PI/400); // sinus de la latitude moyenne
double e2=2.0*f - f*f;
double s2=s*s;
double R1=R0*((double)1.0 - e2)/pow(((double)1.0-e2*s2),(double)(1.5));
double R2=R0/sqrt(1.0-e2*s2);
double RR=sqrt(R1 * R2);
return RR;
}
double TGeomGeo::CalcAzimut(double La, double Ma, double Lb, double Mb )
{
// Dans le triangle PAB où P est le pôle Nord
// relation des sin
double ArcAB=CalcArc(La,Ma,Lb,Mb); //en radians
//fprintf(espion,"CalcAzimut :: ArcAB=%frad\n",ArcAB);
if (ArcAB != 0.0)
{
double Zab=asin(sin((Lb-La)*M_PI/200.0) * cos(Mb*M_PI/200.0) / sin(ArcAB))*200.0/M_PI;
if ((Lb-La) < 0.0) Zab=200.0-Zab;
if (Zab < 0.0) Zab+=400.0;
return Zab;
}
else return 0.0;
}
double TGeomGeo::CalcLongArc(double La, double Ma, double Lb, double Mb )
{
double Arc=CalcArc(La,Ma,Lb,Mb);
double RR=LongRayon(La,Ma,Lb,Mb);
return Arc * RR;
}
double TGeomGeo::CalcArc(double La, double Ma, double Lb, double Mb )
{
// les valeurs ont été modifiées suivant l'unité iUnite2
// Maintenant, on est toujours en degrés décimaux
//fprintf(espion,"\nCalcArc:: La=%f Ma=%f Lb=%f Mb=%f\n",La,Ma,Lb,Mb);
double Arc=acos(sin(Ma*M_PI/200.0)*sin(Mb*M_PI/200.0) +
cos(Ma*M_PI/200.0)*cos(Mb*M_PI/200.0)*cos((La-Lb)*M_PI/200));
// Longitide Greenwitch=0: L
// Latitude Paris = 50; M
//fprintf(espion,"AB calculé =%f Longueur de l'Arc=%f\n",AB,LArc);
return Arc;
} |
Partager