X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FHYDROData%2FHYDROData_Lambert93.cxx;h=c7cc478fce7434a89d412c8b1f567511d676a8f1;hb=148fc44cfb9114849974f209db4c8596e0d97507;hp=786eb7e6ce1ae687af02221ef6616d09595e04e0;hpb=7e60cb9b6e6e46019c188a935a97ef0d61d12abd;p=modules%2Fhydro.git diff --git a/src/HYDROData/HYDROData_Lambert93.cxx b/src/HYDROData/HYDROData_Lambert93.cxx old mode 100755 new mode 100644 index 786eb7e6..c7cc478f --- a/src/HYDROData/HYDROData_Lambert93.cxx +++ b/src/HYDROData/HYDROData_Lambert93.cxx @@ -1,167 +1,232 @@ - -#include -#include - -// 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 ); -} +// 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 +#include + +// 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 ); +}