The Essential GNSS Project: test_geodesy.c Source File

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