geodesy.c

Go to the documentation of this file.
00001 /**
00002 \file    geodesy.c
00003 \brief   GNSS core 'c' function library: geodesy related functions.
00004 \author  Glenn D. MacGougan (GDM)
00005 \date    2007-11-28
00006 \since   2005-08-14
00007 
00008 \b "LICENSE INFORMATION" \n
00009 Copyright (c) 2007, refer to 'author' doxygen tags \n
00010 All rights reserved. \n
00011 
00012 Redistribution and use in source and binary forms, with or without
00013 modification, are permitted provided the following conditions are met: \n
00014 
00015 - Redistributions of source code must retain the above copyright
00016  notice, this list of conditions and the following disclaimer. \n
00017 - Redistributions in binary form must reproduce the above copyright
00018  notice, this list of conditions and the following disclaimer in the
00019  documentation and/or other materials provided with the distribution. \n
00020 - The name(s) of the contributor(s) may not be used to endorse or promote 
00021  products derived from this software without specific prior written 
00022  permission. \n
00023 
00024 THIS SOFTWARE IS PROVIDED BY THE CONTRIBUTORS ``AS IS'' AND ANY EXPRESS 
00025 OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED 
00026 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
00027 DISCLAIMED. IN NO EVENT SHALL THE CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
00028 INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
00029 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 
00030 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 
00031 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 
00032 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 
00033 OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 
00034 SUCH DAMAGE.
00035 */
00036 
00037 #include <math.h>
00038 #include "gnss_error.h"
00039 #include "geodesy.h"
00040 #include "constants.h"
00041 
00042 
00043 /*************************************************************************************************/
00044 // static function definitions (for functions used in this file only)
00045 // static functions are functions that are only visable to other functions in the same file.
00046 
00047 // Many functions contained herein only need a, and e2.
00048 static BOOL GEODESY_GetReferenceEllipseParameters_A_E2( 
00049   const GEODESY_enumReferenceEllipse ellipse,
00050   double* a,      // semi-major axis of the reference ellipse                     [m]
00051   double* e2      // eccentricity of the reference ellipse (e2 = (a*a-b*b)/(a*a)) [] 
00052   );
00053 
00054 // Many functions contained herein only need a, b and e2.
00055 static BOOL GEODESY_GetReferenceEllipseParameters_A_B_E2( 
00056   const GEODESY_enumReferenceEllipse ellipse,
00057   double* a,      // semi-major axis of the reference ellipse                     [m]
00058   double* b,      // semi-minor axis of the reference ellipse                     [m]
00059   double* e2      // eccentricity of the reference ellipse (e2 = (a*a-b*b)/(a*a)) [] 
00060   );
00061 
00062 // return TRUE(1) if valid, FALSE(0) otherwise.
00063 static BOOL GEODESY_IsLatitudeValid( 
00064   const double latitude //!< expecting a value -pi/2 to pi/2 [rad]
00065   );
00066 /*************************************************************************************************/
00067 
00068 
00069 
00070 
00071 BOOL GEODESY_GetReferenceEllipseParameters( 
00072   const GEODESY_enumReferenceEllipse ellipse, //!< reference ellipse enumerated []
00073   double* a,      //!< semi-major axis of the reference ellipse                     [m]
00074   double* b,      //!< semi-minor axis of the reference ellipse (b = a - a*f_inv)   [m] 
00075   double* f_inv,  //!< inverse of the flattening of the reference ellipse           []
00076   double* e2      //!< eccentricity of the reference ellipse (e2 = (a*a-b*b)/(a*a)) [] 
00077   )
00078 {
00079   switch( ellipse )
00080   {
00081   case GEODESY_REFERENCE_ELLIPSE_AIRY_1830: 
00082     *a     = GEODESY_REFERENCE_ELLIPSE_AIRY_1830_A;
00083     *f_inv = GEODESY_REFERENCE_ELLIPSE_AIRY_1830_F_INV;
00084     *b     = GEODESY_REFERENCE_ELLIPSE_AIRY_1830_B;
00085     *e2    = GEODESY_REFERENCE_ELLIPSE_AIRY_1830_E2;
00086     break;
00087     
00088   case GEODESY_REFERENCE_ELLIPSE_MODIFED_AIRY:
00089     *a     = GEODESY_REFERENCE_ELLIPSE_MODIFED_AIRY_A;
00090     *f_inv = GEODESY_REFERENCE_ELLIPSE_MODIFED_AIRY_F_INV;
00091     *b     = GEODESY_REFERENCE_ELLIPSE_MODIFED_AIRY_B;
00092     *e2    = GEODESY_REFERENCE_ELLIPSE_MODIFED_AIRY_E2;
00093     break;
00094     
00095   case GEODESY_REFERENCE_ELLIPSE_AUSTRALIAN_NATIONAL:
00096     *a     = GEODESY_REFERENCE_ELLIPSE_AUSTRALIAN_NATIONAL_A;
00097     *f_inv = GEODESY_REFERENCE_ELLIPSE_AUSTRALIAN_NATIONAL_F_INV;
00098     *b     = GEODESY_REFERENCE_ELLIPSE_AUSTRALIAN_NATIONAL_B;
00099     *e2    = GEODESY_REFERENCE_ELLIPSE_AUSTRALIAN_NATIONAL_E2;
00100     break;
00101     
00102   case GEODESY_REFERENCE_ELLIPSE_BESSEL_1841:
00103     *a     = GEODESY_REFERENCE_ELLIPSE_BESSEL_1841_A;
00104     *f_inv = GEODESY_REFERENCE_ELLIPSE_BESSEL_1841_F_INV;
00105     *b     = GEODESY_REFERENCE_ELLIPSE_BESSEL_1841_B;
00106     *e2    = GEODESY_REFERENCE_ELLIPSE_BESSEL_1841_E2;
00107     break;
00108     
00109   case GEODESY_REFERENCE_ELLIPSE_CLARKE_1866:
00110     *a     = GEODESY_REFERENCE_ELLIPSE_CLARKE_1866_A;
00111     *f_inv = GEODESY_REFERENCE_ELLIPSE_CLARKE_1866_F_INV;
00112     *b     = GEODESY_REFERENCE_ELLIPSE_CLARKE_1866_B;
00113     *e2    = GEODESY_REFERENCE_ELLIPSE_CLARKE_1866_E2;
00114     break;
00115     
00116   case GEODESY_REFERENCE_ELLIPSE_CLARKE_1880:
00117     *a     = GEODESY_REFERENCE_ELLIPSE_CLARKE_1880_A;
00118     *f_inv = GEODESY_REFERENCE_ELLIPSE_CLARKE_1880_F_INV;
00119     *b     = GEODESY_REFERENCE_ELLIPSE_CLARKE_1880_B;
00120     *e2    = GEODESY_REFERENCE_ELLIPSE_CLARKE_1880_E2;
00121     break;
00122     
00123   case GEODESY_REFERENCE_ELLIPSE_EVEREST_INDIA_1830:
00124     *a     = GEODESY_REFERENCE_ELLIPSE_EVEREST_INDIA_1830_A;
00125     *f_inv = GEODESY_REFERENCE_ELLIPSE_EVEREST_INDIA_1830_F_INV;
00126     *b     = GEODESY_REFERENCE_ELLIPSE_EVEREST_INDIA_1830_B;
00127     *e2    = GEODESY_REFERENCE_ELLIPSE_EVEREST_INDIA_1830_E2;
00128     break;
00129     
00130   case GEODESY_REFERENCE_ELLIPSE_EVEREST_BRUNEI_E_MALAYSIA:
00131     *a     = GEODESY_REFERENCE_ELLIPSE_EVEREST_BRUNEI_E_MALAYSIA_A;
00132     *f_inv = GEODESY_REFERENCE_ELLIPSE_EVEREST_BRUNEI_E_MALAYSIA_F_INV;
00133     *b     = GEODESY_REFERENCE_ELLIPSE_EVEREST_BRUNEI_E_MALAYSIA_B;
00134     *e2    = GEODESY_REFERENCE_ELLIPSE_EVEREST_BRUNEI_E_MALAYSIA_E2;
00135     break;
00136     
00137   case GEODESY_REFERENCE_ELLIPSE_EVEREST_W_MALAYSIA_SINGAPORE:
00138     *a     = GEODESY_REFERENCE_ELLIPSE_EVEREST_W_MALAYSIA_SINGAPORE_A;
00139     *f_inv = GEODESY_REFERENCE_ELLIPSE_EVEREST_W_MALAYSIA_SINGAPORE_F_INV;
00140     *b     = GEODESY_REFERENCE_ELLIPSE_EVEREST_W_MALAYSIA_SINGAPORE_B;
00141     *e2    = GEODESY_REFERENCE_ELLIPSE_EVEREST_W_MALAYSIA_SINGAPORE_E2;
00142     break;
00143     
00144   case GEODESY_REFERENCE_ELLIPSE_GRS_1980:
00145     *a     = GEODESY_REFERENCE_ELLIPSE_GRS_1980_A;
00146     *f_inv = GEODESY_REFERENCE_ELLIPSE_GRS_1980_F_INV;
00147     *b     = GEODESY_REFERENCE_ELLIPSE_GRS_1980_B;
00148     *e2    = GEODESY_REFERENCE_ELLIPSE_GRS_1980_E2;
00149     break;
00150     
00151   case GEODESY_REFERENCE_ELLIPSE_HELMERT_1906:
00152     *a     = GEODESY_REFERENCE_ELLIPSE_HELMERT_1906_A;
00153     *f_inv = GEODESY_REFERENCE_ELLIPSE_HELMERT_1906_F_INV;
00154     *b     = GEODESY_REFERENCE_ELLIPSE_HELMERT_1906_B;
00155     *e2    = GEODESY_REFERENCE_ELLIPSE_HELMERT_1906_E2;
00156     break;
00157     
00158   case GEODESY_REFERENCE_ELLIPSE_HOUGH_1960:
00159     *a     = GEODESY_REFERENCE_ELLIPSE_HOUGH_1960_A;
00160     *f_inv = GEODESY_REFERENCE_ELLIPSE_HOUGH_1960_F_INV;
00161     *b     = GEODESY_REFERENCE_ELLIPSE_HOUGH_1960_B;
00162     *e2    = GEODESY_REFERENCE_ELLIPSE_HOUGH_1960_E2;
00163     break;
00164     
00165   case GEODESY_REFERENCE_ELLIPSE_INTERNATIONAL_1924:
00166     *a     = GEODESY_REFERENCE_ELLIPSE_INTERNATIONAL_1924_A;
00167     *f_inv = GEODESY_REFERENCE_ELLIPSE_INTERNATIONAL_1924_F_INV;
00168     *b     = GEODESY_REFERENCE_ELLIPSE_INTERNATIONAL_1924_B;
00169     *e2    = GEODESY_REFERENCE_ELLIPSE_INTERNATIONAL_1924_E2;
00170     break;
00171     
00172   case GEODESY_REFERENCE_ELLIPSE_SOUTH_AMERICAN_1969:
00173     *a     = GEODESY_REFERENCE_ELLIPSE_SOUTH_AMERICAN_1969_A;
00174     *f_inv = GEODESY_REFERENCE_ELLIPSE_SOUTH_AMERICAN_1969_F_INV;
00175     *b     = GEODESY_REFERENCE_ELLIPSE_SOUTH_AMERICAN_1969_B;
00176     *e2    = GEODESY_REFERENCE_ELLIPSE_SOUTH_AMERICAN_1969_E2;
00177     break;
00178     
00179   case GEODESY_REFERENCE_ELLIPSE_WGS72:
00180     *a     = GEODESY_REFERENCE_ELLIPSE_WGS72_A;
00181     *f_inv = GEODESY_REFERENCE_ELLIPSE_WGS72_F_INV;
00182     *b     = GEODESY_REFERENCE_ELLIPSE_WGS72_B;
00183     *e2    = GEODESY_REFERENCE_ELLIPSE_WGS72_E2;
00184     break;
00185     
00186   case GEODESY_REFERENCE_ELLIPSE_WGS84:
00187     *a     = GEODESY_REFERENCE_ELLIPSE_WGS84_A;
00188     *f_inv = GEODESY_REFERENCE_ELLIPSE_WGS84_F_INV;
00189     *b     = GEODESY_REFERENCE_ELLIPSE_WGS84_B;
00190     *e2    = GEODESY_REFERENCE_ELLIPSE_WGS84_E2;
00191     break;
00192 
00193   default:
00194     GNSS_ERROR_MSG( "Unexpected default case." );
00195     return FALSE;
00196     break;  
00197   }
00198   return TRUE;
00199 }
00200 
00201 
00202 // static
00203 BOOL GEODESY_GetReferenceEllipseParameters_A_E2( 
00204   const GEODESY_enumReferenceEllipse ellipse,
00205   double* a,      // semi-major axis of the reference ellipse                     [m]
00206   double* e2      // eccentricity of the reference ellipse (e2 = (a*a-b*b)/(a*a)) [] 
00207   )
00208 {  
00209   double b;
00210   double f_inv;  
00211   BOOL result;
00212 
00213   result = GEODESY_GetReferenceEllipseParameters( 
00214     ellipse,
00215     a,
00216     &b,
00217     &f_inv,
00218     e2 );
00219 
00220   return result;
00221 }
00222 
00223 
00224 // static
00225 BOOL GEODESY_GetReferenceEllipseParameters_A_B_E2( 
00226   const GEODESY_enumReferenceEllipse ellipse,
00227   double* a,      // semi-major axis of the reference ellipse                     [m]
00228   double* b,      // semi-minor axis of the reference ellipse                     [m]
00229   double* e2      // eccentricity of the reference ellipse (e2 = (a*a-b*b)/(a*a)) [] 
00230   )
00231 {  
00232   double f_inv;  
00233   BOOL result;
00234 
00235   result = GEODESY_GetReferenceEllipseParameters( 
00236     ellipse,
00237     a,
00238     b,
00239     &f_inv,
00240     e2 );
00241 
00242   return result;
00243 }
00244 
00245 // static
00246 BOOL GEODESY_IsLatitudeValid( 
00247   const double latitude //!< expecting a value -pi/2 to pi/2 [rad]
00248   )
00249 {
00250   // check for valid latitude out of range
00251   if( latitude > HALFPI || latitude < -HALFPI )  
00252   {
00253     GNSS_ERROR_MSG( "if( latitude > HALFPI || latitude < -HALFPI )" );
00254     return FALSE;
00255   }
00256   else
00257   {
00258     return TRUE;
00259   }
00260 }
00261 
00262 
00263 
00264 
00265 
00266 
00267 
00268 BOOL GEODESY_ConvertGeodeticCurvilinearToEarthFixedCartesianCoordinates(
00269   const GEODESY_enumReferenceEllipse  referenceEllipse,  //!< reference ellipse enumerated []
00270   const double latitude,   //!< geodetic latitude                [rad]
00271   const double longitude,  //!< geodetic longitude               [rad]
00272   const double height,     //!< geodetic height                  [m]
00273   double *x,               //!< earth fixed cartesian coordinate [m]
00274   double *y,               //!< earth fixed cartesian coordinate [m]
00275   double *z                //!< earth fixed cartesian coordinate [m]
00276   )
00277 {  
00278   double a;      // semi-major axis of reference ellipse [m]
00279   double e2;     // first eccentricity of reference ellipse []
00280   double N;      // prime vertical radius of curvature [m]
00281   double sinlat; // sin of the latitude
00282   double dtmp;   // temp
00283   BOOL result;
00284 
00285   // get necessary reference ellipse parameters
00286   result = GEODESY_GetReferenceEllipseParameters_A_E2( referenceEllipse, &a, &e2 );
00287   if( result == FALSE )
00288   {
00289     *x = 0.0;
00290     *y = 0.0;
00291     *z = 0.0;
00292     GNSS_ERROR_MSG( "Reference ellipse invalid." );
00293     return FALSE;
00294   }
00295   
00296   // check for valid latitude out of range
00297   result = GEODESY_IsLatitudeValid( latitude );
00298   if( result == FALSE )
00299   {
00300     *x = 0.0;
00301     *y = 0.0;
00302     *z = 0.0;
00303     GNSS_ERROR_MSG( "Input latitude is invalid." );
00304     return FALSE;
00305   }
00306 
00307   sinlat = sin( latitude );             
00308   N = a / sqrt( 1.0 - e2 * sinlat*sinlat );      
00309   dtmp = (N + height) * cos(latitude);
00310 
00311   *x = dtmp * cos(longitude);
00312   *y = dtmp * sin(longitude);
00313   *z = ( (1.0 - e2)*N + height ) * sinlat;
00314 
00315   return TRUE;
00316 }
00317 
00318 
00319 
00320 
00321 BOOL GEODESY_ConvertEarthFixedCartesianToGeodeticCurvilinearCoordinates(
00322   const GEODESY_enumReferenceEllipse  referenceEllipse,  //!< reference ellipse enumerated []
00323   const double x,              // earth fixed cartesian coordinate [m]
00324   const double y,              // earth fixed cartesian coordinate [m]
00325   const double z,              // earth fixed cartesian coordinate [m]  
00326   double *latitude,            // geodetic latitude              [rad]
00327   double *longitude,           // geodetic longitude             [rad]
00328   double *height               // geodetic height                [m]
00329   )
00330 {
00331   double a;      // semi-major axis of reference ellipse [m]
00332   double b;      // semi-minor axis of reference ellipse [m]
00333   double e2;     // first eccentricity of reference ellipse []
00334   double N;      // prime vertical radius of curvature [m]
00335   double p;      // sqrt( x^2 + y^2 ) [m]
00336   double dtmp;   // temp
00337   double sinlat; // sin(lat)
00338   double lat;    // temp geodetic latitude  [rad]
00339   double lon;    // temp geodetic longitude [rad]
00340   double hgt;    // temp geodetic height    [m]
00341   BOOL result;
00342   
00343   // get necessary reference ellipse parameters
00344   result = GEODESY_GetReferenceEllipseParameters_A_B_E2( referenceEllipse, &a, &b, &e2 );
00345   if( result == FALSE )
00346   {
00347     *latitude  = 0;
00348     *longitude = 0;  
00349     *height    = 0;  
00350     GNSS_ERROR_MSG( "Reference ellipse invalid." );    
00351     return FALSE;
00352   }
00353   
00354   if( x == 0.0 && y == 0.0 ) 
00355   {
00356     // at a pole    
00357     // most likely to happen while using a simulator
00358     
00359     // longitude is really unknown
00360     lon = 0.0; 
00361     
00362     if( z < 0 )
00363     {
00364       hgt = -z - b;
00365       lat = -HALFPI;
00366     }
00367     else
00368     {
00369       hgt = z - b;
00370       lat = HALFPI;
00371     }
00372   }
00373   else
00374   {
00375     p = sqrt( x*x + y*y );
00376 
00377     // unique solution for longitude
00378     // best formula for any longitude and applies well near the poles
00379     // pp. 178 reference [2]
00380     lon = 2.0 * atan2( y , ( x + p ) );
00381     
00382     // set approximate initial latitude assuming a height of 0.0
00383     lat = atan( z / (p * (1.0 - e2)) );
00384     hgt = 0.0;
00385     do
00386     { 
00387       dtmp = hgt;
00388       sinlat = sin(lat);
00389       N   = a / sqrt( 1.0 - e2*sinlat*sinlat );
00390       hgt = p / cos(lat) - N;
00391       lat = atan( z / (p * ( 1.0 - e2*N/(N + hgt) )) );      
00392 
00393     } while( fabs( hgt - dtmp ) > 0.0001 );  // 0.1 mm convergence for height
00394   }
00395 
00396   *latitude  = lat;
00397   *longitude = lon;  
00398   *height    = hgt;
00399   return TRUE;
00400 }
00401 
00402 
00403 BOOL GEODESY_ComputeNorthingEastingVertical(
00404   const GEODESY_enumReferenceEllipse  referenceEllipse,  //!< reference ellipse enumerated []
00405   const double referenceLatitude,  //!< datum geodetic latitude  [rad]
00406   const double referenceLongitude, //!< datum geodetic longitude [rad]
00407   const double referenceHeight,    //!< datum geodetic height    [m]
00408   const double latitude,           //!< geodetic latitude        [rad]
00409   const double longitude,          //!< geodetic longitude       [rad]
00410   const double height,             //!< geodetic height          [m]
00411   double *northing,                //!< local geodetic northing  [m]
00412   double *easting,                 //!< local geodetic easting   [m]
00413   double *vertical                 //!< local geodetic vertical  [m]
00414   )
00415 {
00416   double x_ref;
00417   double y_ref;
00418   double z_ref;
00419   double x;
00420   double y;
00421   double z;
00422   double dx;
00423   double dy;
00424   double dz;
00425   double A;   // rotation angle [rad]
00426   double B;   // rotation angle [rad]
00427   double cosA;
00428   double sinA;
00429   double cosB;
00430   double sinB;
00431   BOOL result;
00432 
00433   *northing = 0;
00434   *easting  = 0;
00435   *vertical = 0;
00436 
00437   result = GEODESY_IsLatitudeValid( referenceLatitude );
00438   if( result == FALSE )
00439   {
00440     GNSS_ERROR_MSG( "Input reference latitude is invalid" );
00441     return FALSE;  
00442   }
00443   result = GEODESY_IsLatitudeValid( latitude );  
00444   if( result == FALSE )
00445   {
00446     GNSS_ERROR_MSG( "Input latitude is invalid." );
00447     return FALSE;  
00448   }
00449   
00450   result = GEODESY_ConvertGeodeticCurvilinearToEarthFixedCartesianCoordinates(
00451     referenceEllipse,
00452     referenceLatitude,
00453     referenceLongitude,
00454     referenceHeight,
00455     &x_ref,
00456     &y_ref,
00457     &z_ref );
00458   if( result == FALSE )
00459   {
00460     GNSS_ERROR_MSG( "GEODESY_ConvertGeodeticCurvilinearToEarthFixedCartesianCoordinates return FALSE." );
00461     return FALSE;
00462   }
00463   
00464   result = GEODESY_ConvertGeodeticCurvilinearToEarthFixedCartesianCoordinates(
00465     referenceEllipse,
00466     latitude,
00467     longitude,
00468     height,
00469     &x,
00470     &y,
00471     &z );
00472   if( result == FALSE )
00473   {
00474     GNSS_ERROR_MSG( "GEODESY_ConvertGeodeticCurvilinearToEarthFixedCartesianCoordinates return FALSE." );
00475     return FALSE;
00476   }
00477 
00478   // A and B are rotation angles
00479   A = referenceLatitude - HALFPI;
00480   B = referenceLongitude - PI;
00481 
00482   cosA = cos(A);
00483   sinA = sin(A);
00484   cosB = cos(B);
00485   sinB = sin(B);
00486 
00487   // the cartesian vector between the two points in the geodetic 
00488   // frame is rotated to the local geodetic frame
00489   dx = x - x_ref;
00490   dy = y - y_ref;
00491   dz = z - z_ref;   
00492 
00493   *northing = cosA*cosB * dx  +  cosA*sinB * dy  -  sinA*dz;
00494   *easting  = sinB      * dx  -  cosB      * dy;
00495   *vertical = sinA*cosB * dx  +  sinA*sinB * dy  +  cosA*dz;   
00496   return TRUE;
00497 }
00498 
00499 
00500 BOOL GEODESY_ComputePositionDifference(
00501   const GEODESY_enumReferenceEllipse  referenceEllipse,  //!< reference ellipse enumerated []
00502   const double referenceLatitude,  //!< reference point geodetic latitude  [rad]
00503   const double referenceLongitude, //!< reference point geodetic longitude [rad]
00504   const double referenceHeight,    //!< reference point geodetic height    [m]
00505   const double latitude,           //!< geodetic latitude        [rad]
00506   const double longitude,          //!< geodetic longitude       [rad]
00507   const double height,             //!< geodetic height          [m]
00508   double *difference_northing,     //!< difference in northing   [m]  (+2 m, means 2 m North of the reference)
00509   double *difference_easting,      //!< difference in easting    [m]  (+2 m, means 2 m East  of the reference)
00510   double *difference_vertical      //!< difference in vertical   [m]  (+2 m, means 2 m above of the reference)
00511   )
00512 {
00513   BOOL result;
00514   result = GEODESY_ComputeNorthingEastingVertical(
00515     referenceEllipse,
00516     referenceLatitude, 
00517     referenceLongitude,
00518     referenceHeight,   
00519     latitude,          
00520     longitude,         
00521     height,            
00522     difference_northing,
00523     difference_easting,
00524     difference_vertical );  
00525   return result;
00526 }
00527 
00528 
00529 
00530 
00531 BOOL GEODESY_ComputeMeridianRadiusOfCurvature(
00532   const GEODESY_enumReferenceEllipse  referenceEllipse, //!< reference ellipse enumerated []
00533   const double latitude,  //!< geodetic latitude                     [rad]
00534   double*  M              //!< computed meridian radius of curvature [m]
00535   )
00536 {
00537   double a;  // semi-major axis of reference ellipse    [m]
00538   double e2; // first eccentricity of reference ellipse []
00539   double dtmp;
00540   BOOL result;
00541   
00542   // get necessary reference ellipse parameters
00543   result = GEODESY_GetReferenceEllipseParameters_A_E2( referenceEllipse, &a, &e2 );
00544   if( result == FALSE )
00545   {
00546     *M = 0;
00547     GNSS_ERROR_MSG( "Reference ellipse invalid." );    
00548     return result;
00549   }
00550     
00551   dtmp = sin(latitude);
00552   dtmp = sqrt( 1.0 - e2 * dtmp * dtmp );  // W
00553   dtmp = dtmp*dtmp*dtmp;                  // W^3    
00554   
00555   *M = a * ( 1.0 - e2 ) / dtmp;
00556   return TRUE;
00557 }
00558 
00559 
00560 
00561 
00562 
00563 BOOL GEODESY_ComputePrimeVerticalRadiusOfCurvature(
00564   const GEODESY_enumReferenceEllipse  referenceEllipse, //!< reference ellipse enumerated []
00565   const double latitude,  //!< geodetic latitude                           [rad]
00566   double*  N              //!< computed prime vertical radius of curvature [m]
00567   )
00568 {
00569   double a;  // semi-major axis of reference ellipse [m]
00570   double e2; // first eccentricity of reference ellipse []
00571   double W;
00572   BOOL result;
00573   
00574   // get necessary reference ellipse parameters
00575   result = GEODESY_GetReferenceEllipseParameters_A_E2( referenceEllipse, &a, &e2 );
00576   if( result == FALSE )
00577   {
00578     *N = 0;
00579     GNSS_ERROR_MSG( "Reference ellipse invalid." );    
00580     return result;  
00581   }
00582   
00583   W = sin(latitude);
00584   W = sqrt( 1.0 - e2 * W * W );      
00585   
00586   *N = a / W;
00587   return TRUE;
00588 }
00589 
00590 
00591 BOOL GEODESY_ComputeMeridianArcBetweenTwoLatitudes(
00592   const GEODESY_enumReferenceEllipse  referenceEllipse, //!< reference ellipse enumerated []
00593   const double referenceLatitude,  //!< datum geodetic latitude  [rad]
00594   const double latitude,           //!< geodetic latitude        [rad]
00595   double*      arc                 //!< computed meridian arc, North +ve, South -ve [m]
00596   )
00597 {
00598   double a;  // semi-major axis of reference ellipse [m]
00599   double e2; // first eccentricity of reference ellipse []
00600   double e4; 
00601   double e6;
00602   double e8;
00603   double dtmp;
00604   double A;
00605   double B;
00606   double C;
00607   double D; 
00608   double E;
00609   double arc_ref; // arc from equator for the reference lat [m]
00610   double arc_p;   // arc from eqautor for point 'P' [m]
00611   BOOL result;
00612 
00613   *arc = 0;
00614 
00615   result = GEODESY_IsLatitudeValid( referenceLatitude );
00616   if( result == FALSE )
00617   {
00618     GNSS_ERROR_MSG( "Reference latitude is invalid." );
00619     return result;
00620   }
00621   
00622   // get necessary reference ellipse parameters
00623   result = GEODESY_GetReferenceEllipseParameters_A_E2( referenceEllipse, &a, &e2 );
00624   if( result == FALSE )
00625   {
00626     GNSS_ERROR_MSG( "Reference ellipse invalid." );    
00627     return result;  
00628   }
00629     
00630   e4 = e2*e2;
00631   e6 = e4*e2;
00632   e8 = e6*e2;
00633   dtmp = a*(1.0-e2);
00634 
00635   A =  dtmp * ( 1.0 + 0.75  * e2 + 0.703125   * e4 + 0.68359375     * e6 + 0.67291259765625 * e8 ); //  dtmp * (1.0 + 3.0/4.0*e2 + 45.0/64.0*e4  + 175.0/256.0*e6  + 11025.0/16384.0*e8 );
00636   B = -dtmp * (       0.375 * e2 + 0.46875    * e4 + 0.5126953125   * e6 + 0.538330078125   * e8 ); // -dtmp * (      3.0/8.0*e2 + 15.0/32.0*e4  + 525.0/1024.0*e6 + 2205.0/4096.0*e8 );
00637   C =  dtmp * (                    0.05859375 * e4 + 0.1025390625   * e6 + 0.13458251953125 * e8 ); // -dtmp * (                   15.0/256.0*e4 + 105.0/1024.0*e6 + 2205.0/16384.0*e8 );  
00638   D = -dtmp * (                                      0.011393229167 * e6 + 0.025634765625   * e8 ); // -dtmp * (                                   35.0/3072.0*e6  + 105.0/4096.0*e8 );  
00639   E =  dtmp * ( 2.4032593e-03 * e8 );                                                             
00640 
00641   arc_ref = A*referenceLatitude + B*sin(2.0*referenceLatitude) + C*sin(4.0*referenceLatitude) + D*sin(6.0*referenceLatitude) + E*sin(8.0*referenceLatitude);
00642   arc_p = A*latitude + B*sin(2.0*latitude) + C*sin(4.0*latitude) + D*sin(6.0*latitude) + E*sin(8.0*latitude);
00643 
00644   *arc = arc_p - arc_ref;
00645   return TRUE;
00646 }
00647 
00648 
00649 
00650 
00651 BOOL GEODESY_ComputeParallelArcBetweenTwoLongitudes(
00652   const GEODESY_enumReferenceEllipse  referenceEllipse, //!< reference ellipse enumerated []
00653   const double referenceLatitude,  //!< reference geodetic latitude  [rad]
00654   const double referenceLongitude, //!< reference geodetic longitude [rad]
00655   const double longitude,          //!< geodetic longitude           [rad]
00656   double*      arc                 //!< computed parallel arc, East +ve, West -ve [m]
00657   )
00658 {
00659   double a;  // semi-major axis of reference ellipse [m]
00660   double e2; // first eccentricity of reference ellipse []
00661   double N;  // computed prime vertical radius of curvature [m]
00662   BOOL result;
00663 
00664   *arc = 0;
00665   
00666   // get necessary reference ellipse parameters
00667   result = GEODESY_GetReferenceEllipseParameters_A_E2( referenceEllipse, &a, &e2 );
00668   if( result == FALSE )
00669   {
00670     GNSS_ERROR_MSG( "Reference ellipse invalid." );    
00671     return result;    
00672   }
00673 
00674   result = GEODESY_IsLatitudeValid( referenceLatitude );
00675   if( result == FALSE )
00676   {
00677     GNSS_ERROR_MSG( "Reference latitude is invalid." );
00678     return result;    
00679   }
00680   
00681   N = sin(referenceLatitude);
00682   N = a / sqrt( 1.0 - e2 * N * N );
00683 
00684   *arc = N * cos(referenceLatitude) * (longitude - referenceLongitude);  
00685   return TRUE;
00686 }
00687 
00688 
00689 
00690 
00691 BOOL GEODESY_RotateVectorFromLocalGeodeticFrameToEarthFixedFrame(
00692   const double referenceLatitude,  //!< reference geodetic latitude                 [rad]
00693   const double referenceLongitude, //!< reference geodetic longitude                [rad]
00694   const double dN,                 //!< local geodetic northing vector component    [m]
00695   const double dE,                 //!< local geodetic easting  vector component    [m]
00696   const double dUp,                //!< local geodetic vertical vector component    [m]
00697   double* dX,                      //!< earth centered earth fixed vector component [m]
00698   double* dY,                      //!< earth centered earth fixed vector component [m]
00699   double* dZ                       //!< earth centered earth fixed vector component [m]
00700   )
00701 {
00702   double sinlat;
00703   double coslat;
00704   double sinlon;
00705   double coslon;
00706   BOOL result;
00707 
00708   result = GEODESY_IsLatitudeValid( referenceLatitude );
00709   if( result == FALSE )
00710   {
00711     *dX = 0;
00712     *dY = 0;
00713     *dZ = 0;
00714     GNSS_ERROR_MSG( "Reference latitude is invalid." );
00715     return result;
00716   }
00717 
00718   sinlat = sin(referenceLatitude);
00719   coslat = cos(referenceLatitude);
00720   sinlon = sin(referenceLongitude);
00721   coslon = cos(referenceLongitude);
00722   
00723   *dX = -sinlat*coslon * dN  -  sinlon * dE  +  coslat*coslon * dUp;
00724   *dY = -sinlat*sinlon * dN  +  coslon * dE  +  coslat*sinlon * dUp;
00725   *dZ =  coslat        * dN                  +  sinlat        * dUp;
00726   return TRUE;
00727 }
00728 
00729 
00730 
00731 
00732 BOOL GEODESY_RotateVectorFromEarthFixedFrameToLocalGeodeticFrame(
00733   const double referenceLatitude,  //!< reference geodetic latitude                 [rad]
00734   const double referenceLongitude, //!< reference geodetic longitude                [rad]
00735   const double dX,                 //!< earth centered earth fixed vector component [m]
00736   const double dY,                 //!< earth centered earth fixed vector component [m]
00737   const double dZ,                 //!< earth centered earth fixed vector component [m]
00738   double* dN,                      //!< local geodetic northing vector component    [m]
00739   double* dE,                      //!< local geodetic easting  vector component    [m]
00740   double* dUp                      //!< local geodetic vertical vector component    [m]
00741   )
00742 {
00743   double sinlat;
00744   double coslat;
00745   double sinlon;
00746   double coslon;
00747   BOOL result;
00748 
00749   result = GEODESY_IsLatitudeValid( referenceLatitude );
00750   if( result == FALSE )
00751   {
00752     *dN = 0;
00753     *dE = 0;
00754     *dUp = 0;
00755     GNSS_ERROR_MSG( "Reference latitude is invalid." );
00756     return result;    
00757   }
00758 
00759   sinlat = sin(referenceLatitude);
00760   coslat = cos(referenceLatitude);
00761   sinlon = sin(referenceLongitude);
00762   coslon = cos(referenceLongitude);
00763   
00764   *dN  = -sinlat*coslon * dX  -  sinlat*sinlon * dY  +  coslat * dZ;
00765   *dE  = -sinlon        * dX  +  coslon        * dY;
00766   *dUp =  coslat*coslon * dX  +  coslat*sinlon * dY  +  sinlat * dZ;  
00767 
00768   return TRUE;
00769 }
00770 
00771 
00772 
00773 BOOL GEODESY_ComputeAzimuthAndElevationAnglesBetweenToPointsInTheEarthFixedFrame(
00774   const GEODESY_enumReferenceEllipse  referenceEllipse, //!< reference ellipse enumerated []
00775   const double fromX, //!< earth centered earth fixed vector from point X component [m]
00776   const double fromY, //!< earth centered earth fixed vector from point Y component [m]
00777   const double fromZ, //!< earth centered earth fixed vector from point Z component [m]
00778   const double toX,   //!< earth centered earth fixed vector to point X component   [m]
00779   const double toY,   //!< earth centered earth fixed vector to point Y component   [m]
00780   const double toZ,   //!< earth centered earth fixed vector to point Z component   [m]
00781   double* elevation,  //!< elevation angle [rad]
00782   double* azimuth     //!< azimuth angle   [rad]
00783   )
00784 {
00785   double lat; // reference geodetic latitude  ('from' point) [rad]
00786   double lon; // reference geodetic longitude ('from' point) [rad]
00787   double dX;  // ECEF X vector component between 'from' and 'to' point (m)
00788   double dY;  // ECEF Y vector component between 'from' and 'to' point (m)
00789   double dZ;  // ECEF Z vector component between 'from' and 'to' point (m)
00790   double dN;  // LG northing vector component between 'from' and 'to' point (m)
00791   double dE;  // LG easting  vector component between 'from' and 'to' point (m)
00792   double dUp; // LG vertical vector component between 'from' and 'to' point (m)
00793   double tmp; // temp value
00794   BOOL result;
00795 
00796   *elevation = 0;
00797   *azimuth = 0; 
00798 
00799   // get the reference geodetic curvilinear coordinates from the 'from' point
00800   result = GEODESY_ConvertEarthFixedCartesianToGeodeticCurvilinearCoordinates(
00801     referenceEllipse,
00802     fromX,
00803     fromY,
00804     fromZ,
00805     &lat,
00806     &lon,
00807     &tmp );
00808   if( result == FALSE )
00809   {
00810     GNSS_ERROR_MSG( "GEODESY_ConvertEarthFixedCartesianToGeodeticCurvilinearCoordinates returned FALSE." );
00811     return result;      
00812   }
00813 
00814   // vector between the two points in the earth fixed frame
00815   dX = toX - fromX;
00816   dY = toY - fromY;
00817   dZ = toZ - fromZ;
00818 
00819   // rotate the vector to the local geodetic frame
00820   result = GEODESY_RotateVectorFromEarthFixedFrameToLocalGeodeticFrame(
00821     lat,
00822     lon,
00823     dX,
00824     dY,
00825     dZ,
00826     &dN,
00827     &dE,
00828     &dUp );
00829   if( result == FALSE )
00830   {
00831     GNSS_ERROR_MSG( "GEODESY_RotateVectorFromEarthFixedFrameToLocalGeodeticFrame returned FALSE." );
00832     return result;    
00833   }
00834   
00835   // compute the elevation
00836   tmp = sqrt( dN*dN + dE*dE );
00837   *elevation = atan( dUp / tmp );
00838   
00839   // compute the azimuth
00840   *azimuth = atan2(dE, dN);
00841 
00842   // by convention, azimuth will be between 0 to 2 PI
00843   if( *azimuth < 0.0 )
00844     *azimuth += TWOPI;
00845   
00846   return TRUE;
00847 }
00848 
SourceForge.net Logo