-\r
-#include <HYDROData_Lambert93.h>\r
-#include <math.h>\r
-\r
-// Standard values\r
-const double PI = 3.14159265;\r
-\r
-// Base values of the Lambert-93\r
-const double a = 6378137; // m -- le demi-grand axe\r
-const double f = 1.0 / 298.257222101; // l'aplatissement\r
-\r
-const double phi_1_deg = 44; // deg -- le premier parallèle d'échelle\r
-const double phi_2_deg = 49; // deg -- le deuxième parallèle d'échell\r
-\r
-const double lambda_0_deg = 3; // deg -- la longitude d'origine donnée par le méridien central de Greenwich\r
-const double phi_0_deg = 46.5; // deg -- la latitude d'origine \r
-\r
-const double X_0 = 700000; // m -- la coordonnée à l'origine\r
-const double Y_0 = 6600000; // m -- la coordonnée à l'origine \r
-\r
-// Derived values of the Lambert-93\r
-const double b = a * ( 1 - f ); // m -- le demi-petit axe\r
-const double e = sqrt( a*a - b*b ) / a; // l'excentricité\r
-\r
-const double phi_0 = HYDROData_Lambert93::toRad( phi_0_deg );\r
-const double phi_1 = HYDROData_Lambert93::toRad( phi_1_deg );\r
-const double phi_2 = HYDROData_Lambert93::toRad( phi_2_deg );\r
-const double lambda_0 = HYDROData_Lambert93::toRad( lambda_0_deg );\r
-\r
-\r
-double cot( double x )\r
-{\r
- return cos( x ) / sin( x );\r
-}\r
-\r
-double ln( double x )\r
-{\r
- return log( x );\r
-}\r
-\r
-const double s1 = sin( phi_1 );\r
-const double s2 = sin( phi_2 );\r
-const double c1 = cos( phi_1 );\r
-const double c2 = cos( phi_2 );\r
-\r
-const double n1 = ln( c2/c1 ) + 1.0/2.0 * ln( (1-e*e*s1*s1)/(1-e*e*s2*s2) );\r
-const double n2 = tan( phi_1 / 2 + PI/4 ) * pow( 1-e*s1, e/2 ) * pow( 1+e*s2, e/2 );\r
-const double n3 = tan( phi_2 / 2 + PI/4 ) * pow( 1+e*s1, e/2 ) * pow( 1-e*s2, e/2 );\r
-const double n = n1 / ( ln( n2/n3 ) );\r
-\r
-const double p1 = a * c1 / ( n * sqrt( 1-e*e*s1*s1 ) );\r
-const double p2 = tan( phi_1 / 2 + PI / 4 );\r
-const double p3 = pow( (1-e*s1)/(1+e*s1), e/2 );\r
-const double p_0 = p1 * pow( p2*p3, n );\r
-\r
-double HYDROData_Lambert93::toRad( double theDeg )\r
-{\r
- return theDeg * PI / 180.0;\r
-}\r
-\r
-double HYDROData_Lambert93::toDeg( double theRad )\r
-{\r
- return theRad / PI * 180.0;\r
-}\r
-\r
-double HYDROData_Lambert93::calc_rho( double phi )\r
-{\r
- double c1 = cot( phi/2 + PI/4 );\r
- double c2 = ( 1 + e * sin( phi ) ) / ( 1 - e * sin( phi ) );\r
- double rho = p_0 * pow( c1 * pow( c2, e/2 ), n );\r
- return rho;\r
-}\r
-\r
-void HYDROData_Lambert93::toXY( double theLatitudeDeg, double theLongitudeDeg,\r
- double& theX, double& theY )\r
-{\r
- double phi = toRad( theLatitudeDeg );\r
- double lambda = toRad( theLongitudeDeg );\r
-\r
- double rho = calc_rho( phi );\r
- double rho_0 = calc_rho( phi_0 );\r
- double theta = n * ( lambda - lambda_0 );\r
-\r
- theX = X_0 + rho * sin( theta );\r
- theY = Y_0 + rho_0 - rho * cos( theta );\r
-}\r
-\r
-double arctan( double x )\r
-{\r
- return atan( x );\r
-}\r
-\r
-typedef double (*FUNC)( double );\r
-double solve( FUNC f, double c, double x1, double x2, double eps )\r
-{\r
- double f1 = f( x1 ) - c;\r
- double f2 = f( x2 ) - c;\r
- while( fabs( x1 - x2 ) > eps )\r
- {\r
- double x = ( x1 + x2 ) / 2;\r
- double fx = f( x ) - c;\r
- bool b1 = f1>=0;\r
- bool b2 = f2>=0;\r
- bool b = fx>=0;\r
- if( b==b1 )\r
- {\r
- x1 = x;\r
- f1 = fx;\r
- }\r
- else\r
- {\r
- x2 = x;\r
- f2 = fx;\r
- }\r
- }\r
- return ( x1 + x2 ) / 2;\r
-}\r
-\r
-double F( double phi )\r
-{\r
- double f1 = tan( phi/2 + PI/4 );\r
- double f2 = ( 1 - e*sin(phi) ) / ( 1 + e*sin(phi) );\r
- return f1 * pow( f2, e/2 );\r
-}\r
-\r
-double Finv( double x, double eps )\r
-{\r
- return solve( F, x, 0, PI/2-eps, eps );\r
-}\r
-\r
-double HYDROData_Lambert93::calc_phi_inv( double rho, double eps )\r
-{\r
- double x = pow( p_0 / rho, 1/n );\r
- double phi = Finv( x, eps );\r
- return phi;\r
-}\r
-\r
-double HYDROData_Lambert93::calc_phi_ign( double rho, double eps )\r
-{\r
- double x = p_0 / rho;\r
- double y = pow( x, 1/n );\r
- double phi_i_1, phi_i = 2*arctan( y ) - PI/2;\r
- while( true )\r
- {\r
- phi_i_1 = phi_i;\r
- double z = y * pow( ( 1 + e*sin(phi_i_1) ) / ( 1 - e*sin(phi_i_1) ), e/2 );\r
- phi_i = 2*arctan( pow( x, 1/n ) * z ) - PI/2;\r
- if( fabs( phi_i - phi_i_1 ) < eps )\r
- return phi_i;\r
- }\r
- return -1;\r
-}\r
-\r
-void HYDROData_Lambert93::toGeo( double theX, double theY,\r
- double& theLatitudeDeg, double& theLongitudeDeg,\r
- double theEps )\r
-{\r
- double rho_0 = calc_rho( phi_0 );\r
- double rho = sqrt( pow( theX - X_0, 2 ) + pow( Y_0 - theY + rho_0, 2 ) );\r
- double theta = 2 * arctan( ( theX - X_0 ) / ( Y_0 - theY + rho_0 + rho ) );\r
-\r
- double lambda = theta / n + lambda_0;\r
- double phi = calc_phi_inv( rho, theEps );\r
-\r
- theLatitudeDeg = toDeg( phi );\r
- theLongitudeDeg = toDeg( lambda );\r
-}\r
+// Copyright (C) 2014-2015 EDF-R&D
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+#include <HYDROData_Lambert93.h>
+#include <math.h>
+
+// Standard values
+const double PI = 3.14159265;
+
+// Base values of the Lambert-93
+const double a = 6378137; // m -- le demi-grand axe
+const double f = 1.0 / 298.257222101; // l'aplatissement
+const double phi_1_deg = 44; // deg -- le premier parallèle d'échelle
+const double phi_2_deg = 49; // deg -- le deuxième parallèle d'échell
+const double lambda_0_deg = 3; // deg -- la longitude d'origine donnée par le méridien central de Greenwich
+const double phi_0_deg = 46.5; // deg -- la latitude d'origine
+const double X_0 = 700000; // m -- la coordonnée à l'origine
+const double Y_0 = 6600000; // m -- la coordonnée à l'origine
+
+// Derived values of the Lambert-93
+const double b = a * ( 1 - f ); // m -- le demi-petit axe
+const double e = sqrt( a*a - b*b ) / a; // l'excentricité
+
+const double phi_0 = HYDROData_Lambert93::toRad( phi_0_deg );
+const double phi_1 = HYDROData_Lambert93::toRad( phi_1_deg );
+const double phi_2 = HYDROData_Lambert93::toRad( phi_2_deg );
+const double lambda_0 = HYDROData_Lambert93::toRad( lambda_0_deg );
+
+double cot( double x )
+{
+ return cos( x ) / sin( x );
+}
+
+double ln( double x )
+{
+ return log( x );
+}
+
+const double s1 = sin( phi_1 );
+const double s2 = sin( phi_2 );
+const double c1 = cos( phi_1 );
+const double c2 = cos( phi_2 );
+
+const double n1 = ln( c2/c1 ) + 1.0/2.0 * ln( (1-e*e*s1*s1)/(1-e*e*s2*s2) );
+const double n2 = tan( phi_1 / 2 + PI/4 ) * pow( 1-e*s1, e/2 ) * pow( 1+e*s2, e/2 );
+const double n3 = tan( phi_2 / 2 + PI/4 ) * pow( 1+e*s1, e/2 ) * pow( 1-e*s2, e/2 );
+const double n = n1 / ( ln( n2/n3 ) );
+
+const double p1 = a * c1 / ( n * sqrt( 1-e*e*s1*s1 ) );
+const double p2 = tan( phi_1 / 2 + PI / 4 );
+const double p3 = pow( (1-e*s1)/(1+e*s1), e/2 );
+const double p_0 = p1 * pow( p2*p3, n );
+
+double HYDROData_Lambert93::toRad( double theDeg )
+{
+ return theDeg * PI / 180.0;
+}
+
+double HYDROData_Lambert93::toDeg( double theRad )
+{
+ return theRad / PI * 180.0;
+}
+
+double HYDROData_Lambert93::calc_rho( double phi )
+{
+ double c1 = cot( phi/2 + PI/4 );
+ double c2 = ( 1 + e * sin( phi ) ) / ( 1 - e * sin( phi ) );
+ double rho = p_0 * pow( c1 * pow( c2, e/2 ), n );
+ return rho;
+}
+
+void HYDROData_Lambert93::toXY( double theLatitudeDeg, double theLongitudeDeg,
+ double& theX, double& theY )
+{
+ double phi = toRad( theLatitudeDeg );
+ double lambda = toRad( theLongitudeDeg );
+
+ double rho = calc_rho( phi );
+ double rho_0 = calc_rho( phi_0 );
+ double theta = n * ( lambda - lambda_0 );
+
+ theX = X_0 + rho * sin( theta );
+ theY = Y_0 + rho_0 - rho * cos( theta );
+}
+
+double arctan( double x )
+{
+ return atan( x );
+}
+
+typedef double (*FUNC)( double );
+
+double solve( FUNC f, double c, double x1, double x2, double eps )
+{
+ double f1 = f( x1 ) - c;
+ double f2 = f( x2 ) - c;
+ while( fabs( x1 - x2 ) > eps )
+ {
+ double x = ( x1 + x2 ) / 2;
+ double fx = f( x ) - c;
+ bool b1 = f1>=0;
+ bool b2 = f2>=0;
+ bool b = fx>=0;
+ if( b==b1 )
+ {
+ x1 = x;
+ f1 = fx;
+ }
+ else
+ {
+ x2 = x;
+ f2 = fx;
+ }
+ }
+ return ( x1 + x2 ) / 2;
+}
+
+double F( double phi )
+{
+ double f1 = tan( phi/2 + PI/4 );
+ double f2 = ( 1 - e*sin(phi) ) / ( 1 + e*sin(phi) );
+ return f1 * pow( f2, e/2 );
+}
+
+double Finv( double x, double eps )
+{
+ return solve( F, x, 0, PI/2-eps, eps );
+}
+
+double HYDROData_Lambert93::calc_phi_inv( double rho, double eps )
+{
+ double x = pow( p_0 / rho, 1/n );
+ double phi = Finv( x, eps );
+ return phi;
+}
+
+double HYDROData_Lambert93::calc_phi_ign( double rho, double eps )
+{
+ double x = p_0 / rho;
+ double y = pow( x, 1/n );
+ double phi_i_1, phi_i = 2*arctan( y ) - PI/2;
+ while( true )
+ {
+ phi_i_1 = phi_i;
+ double z = y * pow( ( 1 + e*sin(phi_i_1) ) / ( 1 - e*sin(phi_i_1) ), e/2 );
+ phi_i = 2*arctan( pow( x, 1/n ) * z ) - PI/2;
+ if( fabs( phi_i - phi_i_1 ) < eps )
+ return phi_i;
+ }
+ return -1;
+}
+
+void HYDROData_Lambert93::toGeo( double theX, double theY,
+ double& theLatitudeDeg, double& theLongitudeDeg,
+ double theEps )
+{
+ double rho_0 = calc_rho( phi_0 );
+ double rho = sqrt( pow( theX - X_0, 2 ) + pow( Y_0 - theY + rho_0, 2 ) );
+ double theta = 2 * arctan( ( theX - X_0 ) / ( Y_0 - theY + rho_0 + rho ) );
+
+ double lambda = theta / n + lambda_0;
+ double phi = calc_phi_inv( rho, theEps );
+
+ theLatitudeDeg = toDeg( phi );
+ theLongitudeDeg = toDeg( lambda );
+}
+
+void HYDROData_Lambert93::DMSToDeg( int theDeg,
+ int theMin,
+ double theSec,
+ double& theDegOut )
+{
+ double aCoef = theDeg > 0 ? 1.0 : -1.0;
+
+ theDegOut = (double)theDeg;
+ theDegOut += aCoef * ( (double)theMin ) / 60.0;
+ theDegOut += aCoef * theSec / 3600.0;
+}
+
+void HYDROData_Lambert93::DMSToSec( int theDeg,
+ int theMin,
+ double theSec,
+ double& theSecOut )
+{
+ double aCoef = theDeg > 0 ? 1.0 : -1.0;
+
+ theSecOut = theDeg * 3600.0;
+ theSecOut += aCoef * theMin * 60.0;
+ theSecOut += aCoef * theSec;
+}
+
+void HYDROData_Lambert93::degToDMS( double theDegIn,
+ int& theDeg,
+ int& theMin,
+ double& theSec )
+{
+ theDeg = int( theDegIn );
+
+ double aRemainder = fabs( theDegIn - theDeg ) * 60.0;
+ theMin = int( aRemainder );
+
+ aRemainder = ( aRemainder - theMin ) * 60.0;
+ theSec = aRemainder;
+}
+
+void HYDROData_Lambert93::secToDMS( double theSecIn,
+ int& theDeg,
+ int& theMin,
+ double& theSec )
+{
+ theDeg = int( theSecIn / 3600.0 );
+
+ double aRemainder = fabs( theSecIn - theDeg * 3600.0 );
+ theMin = int( aRemainder / 60.0 );
+
+ theSec = fabs( aRemainder - theMin * 60.0 );
+}
