- Timestamp:
- Nov 6, 2013, 11:22:15 AM (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/PacpusTools/src/geodesie.cpp
r198 r208 2 2 // This file is part of the PACPUS framework distributed under the 3 3 // CECILL-C License, Version 1.0. 4 /// @author Marek Kurdej <firstname.surname@utc.fr> 4 5 /// @author Jean Laneurit <firstname.surname@utc.fr> 5 6 /// @date April, 2010 … … 10 11 #include <fstream> 11 12 12 //using boost::math::constants::pi; 13 using boost::math::constants::pi; 14 using boost::math::constants::half_pi; 13 15 14 16 #ifdef _MSC_VER … … 16 18 #endif //_MSC_VER 17 19 18 namespace Geodesie { 19 /// //////////////////////////////////////////////////////////////////// 20 Matrice::Matrice(const Matrice & A) { 20 namespace Geodesie 21 { 22 23 /// //////////////////////////////////////////////////////////////////// 24 Matrice::Matrice(const Matrice & A) 25 { 21 26 c0_l0=A.c0_l0;c1_l0=A.c1_l0;c2_l0=A.c2_l0; 22 27 c0_l1=A.c0_l1;c1_l1=A.c1_l1;c2_l1=A.c2_l1; 23 28 c0_l2=A.c0_l2;c1_l2=A.c1_l2;c2_l2=A.c2_l2; 24 29 } 25 /// //////////////////////////////////////////////////////////////////// 26 Matrice::Matrice() { 30 31 /// //////////////////////////////////////////////////////////////////// 32 Matrice::Matrice() 33 { 27 34 c0_l0=0.0;c1_l0=0.0;c2_l0=0.0; 28 35 c0_l1=0.0;c1_l1=0.0;c2_l1=0.0; 29 36 c0_l2=0.0;c1_l2=0.0;c2_l2=0.0; 30 37 } 31 /// //////////////////////////////////////////////////////////////////// 32 void Matrice::Apply(double v0,double v1,double v2, double& Mv0,double& Mv1,double& Mv2) { 38 39 /// //////////////////////////////////////////////////////////////////// 40 void Matrice::Apply(double v0,double v1,double v2, double& Mv0,double& Mv1,double& Mv2) 41 { 33 42 Mv0 = c0_l0*v0 + c1_l0*v1 + c2_l0*v2; 34 43 Mv1 = c0_l1*v0 + c1_l1*v1 + c2_l1*v2; 35 44 Mv2 = c0_l2*v0 + c1_l2*v1 + c2_l2*v2; 36 45 } 37 /// //////////////////////////////////////////////////////////////////// 38 Matrice ProdMat(const Matrice A, const Matrice B) { 46 47 /// //////////////////////////////////////////////////////////////////// 48 Matrice ProdMat(const Matrice A, const Matrice B) 49 { 39 50 Matrice out; 40 51 … … 54 65 55 66 /// //////////////////////////////////////////////////////////////////// 56 Matrice TransMat(const Matrice A) { 67 Matrice TransMat(const Matrice A) 68 { 57 69 Matrice out; 58 70 out.c0_l0=A.c0_l0 ; out.c1_l0 = A.c0_l1 ; out.c2_l0 = A.c0_l2 ; … … 112 124 return true; 113 125 } 114 /// //////////////////////////////////////////////////////////////////// 115 double Raf98::LitGrille(unsigned int c,unsigned int l) const{ 126 127 /// //////////////////////////////////////////////////////////////////// 128 double Raf98::LitGrille(unsigned int c,unsigned int l) const 129 { 116 130 const unsigned int w=421; 117 131 // const unsigned int h=381; 118 132 return m_dvalues.at(c+l*w); 119 133 } 120 /// //////////////////////////////////////////////////////////////////// 121 bool Raf98::Load(const std::string & sin) { 122 std::ifstream in(sin.c_str()); 134 135 /// //////////////////////////////////////////////////////////////////// 136 bool Raf98::Load(const std::string & s) 137 { 138 std::ifstream in(s.c_str()); 123 139 unsigned int w = 421; 124 140 unsigned int h = 381; … … 161 177 /// //////////////////////////////////////////////////////////////////// 162 178 //ALGO0001 163 double Geodesie::LatitueIsometrique(double latitude,double e) { 164 double li; 165 li = log(tan(M_PI_4 + latitude/2.)) + e*log( (1-e*sin(latitude))/(1+e*sin(latitude)) )/2; 179 double Geodesie::LatitueIsometrique(double latitude, double e) 180 { 181 using ::std::log; 182 using ::std::sin; 183 using ::std::tan; 184 185 double li = log(tan(pi<double>() / 4. + latitude/2.)) + e*log( (1-e*sin(latitude))/(1+e*sin(latitude)) )/2; 166 186 return li; 167 187 } … … 169 189 /// //////////////////////////////////////////////////////////////////// 170 190 //ALGO0002 171 double Geodesie::LatitueIsometrique2Lat(double latitude_iso,double e,double epsilon) { 172 double latitude_i=2*atan(exp(latitude_iso)) - M_PI_2; 173 double latitude_ip1=latitude_i+epsilon*2; 174 while (fabs(latitude_i-latitude_ip1)>epsilon) { 175 latitude_i=latitude_ip1; 176 latitude_ip1=2*atan( 191 double Geodesie::LatitueIsometrique2Lat(double latitude_iso,double e,double epsilon) 192 { 193 using ::std::atan; 194 using ::std::exp; 195 using ::std::abs; 196 using ::std::log; 197 using ::std::sin; 198 199 double latitude_i = 2 * atan(exp(latitude_iso)) - half_pi<double>(); 200 double latitude_ip1 = latitude_i + epsilon * 2; 201 while (abs(latitude_i-latitude_ip1) > epsilon) { 202 latitude_i = latitude_ip1; 203 latitude_ip1 = 2 * atan( 177 204 exp(e*0.5* 178 205 log( … … 181 208 ) 182 209 *exp(latitude_iso) 183 ) - M_PI_2;210 ) - half_pi<double>(); 184 211 } 185 212 return latitude_ip1; … … 222 249 223 250 //////////////////////////////////////////////////////////////////// 224 void Geodesie::Geographique_2_Lambert93(const Raf98& raf98,double lambda,double phi,double he,Matrice in,double& E,double& N,double& h,Matrice& out) { 251 void Geodesie::Geographique_2_Lambert93(const Raf98& raf98,double lambda,double phi,double he,Matrice in,double& E,double& N,double& h,Matrice& out) 252 { 225 253 Matrice passage; 226 254 double conv=Geodesie::ConvMerApp(lambda); … … 263 291 h=he-diff_h; 264 292 } 265 /** 266 Converts Lambert93 coordinates (East, North, Height) into geographical coordinates in radians (Longitude = Rad2Deg(lambda), Latitude = Rad2Deg(phi), Height) 267 */ 293 294 /// Converts Lambert93 coordinates (East, North, Height) into geographical coordinates in radians (Longitude = Rad2Deg(lambda), Latitude = Rad2Deg(phi), Height) 268 295 void Geodesie::Lambert93_2_Geographique(const Raf98& raf98,double E,double N,double h,double& lambda,double& phi,double& he) { 269 296 Geodesie::Proj2GeoLambert( … … 278 305 he=h+diff_h; 279 306 } 307 280 308 //////////////////////////////////////////////////////////////////////// 281 309 void Geodesie::Lambert93_2_Geographique(const Raf98& raf98,double E,double N,double h,Matrice in,double& lambda,double& phi,double& he,Matrice& out) { … … 348 376 } 349 377 350 QMatrix4x4 Geodesie::yprenuToMatrix(QVector3D angle, QVector3D position) {351 378 QMatrix4x4 Geodesie::yprenuToMatrix(QVector3D angle, QVector3D position) 379 { 352 380 float c1 = cos(angle.x()); 353 381 float c2 = cos(angle.y());
Note:
See TracChangeset
for help on using the changeset viewer.