src/HYDROData/HYDROData_Lambert93.cxx

   1 \r
   2 #include <HYDROData_Lambert93.h>\r
   3 #include <math.h>\r
   4 \r
   5 // Standard values\r
   6 const double PI = 3.14159265;\r
   7 \r
   8 // Base values of the Lambert-93\r
   9 const double a = 6378137;               //  m  -- le demi-grand axe\r
  10 const double f = 1.0 / 298.257222101;   // l'aplatissement\r
  11 \r
  12 const double phi_1_deg = 44;            // deg -- le premier parallèle d'échelle\r
  13 const double phi_2_deg = 49;            // deg -- le deuxième parallèle d'échell\r
  14 \r
  15 const double lambda_0_deg = 3;          // deg -- la longitude d'origine donnée par le méridien central de Greenwich\r
  16 const double phi_0_deg = 46.5;          // deg -- la latitude d'origine \r
  17 \r
  18 const double X_0 = 700000;              //  m  -- la coordonnée à l'origine\r
  19 const double Y_0 = 6600000;             //  m  -- la coordonnée à l'origine \r
  20 \r
  21 // Derived values of the Lambert-93\r
  22 const double b = a * ( 1 - f );         //  m  -- le demi-petit axe\r
  23 const double e = sqrt( a*a - b*b ) / a; // l'excentricité\r
  24 \r
  25 const double phi_0 = HYDROData_Lambert93::toRad( phi_0_deg );\r
  26 const double phi_1 = HYDROData_Lambert93::toRad( phi_1_deg );\r
  27 const double phi_2 = HYDROData_Lambert93::toRad( phi_2_deg );\r
  28 const double lambda_0 = HYDROData_Lambert93::toRad( lambda_0_deg );\r
  29 \r
  30 \r
  31 double cot( double x )\r
  32 {\r
  33   return cos( x ) / sin( x );\r
  34 }\r
  35 \r
  36 double ln( double x )\r
  37 {\r
  38   return log( x );\r
  39 }\r
  40 \r
  41 const double s1 = sin( phi_1 );\r
  42 const double s2 = sin( phi_2 );\r
  43 const double c1 = cos( phi_1 );\r
  44 const double c2 = cos( phi_2 );\r
  45 \r
  46 const double n1 = ln( c2/c1 ) + 1.0/2.0 * ln( (1-e*e*s1*s1)/(1-e*e*s2*s2) );\r
  47 const double n2 = tan( phi_1 / 2 + PI/4 ) * pow( 1-e*s1, e/2 ) * pow( 1+e*s2, e/2 );\r
  48 const double n3 = tan( phi_2 / 2 + PI/4 ) * pow( 1+e*s1, e/2 ) * pow( 1-e*s2, e/2 );\r
  49 const double n = n1 / ( ln( n2/n3 ) );\r
  50 \r
  51 const double p1 = a * c1 / ( n * sqrt( 1-e*e*s1*s1 ) );\r
  52 const double p2 = tan( phi_1 / 2 + PI / 4 );\r
  53 const double p3 = pow( (1-e*s1)/(1+e*s1), e/2 );\r
  54 const double p_0 = p1 * pow( p2*p3, n );\r
  55 \r
  56 double HYDROData_Lambert93::toRad( double theDeg )\r
  57 {\r
  58   return theDeg * PI / 180.0;\r
  59 }\r
  60 \r
  61 double HYDROData_Lambert93::toDeg( double theRad )\r
  62 {\r
  63   return theRad / PI * 180.0;\r
  64 }\r
  65 \r
  66 double HYDROData_Lambert93::calc_rho( double phi )\r
  67 {\r
  68   double c1 = cot( phi/2 + PI/4 );\r
  69   double c2 = ( 1 + e * sin( phi ) ) / ( 1 - e * sin( phi ) );\r
  70   double rho = p_0 * pow( c1 * pow( c2, e/2 ), n );\r
  71   return rho;\r
  72 }\r
  73 \r
  74 void HYDROData_Lambert93::toXY( double theLatitudeDeg, double theLongitudeDeg,\r
  75                                 double& theX, double& theY )\r
  76 {\r
  77   double phi = toRad( theLatitudeDeg );\r
  78   double lambda = toRad( theLongitudeDeg );\r
  79 \r
  80   double rho = calc_rho( phi );\r
  81   double rho_0 = calc_rho( phi_0 );\r
  82   double theta = n * ( lambda - lambda_0 );\r
  83 \r
  84   theX = X_0 + rho * sin( theta );\r
  85   theY = Y_0 + rho_0 - rho * cos( theta );\r
  86 }\r
  87 \r
  88 double arctan( double x )\r
  89 {\r
  90   return atan( x );\r
  91 }\r
  92 \r
  93 typedef double (*FUNC)( double );\r
  94 double solve( FUNC f, double c, double x1, double x2, double eps )\r
  95 {\r
  96   double f1 = f( x1 ) - c;\r
  97   double f2 = f( x2 ) - c;\r
  98   while( fabs( x1 - x2 ) > eps )\r
  99   {\r
 100     double x = ( x1 + x2 ) / 2;\r
 101     double fx = f( x ) - c;\r
 102     bool b1 = f1>=0;\r
 103     bool b2 = f2>=0;\r
 104     bool b = fx>=0;\r
 105     if( b==b1 )\r
 106     {\r
 107       x1 = x;\r
 108       f1 = fx;\r
 109     }\r
 110     else\r
 111     {\r
 112       x2 = x;\r
 113       f2 = fx;\r
 114     }\r
 115   }\r
 116   return ( x1 + x2 ) / 2;\r
 117 }\r
 118 \r
 119 double F( double phi )\r
 120 {\r
 121   double f1 = tan( phi/2 + PI/4 );\r
 122   double f2 = ( 1 - e*sin(phi) ) / ( 1 + e*sin(phi) );\r
 123   return f1 * pow( f2, e/2 );\r
 124 }\r
 125 \r
 126 double Finv( double x, double eps )\r
 127 {\r
 128   return solve( F, x, 0, PI/2-eps, eps );\r
 129 }\r
 130 \r
 131 double HYDROData_Lambert93::calc_phi_inv( double rho, double eps )\r
 132 {\r
 133   double x = pow( p_0 / rho, 1/n );\r
 134   double phi = Finv( x, eps );\r
 135   return phi;\r
 136 }\r
 137 \r
 138 double HYDROData_Lambert93::calc_phi_ign( double rho, double eps )\r
 139 {\r
 140   double x = p_0 / rho;\r
 141   double y = pow( x, 1/n );\r
 142   double phi_i_1, phi_i = 2*arctan( y ) - PI/2;\r
 143   while( true )\r
 144   {\r
 145     phi_i_1 = phi_i;\r
 146     double z = y * pow( ( 1 + e*sin(phi_i_1) ) / ( 1 - e*sin(phi_i_1) ), e/2 );\r
 147     phi_i = 2*arctan( pow( x, 1/n ) * z ) - PI/2;\r
 148     if( fabs( phi_i - phi_i_1 ) < eps )\r
 149       return phi_i;\r
 150   }\r
 151   return -1;\r
 152 }\r
 153 \r
 154 void HYDROData_Lambert93::toGeo( double theX, double theY,\r
 155                        double& theLatitudeDeg, double& theLongitudeDeg,\r
 156                        double theEps )\r
 157 {\r
 158   double rho_0 = calc_rho( phi_0 );\r
 159   double rho = sqrt( pow( theX - X_0, 2 ) + pow( Y_0 - theY + rho_0, 2 ) );\r
 160   double theta = 2 * arctan( ( theX - X_0 ) / ( Y_0 - theY + rho_0 + rho ) );\r
 161 \r
 162   double lambda = theta / n + lambda_0;\r
 163   double phi = calc_phi_inv( rho, theEps );\r
 164 \r
 165   theLatitudeDeg = toDeg( phi );\r
 166   theLongitudeDeg = toDeg( lambda );\r
 167 }\r