Merge remote-tracking branch 'origin/BR_IMPROVEMENTS' into BR_v14_rc

[modules/hydro.git] / src / HYDROData / HYDROData_Lambert93.cxx

// 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 );
}