00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include <math.h>
00038 #include "gnss_error.h"
00039 #include "geodesy.h"
00040 #include "constants.h"
00041
00042
00043
00044
00045
00046
00047
00048 static BOOL GEODESY_GetReferenceEllipseParameters_A_E2(
00049 const GEODESY_enumReferenceEllipse ellipse,
00050 double* a,
00051 double* e2
00052 );
00053
00054
00055 static BOOL GEODESY_GetReferenceEllipseParameters_A_B_E2(
00056 const GEODESY_enumReferenceEllipse ellipse,
00057 double* a,
00058 double* b,
00059 double* e2
00060 );
00061
00062
00063 static BOOL GEODESY_IsLatitudeValid(
00064 const double latitude
00065 );
00066
00067
00068
00069
00070
00071 BOOL GEODESY_GetReferenceEllipseParameters(
00072 const GEODESY_enumReferenceEllipse ellipse,
00073 double* a,
00074 double* b,
00075 double* f_inv,
00076 double* e2
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
00203 BOOL GEODESY_GetReferenceEllipseParameters_A_E2(
00204 const GEODESY_enumReferenceEllipse ellipse,
00205 double* a,
00206 double* e2
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
00225 BOOL GEODESY_GetReferenceEllipseParameters_A_B_E2(
00226 const GEODESY_enumReferenceEllipse ellipse,
00227 double* a,
00228 double* b,
00229 double* e2
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
00246 BOOL GEODESY_IsLatitudeValid(
00247 const double latitude
00248 )
00249 {
00250
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,
00270 const double latitude,
00271 const double longitude,
00272 const double height,
00273 double *x,
00274 double *y,
00275 double *z
00276 )
00277 {
00278 double a;
00279 double e2;
00280 double N;
00281 double sinlat;
00282 double dtmp;
00283 BOOL result;
00284
00285
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
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,
00323 const double x,
00324 const double y,
00325 const double z,
00326 double *latitude,
00327 double *longitude,
00328 double *height
00329 )
00330 {
00331 double a;
00332 double b;
00333 double e2;
00334 double N;
00335 double p;
00336 double dtmp;
00337 double sinlat;
00338 double lat;
00339 double lon;
00340 double hgt;
00341 BOOL result;
00342
00343
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
00357
00358
00359
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
00378
00379
00380 lon = 2.0 * atan2( y , ( x + p ) );
00381
00382
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 );
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,
00405 const double referenceLatitude,
00406 const double referenceLongitude,
00407 const double referenceHeight,
00408 const double latitude,
00409 const double longitude,
00410 const double height,
00411 double *northing,
00412 double *easting,
00413 double *vertical
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;
00426 double B;
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
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
00488
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,
00502 const double referenceLatitude,
00503 const double referenceLongitude,
00504 const double referenceHeight,
00505 const double latitude,
00506 const double longitude,
00507 const double height,
00508 double *difference_northing,
00509 double *difference_easting,
00510 double *difference_vertical
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,
00533 const double latitude,
00534 double* M
00535 )
00536 {
00537 double a;
00538 double e2;
00539 double dtmp;
00540 BOOL result;
00541
00542
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 );
00553 dtmp = dtmp*dtmp*dtmp;
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,
00565 const double latitude,
00566 double* N
00567 )
00568 {
00569 double a;
00570 double e2;
00571 double W;
00572 BOOL result;
00573
00574
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,
00593 const double referenceLatitude,
00594 const double latitude,
00595 double* arc
00596 )
00597 {
00598 double a;
00599 double e2;
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;
00610 double arc_p;
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
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 );
00636 B = -dtmp * ( 0.375 * e2 + 0.46875 * e4 + 0.5126953125 * e6 + 0.538330078125 * e8 );
00637 C = dtmp * ( 0.05859375 * e4 + 0.1025390625 * e6 + 0.13458251953125 * e8 );
00638 D = -dtmp * ( 0.011393229167 * e6 + 0.025634765625 * 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,
00653 const double referenceLatitude,
00654 const double referenceLongitude,
00655 const double longitude,
00656 double* arc
00657 )
00658 {
00659 double a;
00660 double e2;
00661 double N;
00662 BOOL result;
00663
00664 *arc = 0;
00665
00666
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,
00693 const double referenceLongitude,
00694 const double dN,
00695 const double dE,
00696 const double dUp,
00697 double* dX,
00698 double* dY,
00699 double* dZ
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,
00734 const double referenceLongitude,
00735 const double dX,
00736 const double dY,
00737 const double dZ,
00738 double* dN,
00739 double* dE,
00740 double* dUp
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,
00775 const double fromX,
00776 const double fromY,
00777 const double fromZ,
00778 const double toX,
00779 const double toY,
00780 const double toZ,
00781 double* elevation,
00782 double* azimuth
00783 )
00784 {
00785 double lat;
00786 double lon;
00787 double dX;
00788 double dY;
00789 double dZ;
00790 double dN;
00791 double dE;
00792 double dUp;
00793 double tmp;
00794 BOOL result;
00795
00796 *elevation = 0;
00797 *azimuth = 0;
00798
00799
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
00815 dX = toX - fromX;
00816 dY = toY - fromY;
00817 dZ = toZ - fromZ;
00818
00819
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
00836 tmp = sqrt( dN*dN + dE*dE );
00837 *elevation = atan( dUp / tmp );
00838
00839
00840 *azimuth = atan2(dE, dN);
00841
00842
00843 if( *azimuth < 0.0 )
00844 *azimuth += TWOPI;
00845
00846 return TRUE;
00847 }
00848