source: pacpusframework/trunk/include/Pacpus/PacpusTools/geodesie.h@ 18

Last change on this file since 18 was 3, checked in by sgosseli, 12 years ago
  • Add the existing Pacpus files from pacpusdev and pacpuscore.
  • Provide a clean build system based on multiple CMake files.
File size: 5.2 KB
Line 
1#ifndef GEODESIE_H
2#define GEODESIE_H
3
4#include <cmath>
5#include <iostream>
6#include <vector>
7
8namespace Geodesie {
9
10#ifndef M_PI
11# define M_PI 3.14159265358979323846
12#endif
13#ifndef M_PI_2
14# define M_PI_2 1.57079632679489661923
15#endif
16#ifndef M_PI_4
17# define M_PI_4 0.78539816339744830962
18#endif
19
20////////////////////////////////////////////////////////////////////////
21struct Matrice {
22 Matrice(const Matrice & A);
23 Matrice();
24 void Apply(double v0, double v1, double v2, double & Mv0, double & Mv1, double & Mv2);
25 double c0_l0;double c1_l0;double c2_l0;
26 double c0_l1;double c1_l1;double c2_l1;
27 double c0_l2;double c1_l2;double c2_l2;
28}; // class
29
30Matrice TransMat(const Matrice A);
31
32Matrice ProdMat(const Matrice A,const Matrice B);
33void Write(const Matrice A,std::ostream& out);
34
35////////////////////////////////////////////////////////////////////////
36class Raf98 {
37private :
38 std::vector<double> m_dvalues;
39 double LitGrille(unsigned int c,unsigned int l) const;
40public :
41 ~Raf98();
42 Raf98() {}
43 bool Load(const std::string & s);
44 bool Interpol(double longitude/*deg*/, double latitude/*deg*/, double* Hwgs84) const;
45}; // class
46////////////////////////////////////////////////////////////////////////
47
48////////////////////////////////////////////////////////////////////////
49inline double Deg2Rad(double deg) {return deg*M_PI/180.0;}
50inline double Rad2Deg(double rad) {return rad*180.0/M_PI;}
51////////////////////////////////////////////////////////////////////////
52
53const double a_Lambert93=6378137;
54const double f_Lambert93=1 / 298.257222101;
55const double e_Lambert93=sqrt(f_Lambert93*(2-f_Lambert93));
56const double lambda0_Lambert93=Deg2Rad(3.0);//degres
57const double phi0_Lambert93=Deg2Rad(46.5);
58const double phi1_Lambert93=Deg2Rad(44.0);
59const double phi2_Lambert93=Deg2Rad(49.0);//degres
60const double X0_Lambert93=700000;//
61const double Y0_Lambert93=6600000;//
62const double n_Lambert93 = 0.7256077650;
63const double c_Lambert93 = 11754255.426;
64const double xs_Lambert93 = 700000;
65const double ys_Lambert93 = 12655612.050;
66
67const double GRS_a = 6378137;
68const double GRS_f = 1/298.257222101;
69const double GRS_b = GRS_a*(1-GRS_f);
70const double GRS_e = sqrt((pow(GRS_a,2) - pow(GRS_b,2)) / pow(GRS_a,2));
71
72////////////////////////////////////////////////////////////////////////
73void Geographique_2_Lambert93(const Raf98& raf98,double lambda,double phi,double he,Matrice in,double& E,double& N,double& h,Matrice& out);
74void Geographique_2_Lambert93(const Raf98& raf98,double lambda,double phi,double he,double& E,double& N,double& h);
75void Lambert93_2_Geographique(const Raf98& raf98,double E,double N,double h,double& lambda,double& phi,double& he);
76void Lambert93_2_Geographique(const Raf98& raf98,double E,double N,double h,Matrice in,double& lambda,double& phi,double& he,Matrice& out);
77/** Convert from geographique to ECEF.
78 * @param[in] longitude Longitude in radian.
79 * @param[in] latitude Latitude in radian.
80 * @param[in] he Height in meter.
81 */
82void Geographique_2_ECEF(double longitude, double latitude, double he, double& x, double& y, double& z);
83/** Convert from ECEF two ENU.
84 * @param[in] lon0 Longitude of the origin in radian.
85 * @param[in] lat0 Latitude of the origin in radian.
86 * @param[in] he0 Height of the origin in radian.
87 */
88void ECEF_2_ENU(double x,double y,double z,double& e,double& n,double& u,double lon0,double lat0,double he0);
89////////////////////////////////////////////////////////////////////////
90
91//ALGO0001
92double LatitueIsometrique(double latitude,double e);
93//ALGO0002
94double LatitueIsometrique2Lat(double latitude_iso,double e,double epsilon);
95
96//ALGO0003
97void Geo2ProjLambert(
98 double lambda,double phi,
99 double n, double c,double e,
100 double lambdac,double xs,double ys,
101 double& X,double& Y);
102//ALGO0004
103void Proj2GeoLambert(
104 double X,double Y,
105 double n, double c,double e,
106 double lambdac,double xs,double ys,
107 double epsilon,
108 double& lambda,double& phi);
109
110double ConvMerApp(double longitude);
111
112/**
113Converts Cartesian (x, y) coordinates to polar coordinates (r, theta)
114*/
115template <typename _T1, typename _T2>
116void cartesianToPolar(const _T1 x, const _T1 y, _T2 & r, _T2 & theta) {
117 r = std::sqrt(x*x + y*y);
118 theta = std::atan2(x, y);
119}
120
121/**
122Converts polar coordinates (r, theta) to Cartesian (x, y) coordinates
123*/
124template <typename _T1, typename _T2>
125void polarToCartesian(const _T1 r, const _T1 theta, _T2 & x, _T2 & y) {
126 x = r * std::cos(theta);
127 y = r * std::sin(theta);
128}
129
130/**
131Converts Cartesian (x, y, z) coordinates to spherical coordinates (r, theta, phi)
132Angles expressed in radians.
133*/
134template <typename _T1, typename _T2>
135void cartesianToSpherical(const _T1 x, const _T1 y, const _T1 z, _T2 & r, _T2 & theta, _T2 & phi) {
136 r = std::sqrt(x*x + y*y + z*z);
137 theta = std::acos(z / r);
138 phi = std::atan2(y, x);
139}
140
141/**
142Converts spherical coordinates (r, theta, phi) to Cartesian (x, y, z) coordinates.
143Angles expressed in radians.
144*/
145template <typename _T1, typename _T2>
146void sphericalToCartesian(const _T1 r, const _T1 theta, const _T1 phi, _T2 & x, _T2 & y, _T2 & z) {
147 x = r * std::sin(theta) * std::cos(phi);
148 y = r * std::sin(theta) * std::sin(phi);
149 z = r * std::cos(theta);
150}
151
152} // namespace Geodesie
153
154#endif // GEODESIE_H
Note: See TracBrowser for help on using the repository browser.