3 #include <HYDROData_Lambert93.h>
11 const double PI = 3.14159265;
15 // Base values of the Lambert-93
17 const double a = 6378137; // m -- le demi-grand axe
19 const double f = 1.0 / 298.257222101; // l'aplatissement
23 const double phi_1_deg = 44; // deg -- le premier parallèle d'échelle
25 const double phi_2_deg = 49; // deg -- le deuxième parallèle d'échell
29 const double lambda_0_deg = 3; // deg -- la longitude d'origine donnée par le méridien central de Greenwich
31 const double phi_0_deg = 46.5; // deg -- la latitude d'origine
35 const double X_0 = 700000; // m -- la coordonnée à l'origine
37 const double Y_0 = 6600000; // m -- la coordonnée à l'origine
41 // Derived values of the Lambert-93
43 const double b = a * ( 1 - f ); // m -- le demi-petit axe
45 const double e = sqrt( a*a - b*b ) / a; // l'excentricité
49 const double phi_0 = HYDROData_Lambert93::toRad( phi_0_deg );
51 const double phi_1 = HYDROData_Lambert93::toRad( phi_1_deg );
53 const double phi_2 = HYDROData_Lambert93::toRad( phi_2_deg );
55 const double lambda_0 = HYDROData_Lambert93::toRad( lambda_0_deg );
61 double cot( double x )
65 return cos( x ) / sin( x );
81 const double s1 = sin( phi_1 );
83 const double s2 = sin( phi_2 );
85 const double c1 = cos( phi_1 );
87 const double c2 = cos( phi_2 );
91 const double n1 = ln( c2/c1 ) + 1.0/2.0 * ln( (1-e*e*s1*s1)/(1-e*e*s2*s2) );
93 const double n2 = tan( phi_1 / 2 + PI/4 ) * pow( 1-e*s1, e/2 ) * pow( 1+e*s2, e/2 );
95 const double n3 = tan( phi_2 / 2 + PI/4 ) * pow( 1+e*s1, e/2 ) * pow( 1-e*s2, e/2 );
97 const double n = n1 / ( ln( n2/n3 ) );
101 const double p1 = a * c1 / ( n * sqrt( 1-e*e*s1*s1 ) );
103 const double p2 = tan( phi_1 / 2 + PI / 4 );
105 const double p3 = pow( (1-e*s1)/(1+e*s1), e/2 );
107 const double p_0 = p1 * pow( p2*p3, n );
111 double HYDROData_Lambert93::toRad( double theDeg )
115 return theDeg * PI / 180.0;
121 double HYDROData_Lambert93::toDeg( double theRad )
125 return theRad / PI * 180.0;
131 double HYDROData_Lambert93::calc_rho( double phi )
135 double c1 = cot( phi/2 + PI/4 );
137 double c2 = ( 1 + e * sin( phi ) ) / ( 1 - e * sin( phi ) );
139 double rho = p_0 * pow( c1 * pow( c2, e/2 ), n );
147 void HYDROData_Lambert93::toXY( double theLatitudeDeg, double theLongitudeDeg,
149 double& theX, double& theY )
153 double phi = toRad( theLatitudeDeg );
155 double lambda = toRad( theLongitudeDeg );
159 double rho = calc_rho( phi );
161 double rho_0 = calc_rho( phi_0 );
163 double theta = n * ( lambda - lambda_0 );
167 theX = X_0 + rho * sin( theta );
169 theY = Y_0 + rho_0 - rho * cos( theta );
175 double arctan( double x )
185 typedef double (*FUNC)( double );
187 double solve( FUNC f, double c, double x1, double x2, double eps )
191 double f1 = f( x1 ) - c;
193 double f2 = f( x2 ) - c;
195 while( fabs( x1 - x2 ) > eps )
199 double x = ( x1 + x2 ) / 2;
201 double fx = f( x ) - c;
231 return ( x1 + x2 ) / 2;
237 double F( double phi )
241 double f1 = tan( phi/2 + PI/4 );
243 double f2 = ( 1 - e*sin(phi) ) / ( 1 + e*sin(phi) );
245 return f1 * pow( f2, e/2 );
251 double Finv( double x, double eps )
255 return solve( F, x, 0, PI/2-eps, eps );
261 double HYDROData_Lambert93::calc_phi_inv( double rho, double eps )
265 double x = pow( p_0 / rho, 1/n );
267 double phi = Finv( x, eps );
275 double HYDROData_Lambert93::calc_phi_ign( double rho, double eps )
279 double x = p_0 / rho;
281 double y = pow( x, 1/n );
283 double phi_i_1, phi_i = 2*arctan( y ) - PI/2;
291 double z = y * pow( ( 1 + e*sin(phi_i_1) ) / ( 1 - e*sin(phi_i_1) ), e/2 );
293 phi_i = 2*arctan( pow( x, 1/n ) * z ) - PI/2;
295 if( fabs( phi_i - phi_i_1 ) < eps )
307 void HYDROData_Lambert93::toGeo( double theX, double theY,
309 double& theLatitudeDeg, double& theLongitudeDeg,
315 double rho_0 = calc_rho( phi_0 );
317 double rho = sqrt( pow( theX - X_0, 2 ) + pow( Y_0 - theY + rho_0, 2 ) );
319 double theta = 2 * arctan( ( theX - X_0 ) / ( Y_0 - theY + rho_0 + rho ) );
323 double lambda = theta / n + lambda_0;
325 double phi = calc_phi_inv( rho, theEps );
329 theLatitudeDeg = toDeg( phi );
331 theLongitudeDeg = toDeg( lambda );
335 void HYDROData_Lambert93::DMSToDeg( int theDeg,
340 double aCoef = theDeg > 0 ? 1.0 : -1.0;
342 theDegOut = (double)theDeg;
343 theDegOut += aCoef * ( (double)theMin ) / 60.0;
344 theDegOut += aCoef * theSec / 3600.0;
347 void HYDROData_Lambert93::DMSToSec( int theDeg,
352 double aCoef = theDeg > 0 ? 1.0 : -1.0;
354 theSecOut = theDeg * 3600.0;
355 theSecOut += aCoef * theMin * 60.0;
356 theSecOut += aCoef * theSec;
359 void HYDROData_Lambert93::degToDMS( double theDegIn,
364 theDeg = int( theDegIn );
366 double aRemainder = fabs( theDegIn - theDeg ) * 60.0;
367 theMin = int( aRemainder );
369 aRemainder = ( aRemainder - theMin ) * 60.0;
373 void HYDROData_Lambert93::secToDMS( double theSecIn,
378 theDeg = int( theSecIn / 3600.0 );
380 double aRemainder = fabs( theSecIn - theDeg * 3600.0 );
381 theMin = int( aRemainder / 60.0 );
383 theSec = fabs( aRemainder - theMin * 60.0 );