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

Last change on this file since 36 was 31, checked in by sgosseli, 12 years ago

Huge commit: use the new includes style in all the files, add the license header in all the headers, and in some cpp.

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