]> git.mxchange.org Git - simgear.git/blobdiff - simgear/math/sg_geodesy.cxx
Don't bother other develoers with problems caused by MipsPro (version < 7.3) compilers
[simgear.git] / simgear / math / sg_geodesy.cxx
index 3981530885f7e66e07ecaa38dbc00a408c964c88..186907b2febc5b6bc9f747e7291f1f50548aa0b9 100644 (file)
@@ -9,7 +9,7 @@
 
 #include <simgear/compiler.h>
 
-#ifdef FG_HAVE_STD_INCLUDES
+#ifdef SG_HAVE_STD_INCLUDES
 # include <cmath>
 # include <cerrno>
 # include <cstdio>
 #include "localconsts.hxx"
 
 
-#ifndef FG_HAVE_NATIVE_SGI_COMPILERS
-FG_USING_STD(cout);
-#endif
-
-// ONE_SECOND is pi/180/60/60, or about 100 feet at earths' equator
-#define ONE_SECOND 4.848136811E-6
+SG_USING_STD(cout);
 
 
-#define DOMAIN_ERR_DEBUG 1
+// #define DOMAIN_ERR_DEBUG 1
 
 
 // sgGeocToGeod(lat_geoc, radius, *lat_geod, *alt, *sea_level_r)
@@ -50,8 +45,8 @@ FG_USING_STD(cout);
 //                      local vertical (surface normal) of C.G. (meters)
 
 
-void sgGeocToGeod( double lat_geoc, double radius, double
-                  *lat_geod, double *alt, double *sea_level_r )
+void sgGeocToGeod( const double& lat_geoc, const double& radius,
+                   double *lat_geod, double *alt, double *sea_level_r )
 {
 #ifdef DOMAIN_ERR_DEBUG
     errno = 0;                 // start with error zero'd
@@ -60,30 +55,30 @@ void sgGeocToGeod( double lat_geoc, double radius, double
     double t_lat, x_alpha, mu_alpha, delt_mu, r_alpha, l_point, rho_alpha;
     double sin_mu_a, denom,delt_lambda, lambda_sl, sin_lambda_sl;
 
-    if( ( (FG_PI_2 - lat_geoc) < ONE_SECOND )        // near North pole
-       || ( (FG_PI_2 + lat_geoc) < ONE_SECOND ) )   // near South pole
+    if( ( (SGD_PI_2 - lat_geoc) < SG_ONE_SECOND )        // near North pole
+       || ( (SGD_PI_2 + lat_geoc) < SG_ONE_SECOND ) )   // near South pole
     {
        *lat_geod = lat_geoc;
-       *sea_level_r = EQUATORIAL_RADIUS_M*E;
+       *sea_level_r = SG_EQUATORIAL_RADIUS_M*E;
        *alt = radius - *sea_level_r;
     } else {
        // cout << "  lat_geoc = " << lat_geoc << endl;
        t_lat = tan(lat_geoc);
        // cout << "  tan(t_lat) = " << t_lat << endl;
-       x_alpha = E*EQUATORIAL_RADIUS_M/sqrt(t_lat*t_lat + E*E);
+       x_alpha = E*SG_EQUATORIAL_RADIUS_M/sqrt(t_lat*t_lat + E*E);
 #ifdef DOMAIN_ERR_DEBUG
        if ( errno ) {
            perror("fgGeocToGeod()");
-           FG_LOG( FG_GENERAL, FG_ALERT, "sqrt(" << t_lat*t_lat + E*E << ")" );
+           SG_LOG( SG_GENERAL, SG_ALERT, "sqrt(" << t_lat*t_lat + E*E << ")" );
        }
 #endif
        // cout << "  x_alpha = " << x_alpha << endl;
-       double tmp = sqrt(RESQ_M - x_alpha * x_alpha);
+       double tmp = sqrt(SG_EQ_RAD_SQUARE_M - x_alpha * x_alpha);
        if ( tmp < 0.0 ) { tmp = 0.0; }
 #ifdef DOMAIN_ERR_DEBUG
        if ( errno ) {
            perror("fgGeocToGeod()");
-           FG_LOG( FG_GENERAL, FG_ALERT, "sqrt(" << RESQ_M - x_alpha * x_alpha
+           SG_LOG( SG_GENERAL, SG_ALERT, "sqrt(" << SG_EQ_RAD_SQUARE_M - x_alpha * x_alpha
                    << ")" );
        }
 #endif
@@ -99,23 +94,23 @@ void sgGeocToGeod( double lat_geoc, double radius, double
 #ifdef DOMAIN_ERR_DEBUG
        if ( errno ) {
            perror("fgGeocToGeod()");
-           FG_LOG( FG_GENERAL, FG_ALERT, "sqrt(" <<
+           SG_LOG( SG_GENERAL, SG_ALERT, "sqrt(" <<
                    1-EPS*EPS*sin_mu_a*sin_mu_a << ")" );
        }
 #endif
-       rho_alpha = EQUATORIAL_RADIUS_M*(1-EPS)/
+       rho_alpha = SG_EQUATORIAL_RADIUS_M*(1-EPS)/
            (denom*denom*denom);
        delt_mu = atan2(l_point*sin(delt_lambda),rho_alpha + *alt);
        *lat_geod = mu_alpha - delt_mu;
        lambda_sl = atan( E*E * tan(*lat_geod) ); // SL geoc. latitude
        sin_lambda_sl = sin( lambda_sl );
        *sea_level_r = 
-           sqrt(RESQ_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl));
+           sqrt(SG_EQ_RAD_SQUARE_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl));
 #ifdef DOMAIN_ERR_DEBUG
        if ( errno ) {
            perror("fgGeocToGeod()");
-           FG_LOG( FG_GENERAL, FG_ALERT, "sqrt(" <<
-                   RESQ_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl)
+           SG_LOG( SG_GENERAL, SG_ALERT, "sqrt(" <<
+                   SG_EQ_RAD_SQUARE_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl)
                    << ")" );
        }
 #endif
@@ -137,23 +132,27 @@ void sgGeocToGeod( double lat_geoc, double radius, double
 //
 
 
-void sgGeodToGeoc( double lat_geod, double alt, double *sl_radius,
+void sgGeodToGeoc( const double& lat_geod, const double& alt, double *sl_radius,
                      double *lat_geoc )
 {
     double lambda_sl, sin_lambda_sl, cos_lambda_sl, sin_mu, cos_mu, px, py;
     
+#ifdef DOMAIN_ERR_DEBUG
+    errno = 0;
+#endif
+
     lambda_sl = atan( E*E * tan(lat_geod) ); // sea level geocentric latitude
     sin_lambda_sl = sin( lambda_sl );
     cos_lambda_sl = cos( lambda_sl );
     sin_mu = sin(lat_geod);                  // Geodetic (map makers') latitude
     cos_mu = cos(lat_geod);
     *sl_radius = 
-       sqrt(RESQ_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl));
+       sqrt(SG_EQ_RAD_SQUARE_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl));
 #ifdef DOMAIN_ERR_DEBUG
        if ( errno ) {
            perror("fgGeodToGeoc()");
-           FG_LOG( FG_GENERAL, FG_ALERT, "sqrt(" <<
-                   RESQ_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl)
+           SG_LOG( SG_GENERAL, SG_ALERT, "sqrt(" <<
+                   SG_EQ_RAD_SQUARE_M / (1 + ((1/(E*E))-1)*sin_lambda_sl*sin_lambda_sl)
                    << ")" );
        }
 #endif
@@ -176,14 +175,14 @@ void sgGeodToGeoc( double lat_geod, double alt, double *sl_radius,
 //
 // modified for FlightGear to use WGS84 only -- Norman Vine
 
-#define GEOD_INV_PI FG_PI
+#define GEOD_INV_PI SGD_PI
 
 // s == distance
 // az = azimuth
 
 // for WGS_84 a = 6378137.000, rf = 298.257223563;
 
-static double M0( double e2 ) {
+static inline double M0( double e2 ) {
     //double e4 = e2*e2;
     return GEOD_INV_PI*(1.0 - e2*( 1.0/4.0 + e2*( 3.0/64.0 + 
                                                  e2*(5.0/256.0) )))/2.0;
@@ -192,8 +191,10 @@ static double M0( double e2 ) {
 
 // given, alt, lat1, lon1, az1 and distance (s), calculate lat2, lon2
 // and az2.  Lat, lon, and azimuth are in degrees.  distance in meters
-int geo_direct_wgs_84 ( double alt, double lat1, double lon1, double az1, 
-                       double s, double *lat2, double *lon2,  double *az2 )
+int geo_direct_wgs_84 ( const double& alt, const double& lat1,
+                        const double& lon1, const double& az1,
+                        const double& s, double *lat2, double *lon2,
+                        double *az2 )
 {
     double a = 6378137.000, rf = 298.257223563;
     double RADDEG = (GEOD_INV_PI)/180.0, testv = 1.0E-10;
@@ -274,7 +275,8 @@ int geo_direct_wgs_84 ( double alt, double lat1, double lon1, double az1,
     } else {                   // phi1 == 90 degrees, polar origin
        double dM = a*M0(e2) - s;
        double paz = ( phi1 < 0.0 ? 180.0 : 0.0 );
-       return geo_direct_wgs_84( alt, 0.0, lon1, paz, dM,lat2,lon2,az2 );
+        double zero = 0.0f;
+       return geo_direct_wgs_84( alt, zero, lon1, paz, dM, lat2, lon2, az2 );
     } 
 }
 
@@ -282,8 +284,10 @@ int geo_direct_wgs_84 ( double alt, double lat1, double lon1, double az1,
 // given alt, lat1, lon1, lat2, lon2, calculate starting and ending
 // az1, az2 and distance (s).  Lat, lon, and azimuth are in degrees.
 // distance in meters
-int geo_inverse_wgs_84( double alt, double lat1, double lon1, double lat2,
-                       double lon2, double *az1, double *az2, double *s )
+int geo_inverse_wgs_84( const double& alt, const double& lat1,
+                        const double& lon1, const double& lat2,
+                       const double& lon2, double *az1, double *az2,
+                        double *s )
 {
     double a = 6378137.000, rf = 298.257223563;
     int iter=0;
@@ -310,8 +314,9 @@ int geo_inverse_wgs_84( double alt, double lat1, double lon1, double lat2,
        return 0;
     } else if( fabs(cosphi2) < testv ) {
        // terminal point is polar
-       int k = geo_inverse_wgs_84( alt, lat1,lon1,lat1,lon1+180.0, 
-                                   az1,az2,s );
+        double _lon1 = lon1 + 180.0f;
+       int k = geo_inverse_wgs_84( alt, lat1, lon1, lat1, _lon1, 
+                                   az1, az2, s );
        k = k; // avoid compiler error since return result is unused
        *s /= 2.0;
        *az2 = *az1 + 180.0;