src/HYDROData/HYDROData_Lambert93.cxx

   1
   2
   3 #include <HYDROData_Lambert93.h>
   4
   5 #include <math.h>
   6
   7
   8
   9 // Standard values
  10
  11 const double PI = 3.14159265;
  12
  13
  14
  15 // Base values of the Lambert-93
  16
  17 const double a = 6378137;               //  m  -- le demi-grand axe
  18
  19 const double f = 1.0 / 298.257222101;   // l'aplatissement
  20
  21
  22
  23 const double phi_1_deg = 44;            // deg -- le premier parallèle d'échelle
  24
  25 const double phi_2_deg = 49;            // deg -- le deuxième parallèle d'échell
  26
  27
  28
  29 const double lambda_0_deg = 3;          // deg -- la longitude d'origine donnée par le méridien central de Greenwich
  30
  31 const double phi_0_deg = 46.5;          // deg -- la latitude d'origine
  32
  33
  34
  35 const double X_0 = 700000;              //  m  -- la coordonnée à l'origine
  36
  37 const double Y_0 = 6600000;             //  m  -- la coordonnée à l'origine
  38
  39
  40
  41 // Derived values of the Lambert-93
  42
  43 const double b = a * ( 1 - f );         //  m  -- le demi-petit axe
  44
  45 const double e = sqrt( a*a - b*b ) / a; // l'excentricité
  46
  47
  48
  49 const double phi_0 = HYDROData_Lambert93::toRad( phi_0_deg );
  50
  51 const double phi_1 = HYDROData_Lambert93::toRad( phi_1_deg );
  52
  53 const double phi_2 = HYDROData_Lambert93::toRad( phi_2_deg );
  54
  55 const double lambda_0 = HYDROData_Lambert93::toRad( lambda_0_deg );
  56
  57
  58
  59
  60
  61 double cot( double x )
  62
  63 {
  64
  65   return cos( x ) / sin( x );
  66
  67 }
  68
  69
  70
  71 double ln( double x )
  72
  73 {
  74
  75   return log( x );
  76
  77 }
  78
  79
  80
  81 const double s1 = sin( phi_1 );
  82
  83 const double s2 = sin( phi_2 );
  84
  85 const double c1 = cos( phi_1 );
  86
  87 const double c2 = cos( phi_2 );
  88
  89
  90
  91 const double n1 = ln( c2/c1 ) + 1.0/2.0 * ln( (1-e*e*s1*s1)/(1-e*e*s2*s2) );
  92
  93 const double n2 = tan( phi_1 / 2 + PI/4 ) * pow( 1-e*s1, e/2 ) * pow( 1+e*s2, e/2 );
  94
  95 const double n3 = tan( phi_2 / 2 + PI/4 ) * pow( 1+e*s1, e/2 ) * pow( 1-e*s2, e/2 );
  96
  97 const double n = n1 / ( ln( n2/n3 ) );
  98
  99
 100
 101 const double p1 = a * c1 / ( n * sqrt( 1-e*e*s1*s1 ) );
 102
 103 const double p2 = tan( phi_1 / 2 + PI / 4 );
 104
 105 const double p3 = pow( (1-e*s1)/(1+e*s1), e/2 );
 106
 107 const double p_0 = p1 * pow( p2*p3, n );
 108
 109
 110
 111 double HYDROData_Lambert93::toRad( double theDeg )
 112
 113 {
 114
 115   return theDeg * PI / 180.0;
 116
 117 }
 118
 119
 120
 121 double HYDROData_Lambert93::toDeg( double theRad )
 122
 123 {
 124
 125   return theRad / PI * 180.0;
 126
 127 }
 128
 129
 130
 131 double HYDROData_Lambert93::calc_rho( double phi )
 132
 133 {
 134
 135   double c1 = cot( phi/2 + PI/4 );
 136
 137   double c2 = ( 1 + e * sin( phi ) ) / ( 1 - e * sin( phi ) );
 138
 139   double rho = p_0 * pow( c1 * pow( c2, e/2 ), n );
 140
 141   return rho;
 142
 143 }
 144
 145
 146
 147 void HYDROData_Lambert93::toXY( double theLatitudeDeg, double theLongitudeDeg,
 148
 149                                 double& theX, double& theY )
 150
 151 {
 152
 153   double phi = toRad( theLatitudeDeg );
 154
 155   double lambda = toRad( theLongitudeDeg );
 156
 157
 158
 159   double rho = calc_rho( phi );
 160
 161   double rho_0 = calc_rho( phi_0 );
 162
 163   double theta = n * ( lambda - lambda_0 );
 164
 165
 166
 167   theX = X_0 + rho * sin( theta );
 168
 169   theY = Y_0 + rho_0 - rho * cos( theta );
 170
 171 }
 172
 173
 174
 175 double arctan( double x )
 176
 177 {
 178
 179   return atan( x );
 180
 181 }
 182
 183
 184
 185 typedef double (*FUNC)( double );
 186
 187 double solve( FUNC f, double c, double x1, double x2, double eps )
 188
 189 {
 190
 191   double f1 = f( x1 ) - c;
 192
 193   double f2 = f( x2 ) - c;
 194
 195   while( fabs( x1 - x2 ) > eps )
 196
 197   {
 198
 199     double x = ( x1 + x2 ) / 2;
 200
 201     double fx = f( x ) - c;
 202
 203     bool b1 = f1>=0;
 204
 205     bool b2 = f2>=0;
 206
 207     bool b = fx>=0;
 208
 209     if( b==b1 )
 210
 211     {
 212
 213       x1 = x;
 214
 215       f1 = fx;
 216
 217     }
 218
 219     else
 220
 221     {
 222
 223       x2 = x;
 224
 225       f2 = fx;
 226
 227     }
 228
 229   }
 230
 231   return ( x1 + x2 ) / 2;
 232
 233 }
 234
 235
 236
 237 double F( double phi )
 238
 239 {
 240
 241   double f1 = tan( phi/2 + PI/4 );
 242
 243   double f2 = ( 1 - e*sin(phi) ) / ( 1 + e*sin(phi) );
 244
 245   return f1 * pow( f2, e/2 );
 246
 247 }
 248
 249
 250
 251 double Finv( double x, double eps )
 252
 253 {
 254
 255   return solve( F, x, 0, PI/2-eps, eps );
 256
 257 }
 258
 259
 260
 261 double HYDROData_Lambert93::calc_phi_inv( double rho, double eps )
 262
 263 {
 264
 265   double x = pow( p_0 / rho, 1/n );
 266
 267   double phi = Finv( x, eps );
 268
 269   return phi;
 270
 271 }
 272
 273
 274
 275 double HYDROData_Lambert93::calc_phi_ign( double rho, double eps )
 276
 277 {
 278
 279   double x = p_0 / rho;
 280
 281   double y = pow( x, 1/n );
 282
 283   double phi_i_1, phi_i = 2*arctan( y ) - PI/2;
 284
 285   while( true )
 286
 287   {
 288
 289     phi_i_1 = phi_i;
 290
 291     double z = y * pow( ( 1 + e*sin(phi_i_1) ) / ( 1 - e*sin(phi_i_1) ), e/2 );
 292
 293     phi_i = 2*arctan( pow( x, 1/n ) * z ) - PI/2;
 294
 295     if( fabs( phi_i - phi_i_1 ) < eps )
 296
 297       return phi_i;
 298
 299   }
 300
 301   return -1;
 302
 303 }
 304
 305
 306
 307 void HYDROData_Lambert93::toGeo( double theX, double theY,
 308
 309                        double& theLatitudeDeg, double& theLongitudeDeg,
 310
 311                        double theEps )
 312
 313 {
 314
 315   double rho_0 = calc_rho( phi_0 );
 316
 317   double rho = sqrt( pow( theX - X_0, 2 ) + pow( Y_0 - theY + rho_0, 2 ) );
 318
 319   double theta = 2 * arctan( ( theX - X_0 ) / ( Y_0 - theY + rho_0 + rho ) );
 320
 321
 322
 323   double lambda = theta / n + lambda_0;
 324
 325   double phi = calc_phi_inv( rho, theEps );
 326
 327
 328
 329   theLatitudeDeg = toDeg( phi );
 330
 331   theLongitudeDeg = toDeg( lambda );
 332
 333 }
 334
 335 void HYDROData_Lambert93::DMSToDeg( int theDeg,
 336                                     int theMin,
 337                                     double theSec,
 338                                     double& theDegOut )
 339 {
 340   double aCoef = theDeg > 0 ? 1.0 : -1.0;
 341
 342   theDegOut = (double)theDeg;
 343   theDegOut += aCoef * ( (double)theMin ) / 60.0;
 344   theDegOut += aCoef * theSec / 3600.0;
 345 }
 346
 347 void HYDROData_Lambert93::DMSToSec( int theDeg,
 348                                     int theMin,
 349                                     double theSec,
 350                                     double& theSecOut )
 351 {
 352   double aCoef = theDeg > 0 ? 1.0 : -1.0;
 353
 354   theSecOut = theDeg * 3600.0;
 355   theSecOut += aCoef * theMin * 60.0;
 356   theSecOut += aCoef * theSec;
 357 }
 358
 359 void HYDROData_Lambert93::degToDMS( double theDegIn,
 360                                     int& theDeg,
 361                                     int& theMin,
 362                                     double& theSec )
 363 {
 364   theDeg = int( theDegIn );
 365
 366   double aRemainder = fabs( theDegIn - theDeg ) * 60.0;
 367   theMin = int( aRemainder );
 368
 369   aRemainder = ( aRemainder - theMin ) * 60.0;
 370   theSec = aRemainder;
 371 }
 372
 373 void HYDROData_Lambert93::secToDMS( double theSecIn,
 374                                     int& theDeg,
 375                                     int& theMin,
 376                                     double& theSec )
 377 {
 378   theDeg = int( theSecIn / 3600.0 );
 379
 380   double aRemainder = fabs( theSecIn - theDeg * 3600.0 );
 381   theMin = int( aRemainder / 60.0 );
 382
 383   theSec = fabs( aRemainder - theMin * 60.0 );
 384 }