#include <iostream>
#include <math.h>
#include <cstdlib>
#include<bits/stdc++.h>
using namespace std;
typedef struct ANGLE_DEG;
struct ANGLE_DEG
{
double Deg;
double Min;
double Sec;
};
double Orthonomie(double Long1, double Lat1, double Long2, double Lat2);
void Angle_Rad_Deg(double Ang_Rad, ANGLE_DEG *Angle_Deg);
double Long0_Lat0(double Lat0, double Long_centre, double Lat_centre, double R_Rech);
double Rterre = 6367.4633; // Rayon de la terre en km
// pour 1 minute -> 1 mile marin = 1852.22 metres
double Pi = 4 * atan(1);
int main(int argc, char *argv[])
{
int Lat_Deg, Lat_Min, Lat_Sec;
int Long_Deg, Long_Min, Long_Sec;
double Lat_centre, Long_centre, Lat0, Lat0_max, Lat0_min;
double R_Rech;
double Phi, Pas_Lat0;
ANGLE_DEG Lat; ANGLE_DEG Long;
int i;
double Lat0_Long0(double Long0, double Long_centre, double Lat_centre, double R_Rech);
// Entrée Latitude
cout << "Latitude du centre de recherche en Degres Minute Seconde \n Entrer Valeur degres ? " << endl;
cin >> Lat_Deg;
cout << "\n Entrer Valeur minutes ? " << endl;
cin >> Lat_Min;
cout << "\n Entrer Valeur secondes ? " << endl;
cin >> Lat_Sec;
cout << "Latitude du centre de recherche = " << Lat_Deg << " Degres " << Lat_Min << " Minutes " << Lat_Sec << "Secondes" << endl;
Lat_centre = Lat_Deg + Lat_Min / 60.0 + Lat_Sec / 3600.0 ;
Lat_centre = Lat_centre * Pi / 180.0; // Conversion en radians
cout << "Latitude du centre de recherche en radians = " << Lat_centre << endl;
// Entrée Longitude
cout << "\n Longitude du centre de recherche en Degres Minute Seconde \n Entrer Valeur degres ? " << endl;
cin >> Long_Deg;
cout << "\n Entrer Valeur minutes ? " << endl;
cin >> Long_Min;
cout << "\n Entrer Valeur secondes ? " << endl;
cin >> Long_Sec;
cout << "Longitude du centre de recherche = " << Long_Deg << " Degres " << Long_Min << " Minutes " << Long_Sec << "Secondes" << endl;
Long_centre = Long_Deg + Long_Min / 60.0 + Long_Sec / 3600.0 ;
Long_centre = Long_centre * Pi / 180.0; // Conversion en radians
cout << "Longitude du centre de recherche en radians = " << Long_centre << endl;
cout << " Entrer la Valeur du rayon de recherche en km \n " << endl;
cin >> R_Rech;
cout << " Rayon de recherche en km \n " << R_Rech << endl;
//Long_centre = Pi/2.1;
//Lat_centre = Pi/2.09;
//R_Rech = 250.0; // en km
Phi = R_Rech / Rterre;
cout << " Phi = " << Phi <<endl;
Lat0_max = Lat_centre + Phi;
Lat0_min = Lat_centre - Phi;
cout << " Lat_centre = " << Lat_centre <<endl;
cout << " Long_centre = " << Long_centre <<endl;
cout << " Lat0_min = " << Lat0_min <<endl;
cout << " Lat0_max = " << Lat0_max <<endl;
cout << " Pi/2 = " << Pi/2.0 <<endl;
if (Lat0_min >= Pi/2.0 || Lat0_max >= Pi/2.0)
{
cout << " CALCUL IMPOSSIBLE Lat0_max >= 90°" << endl;
system ("PAUSE");
}
// Verif
cout << " ORTHO LAT0 MAX = " << Orthonomie(Long_centre, Lat_centre,Long_centre, Lat0_max) << endl;
cout << " ORTHO LAT0 MIN = " << Orthonomie(Long_centre, Lat_centre,Long_centre, Lat0_min) << endl;
Pas_Lat0 = (Lat0_max - Lat0_min) / 180.0; // en radians 1 degres
cout << " Pas_Lat0 = " << Pas_Lat0 << endl;
Phi = Pas_Lat0 * Pi / 180.0;
cout << " Phi en degres = " << Phi <<endl;
double Lg0;
for (i = 0; i <= 180; i++)
{
Lat0 = Lat0_min + Pas_Lat0 * i ;
Lg0 = Lat0_Long0(Lat0, Long_centre, Lat_centre, R_Rech);
Angle_Rad_Deg(Lat0, &Lat);
Angle_Rad_Deg(Lg0, &Long); // cor biname Lat0 to Lg0
cout << setprecision(6) << " Lat0 = " << Lat0 << " " << int(Lat.Deg) << " deg " << int(Lat.Min) << "' " <<
int(Lat.Sec) << "'' " << " Long0 = " << Lg0 << " " << int(Long.Deg) << " deg " << int(Long.Min) <<
"' " << int(Long.Sec) << " '' " << " i= " << i << endl;
}
return 0;
}
//******************************************************************************
// Orthonomie
// Calcul la longueur de la route orthonomique entre les points de
// longitude (Long1, Long2) et de latitude (Lat1, Lat2) exprimes en Radians
// Orthonomie en km
//******************************************************************************
double Orthonomie(double Long1, double Lat1, double Long2, double Lat2)
{
double Corde;
Corde = pow((sin(Long1) * cos(Lat1) - sin(Long2) * cos(Lat2)), 2);
Corde = Corde + pow((cos(Long1) * cos(Lat1) - cos(Long2) * cos(Lat2)),2);
Corde = Corde + pow((sin(Lat1) - sin(Lat2)),2);
Corde = sqrt(Corde);
return( 2.0 * Rterre * (asin(Corde / 2.0)));
}
void Angle_Rad_Deg(double Angle_Rad, ANGLE_DEG *Angle_Deg)
{
Angle_Deg-> Deg = Angle_Rad * 180.0/Pi;
Angle_Deg-> Min = (Angle_Deg -> Deg - int(Angle_Deg -> Deg)) * 60.0 ;
Angle_Deg-> Sec = (Angle_Deg-> Min - int(Angle_Deg-> Min)) * 60.0;
}
double Lat0_Long0(double Lat0, double Long_centre, double Lat_centre, double R_Rech)
//**********************************************************************************************************
// Renvoie la longitude Long0 en fonction de la latitude Lat0 et des coordonées du centre de recherches :
// (Long_centre, Lat_centre en radians) et du Rayon de recherche en km
//**********************************************************************************************************
{
int i;
double Ortho0 = 0.0;
double Epsilon = 0.000001; // précision du calcul
double Pas = R_Rech /100.0; // Valeur initiale du pas
double Long0_1;
double Long_Depart = Long_centre ;
// double Ecart0 = 100000.0; double Ecart1 = 0.0;
T10:
for ( i = 0; i < 50; i++)
{
Long0_1 = Long_Depart + Pas * i;
Ortho0 = Orthonomie(Long0_1, Lat0, Long_centre, Lat_centre);
if(abs(Ortho0 - R_Rech )<= Epsilon)
{
return ( Long0_1);
}
if ((Ortho0 - R_Rech) >= 0)
{
Long_Depart = Long0_1 - Pas;
Pas = Pas /2.0;
// system ("PAUSE");
goto T10;
}
else
{
Long_Depart = Long0_1 + Pas;
Pas = Pas /2.0;
// system ("PAUSE");
goto T10;
}
}
}
