src/HYDROData/HYDROData_Lambert93.cxx

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