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 #include <math.h>
  21
  22 // Standard values
  23 const double PI = 3.14159265;
  24
  25 // Base values of the Lambert-93
  26 const double a = 6378137;               //  m  -- le demi-grand axe
  27 const double f = 1.0 / 298.257222101;   // l'aplatissement
  28 const double phi_1_deg = 44;            // deg -- le premier parallèle d'échelle
  29 const double phi_2_deg = 49;            // deg -- le deuxième parallèle d'échell
  30 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
  32 const double X_0 = 700000;              //  m  -- la coordonnée à l'origine
  33 const double Y_0 = 6600000;             //  m  -- la coordonnée à l'origine
  34
  35 // Derived values of the Lambert-93
  36 const double b = a * ( 1 - f );         //  m  -- le demi-petit axe
  37 const double e = sqrt( a*a - b*b ) / a; // l'excentricité
  38
  39 const double phi_0 = HYDROData_Lambert93::toRad( phi_0_deg );
  40 const double phi_1 = HYDROData_Lambert93::toRad( phi_1_deg );
  41 const double phi_2 = HYDROData_Lambert93::toRad( phi_2_deg );
  42 const double lambda_0 = HYDROData_Lambert93::toRad( lambda_0_deg );
  43
  44 double cot( double x )
  45 {
  46   return cos( x ) / sin( x );
  47 }
  48
  49 double ln( double x )
  50 {
  51   return log( x );
  52 }
  53
  54 const double s1 = sin( phi_1 );
  55 const double s2 = sin( phi_2 );
  56 const double c1 = cos( phi_1 );
  57 const double c2 = cos( phi_2 );
  58
  59 const double n1 = ln( c2/c1 ) + 1.0/2.0 * ln( (1-e*e*s1*s1)/(1-e*e*s2*s2) );
  60 const double n2 = tan( phi_1 / 2 + PI/4 ) * pow( 1-e*s1, e/2 ) * pow( 1+e*s2, e/2 );
  61 const double n3 = tan( phi_2 / 2 + PI/4 ) * pow( 1+e*s1, e/2 ) * pow( 1-e*s2, e/2 );
  62 const double n = n1 / ( ln( n2/n3 ) );
  63
  64 const double p1 = a * c1 / ( n * sqrt( 1-e*e*s1*s1 ) );
  65 const double p2 = tan( phi_1 / 2 + PI / 4 );
  66 const double p3 = pow( (1-e*s1)/(1+e*s1), e/2 );
  67 const double p_0 = p1 * pow( p2*p3, n );
  68
  69 double HYDROData_Lambert93::toRad( double theDeg )
  70 {
  71   return theDeg * PI / 180.0;
  72 }
  73
  74 double HYDROData_Lambert93::toDeg( double theRad )
  75 {
  76   return theRad / PI * 180.0;
  77 }
  78
  79 double HYDROData_Lambert93::calc_rho( double phi )
  80 {
  81   double c1 = cot( phi/2 + PI/4 );
  82   double c2 = ( 1 + e * sin( phi ) ) / ( 1 - e * sin( phi ) );
  83   double rho = p_0 * pow( c1 * pow( c2, e/2 ), n );
  84   return rho;
  85 }
  86
  87 void HYDROData_Lambert93::toXY( double theLatitudeDeg, double theLongitudeDeg,
  88                                 double& theX, double& theY )
  89 {
  90   double phi = toRad( theLatitudeDeg );
  91   double lambda = toRad( theLongitudeDeg );
  92
  93   double rho = calc_rho( phi );
  94   double rho_0 = calc_rho( phi_0 );
  95   double theta = n * ( lambda - lambda_0 );
  96
  97   theX = X_0 + rho * sin( theta );
  98   theY = Y_0 + rho_0 - rho * cos( theta );
  99 }
 100
 101 double arctan( double x )
 102 {
 103   return atan( x );
 104 }
 105
 106 typedef double (*FUNC)( double );
 107
 108 double solve( FUNC f, double c, double x1, double x2, double eps )
 109 {
 110   double f1 = f( x1 ) - c;
 111   double f2 = f( x2 ) - c;
 112   while( fabs( x1 - x2 ) > eps )
 113   {
 114     double x = ( x1 + x2 ) / 2;
 115     double fx = f( x ) - c;
 116     bool b1 = f1>=0;
 117     bool b2 = f2>=0;
 118     bool b = fx>=0;
 119     if( b==b1 )
 120     {
 121       x1 = x;
 122       f1 = fx;
 123     }
 124     else
 125     {
 126       x2 = x;
 127       f2 = fx;
 128     }
 129   }
 130   return ( x1 + x2 ) / 2;
 131 }
 132
 133 double F( double phi )
 134 {
 135   double f1 = tan( phi/2 + PI/4 );
 136   double f2 = ( 1 - e*sin(phi) ) / ( 1 + e*sin(phi) );
 137   return f1 * pow( f2, e/2 );
 138 }
 139
 140 double Finv( double x, double eps )
 141 {
 142   return solve( F, x, 0, PI/2-eps, eps );
 143 }
 144
 145 double HYDROData_Lambert93::calc_phi_inv( double rho, double eps )
 146 {
 147   double x = pow( p_0 / rho, 1/n );
 148   double phi = Finv( x, eps );
 149   return phi;
 150 }
 151
 152 double HYDROData_Lambert93::calc_phi_ign( double rho, double eps )
 153 {
 154   double x = p_0 / rho;
 155   double y = pow( x, 1/n );
 156   double phi_i_1, phi_i = 2*arctan( y ) - PI/2;
 157   while( true )
 158   {
 159     phi_i_1 = phi_i;
 160     double z = y * pow( ( 1 + e*sin(phi_i_1) ) / ( 1 - e*sin(phi_i_1) ), e/2 );
 161     phi_i = 2*arctan( pow( x, 1/n ) * z ) - PI/2;
 162     if( fabs( phi_i - phi_i_1 ) < eps )
 163       return phi_i;
 164   }
 165   return -1;
 166 }
 167
 168 void HYDROData_Lambert93::toGeo( double theX, double theY,
 169                                  double& theLatitudeDeg, double& theLongitudeDeg,
 170                                  double theEps )
 171 {
 172   double rho_0 = calc_rho( phi_0 );
 173   double rho = sqrt( pow( theX - X_0, 2 ) + pow( Y_0 - theY + rho_0, 2 ) );
 174   double theta = 2 * arctan( ( theX - X_0 ) / ( Y_0 - theY + rho_0 + rho ) );
 175
 176   double lambda = theta / n + lambda_0;
 177   double phi = calc_phi_inv( rho, theEps );
 178
 179   theLatitudeDeg = toDeg( phi );
 180   theLongitudeDeg = toDeg( lambda );
 181 }
 182
 183 void HYDROData_Lambert93::DMSToDeg( int theDeg,
 184                                     int theMin,
 185                                     double theSec,
 186                                     double& theDegOut )
 187 {
 188   double aCoef = theDeg > 0 ? 1.0 : -1.0;
 189
 190   theDegOut = (double)theDeg;
 191   theDegOut += aCoef * ( (double)theMin ) / 60.0;
 192   theDegOut += aCoef * theSec / 3600.0;
 193 }
 194
 195 void HYDROData_Lambert93::DMSToSec( int theDeg,
 196                                     int theMin,
 197                                     double theSec,
 198                                     double& theSecOut )
 199 {
 200   double aCoef = theDeg > 0 ? 1.0 : -1.0;
 201
 202   theSecOut = theDeg * 3600.0;
 203   theSecOut += aCoef * theMin * 60.0;
 204   theSecOut += aCoef * theSec;
 205 }
 206
 207 void HYDROData_Lambert93::degToDMS( double theDegIn,
 208                                     int& theDeg,
 209                                     int& theMin,
 210                                     double& theSec )
 211 {
 212   theDeg = int( theDegIn );
 213
 214   double aRemainder = fabs( theDegIn - theDeg ) * 60.0;
 215   theMin = int( aRemainder );
 216
 217   aRemainder = ( aRemainder - theMin ) * 60.0;
 218   theSec = aRemainder;
 219 }
 220
 221 void HYDROData_Lambert93::secToDMS( double theSecIn,
 222                                     int& theDeg,
 223                                     int& theMin,
 224                                     double& theSec )
 225 {
 226   theDeg = int( theSecIn / 3600.0 );
 227
 228   double aRemainder = fabs( theSecIn - theDeg * 3600.0 );
 229   theMin = int( aRemainder / 60.0 );
 230
 231   theSec = fabs( aRemainder - theMin * 60.0 );
 232 }