00001 00003 // // 00004 // UNCLASSIFIED UNCLASSIFIED UNCLASSIFIED UNCLASSIFIED // 00005 // // 00006 // Description of this module: // 00007 // // 00008 // Utility software that interpolates EGM 2008 // 00009 // geoid heights from one of NGA's geoid height grids. // 00010 // // 00011 // This base class supports two geoid height interpolators: // 00012 // // 00013 // 1) the first one reads the entire worldwide EGM2008 grid // 00014 // before interpolating any geoid heights from the worldwide grid; // 00015 // 2) the second one does not read geoid height data from // 00016 // file until a user requests that a geoid height be computed; // 00017 // the 2nd interpolator then reads an Area of Interest // 00018 // (AOI) grid from file before interpolating from the AOI grid; // 00019 // the Area of Interest grid is refreshed, if needed, // 00020 // when users submit subsequent geoid height requests. // 00021 // // 00022 // Revision History: // 00023 // Date Name Description // 00024 // ----------- ------------ ----------------------------------------------// 00025 // 19 Nov 2010 RD Craig Release // 00026 // 27 Jan 2011 S. Gillis BAEts28121, Terrain Service rearchitecture // 00027 // 17 May 2011 T. Thompson BAEts27393, let user know when problem // 00028 // is due to undefined MSPCCS_DATA // 00029 // // 00031 00032 // This file contains definitions 00033 // for functions in the Egm2008GeoidGrid class. 00034 00035 // THIS PACKAGE REQUIRES NGA's WORLDWIDE GEOID HEIGHT GRID. THIS GRID'S 00036 // DIRECTORY PATH MUST BE SPECIFIED BY ENVIRONMENT VARIABLE "MSPCCS_DATA". 00037 00038 // Note: Following common usage, this class uses the term 00039 // "geoid height" to mean heights of the geoid above or 00040 // below the earth ellipsoid. The GeoidLibrary class confuses 00041 // the term "geoid height" with height of a point above or below 00042 // the geoid: these heights are properly called "orthometric heights". 00043 00044 #include <fstream> 00045 00046 #ifdef IRIXN32 00047 #include <math.h> 00048 #else 00049 #include <cmath> 00050 #endif 00051 00052 #include "CoordinateConversionException.h" 00053 00054 #include "egm2008_geoid_grid.h" 00055 00056 namespace { 00057 00058 const double EPSILON = 1.0e-15; // ~4 times machine epsilon 00059 00060 const double PI = 3.14159265358979323; 00061 00062 const double PIDIV2 = PI / 2.0; 00063 const double TWOPI = 2.0 * PI; 00064 00065 const double DEG_PER_RAD = 180.0 / PI; 00066 const double RAD_PER_DEG = PI / 180.0; 00067 00068 const double SEMI_MAJOR_AXIS = 6378137.0; // WGS-84 00069 const double FLATTENING = 1.0 / 298.257223563; // WGS-84 00070 } 00071 00072 using namespace MSP; 00073 00074 // *************************** 00075 // ** Public user functions ** 00076 // *************************** 00077 00078 // ******************************** 00079 // * Egm2008GeoidGrid constructor * 00080 // ******************************** 00081 00082 Egm2008GeoidGrid::Egm2008GeoidGrid (void) 00083 { 00084 // November 19, 2010: Version 1.00 00085 00086 // This function implements the 00087 // default Egm2008GeoidGrid constructor. 00088 00089 //_mutex.lock(); // Use CCSThreadMutex function in constructors 00090 00091 int i; 00092 int length; 00093 char* pathName = NULL; 00094 00095 // Get reformatted geoid height grid's file name ..... 00096 pathName = getenv( "MSPCCS_DATA" ); 00097 00098 if ( NULL == pathName ) 00099 { 00100 strcpy( pathName, "../../data" ); 00101 } 00102 00103 length = strlen( pathName ); 00104 for (i = 0; i < length; i++) _gridFname += pathName[i]; 00105 00106 _gridFname += 00107 "/Und_min2.5x2.5_egm2008_WGS84_TideFree_reformatted"; 00108 00109 // BAEts27393 Reverse the change for this DR for plug-in backward 00110 // compatible with MSP 1.1, i.e. do not throw error when plug-in is 00111 // dropped into MSP 1.1 where EGM2008 data file is not available. 00112 00113 //int undefEnvVar = 0; 00114 //FILE* fp = NULL; /* File pointer to file */ 00115 // 00117 //char* fileName = new char[50]; 00118 //strcpy( fileName, "Und_min2.5x2.5_egm2008_WGS84_TideFree_reformatted" ); 00119 // 00121 //char* dataPath = getenv( "MSPCCS_DATA" ); 00122 //if ( NULL == dataPath ) 00123 //{ 00124 // undefEnvVar = 1; 00125 // strcpy( dataPath, "../../data" ); 00126 //} 00127 // 00129 //std::string pathName = dataPath; 00130 //pathName.append("/"); 00131 //pathName.append(fileName); 00132 // 00134 //if ( ( fp = fopen( pathName.c_str(), "r" ) ) == NULL ) 00135 //{ 00136 // ErrorControl errCtrl; 00137 // if ( 1 == undefEnvVar ) 00138 // { 00139 // errCtrl.issueError( "Environment variable undefined: MSPCCS_DATA.", 00140 // "Egm2008GeoidGrid::Egm2008GeoidGrid"); 00141 // } 00142 // else 00143 // { 00144 // std::string errMsg = "Unable to open "; 00145 // errMsg.append(pathName); 00146 // errCtrl.issueError( errMsg, "Egm2008GeoidGrid::Egm2008GeoidGrid" ); 00147 // } 00148 //} 00149 //fclose( fp ); 00150 00152 //_gridFname = pathName.c_str(); 00153 00154 00155 // Initialize base grid parameters ..... 00156 00157 MAX_WSIZE = 20; 00158 00159 _nGridPad = 0; 00160 _nGridCols = 0; 00161 _nGridRows = 0; 00162 _nOrigCols = 0; 00163 _nOrigRows = 0; 00164 00165 _baseLatitude = 0.0; 00166 _baseLongitude = 0.0; 00167 00168 _dLat = 0.0; 00169 _dLon = 0.0; 00170 00171 // The thread will be unlocked 00172 // by a derived class constructor. 00173 00174 } // End of Egm2008GeoidGrid constuctor 00175 00176 00177 // ************************************* 00178 // * Egm2008GeoidGrid copy constructor * 00179 // ************************************* 00180 00181 Egm2008GeoidGrid::Egm2008GeoidGrid 00182 ( const Egm2008GeoidGrid& oldGrid ) // input 00183 { 00184 // November 19, 2010: Version 1.00 00185 00186 // This function implements the 00187 // Egm2008GeoidGrid copy constructor. 00188 00189 // Definition: 00190 00191 // oldGrid: The Egm2008GeoidGrid object being copied. 00192 00193 oldGrid._mutex.lock(); // Use CCSThreadMutex function in copy constructors 00194 00195 // Copy worldwide grid's file name ..... 00196 00197 _gridFname = oldGrid._gridFname; 00198 00199 // Copy grid's parameters ..... 00200 00201 MAX_WSIZE = oldGrid.MAX_WSIZE; 00202 00203 _nGridPad = oldGrid._nGridPad; 00204 00205 _nGridRows = oldGrid._nGridRows; 00206 _nGridCols = oldGrid._nGridCols; 00207 00208 _nOrigRows = oldGrid._nOrigRows; 00209 _nOrigCols = oldGrid._nOrigCols; 00210 00211 _baseLatitude = oldGrid._baseLatitude; 00212 _baseLongitude = oldGrid._baseLongitude; 00213 00214 _dLat = oldGrid._dLat; 00215 _dLon = oldGrid._dLon; 00216 00217 // The oldGrid object will be unlocked 00218 // by a derived class copy constructor. 00219 00220 } // End of Egm2008GeoidGrid copy constuctor 00221 00222 00223 // ******************************* 00224 // * Egm2008GeoidGrid destructor * 00225 // ******************************* 00226 00227 Egm2008GeoidGrid::~Egm2008GeoidGrid (void) 00228 { 00229 // November 19, 2010: Version 1.00 00230 00231 // This function implements 00232 // the Egm2008GeoidGrid destructor. 00233 00234 ; 00235 00236 } // End of Egm2008GeoidGrid destructor 00237 00238 00239 // **************************************** 00240 // * Egm2008GeoidGrid assignment operator * 00241 // **************************************** 00242 00243 Egm2008GeoidGrid& 00244 Egm2008GeoidGrid::operator= ( const Egm2008GeoidGrid& oldGrid ) 00245 { 00246 // November 19, 2010: Version 1.00 00247 00248 // Definition: 00249 00250 // oldGrid: The Egm2008GeoidGrid object being assigned. 00251 00252 // This function implements the 00253 // Egm2008GeoidGrid assignment operator. 00254 00255 _mutex.lock(); // Use CCSThreadMutex function in assignments 00256 oldGrid._mutex.lock(); 00257 00258 if ( this == & oldGrid ) return ( *this ); 00259 00260 // Assign worldwide grid's file name ..... 00261 00262 _gridFname = oldGrid._gridFname; 00263 00264 // Assign grids' parameters ..... 00265 00266 MAX_WSIZE = oldGrid.MAX_WSIZE; 00267 00268 _nGridPad = oldGrid._nGridPad; 00269 00270 _nGridRows = oldGrid._nGridRows; 00271 _nGridCols = oldGrid._nGridCols; 00272 00273 _nOrigRows = oldGrid._nOrigRows; 00274 _nOrigCols = oldGrid._nOrigCols; 00275 00276 _baseLatitude = oldGrid._baseLatitude; 00277 _baseLongitude = oldGrid._baseLongitude; 00278 00279 _dLat = oldGrid._dLat; 00280 _dLon = oldGrid._dLon; 00281 00282 // The object will be unlocked by a 00283 // derived class assignment operator. 00284 00285 return( *this ); 00286 00287 } // End of Egm2008GeoidGrid assignment operator 00288 00289 00290 // ********************** 00291 // ** Hidden functions ** 00292 // ********************** 00293 00294 // ************************************** 00295 // * Retrieve user-specified grid point * 00296 // ************************************** 00297 00298 int 00299 Egm2008GeoidGrid::loadGridCoords( 00300 int i, // input 00301 int j, // input 00302 double& latitude, // output 00303 double& longitude ) // output 00304 { 00305 // November 19, 2010: Version 1.00 00306 00307 // This function retrieves geodetic 00308 // coordinates corresponding to geoid height posts. 00309 00310 // This function allows users to supply 00311 // indices refering to the worldwide grid's pad region. 00312 00313 // Definitions: 00314 00315 // i: Worldwide grid post's row index. 00316 // j: Worldwide grid post's column index. 00317 // latitude: Grid point's geodetic latitude (radians). 00318 // longitude: Grid point's geodetic longitude (radians). 00319 00320 // Thread locks are not needed here, because this function 00321 // can only be invoked from a bicubic spline geoidHeight function; the 00322 // thread locks, if any, reside in the bicubic spline's geoidHeight function. 00323 00324 try { 00325 00326 const int LIMIT1 = _nGridPad; 00327 const int LIMIT2 = _nGridRows - _nGridPad - 1; 00328 00329 int k; 00330 00331 // EDIT THE INPUT AND INITIALIZE ..... 00332 00333 if (( i < 0 ) || ( i >= _nGridRows )) return( 1 ); 00334 00335 while ( j < 0 ) j += _nGridCols; 00336 while ( j >= _nGridCols ) j -= _nGridCols; 00337 00338 // COMPUTE GEODETIC COORDINATES ..... 00339 00340 // Indices refer to Southern pad region ..... 00341 00342 if ( i < LIMIT1 ) 00343 { 00344 latitude = 00345 ( -PIDIV2 - _dLat * double( i - LIMIT1 )); 00346 00347 k = j + ( _nOrigCols / 2 ); 00348 00349 if ( k >= _nGridCols ) k -= _nOrigCols; 00350 00351 longitude = 00352 ( _baseLongitude + _dLon * double( k )); 00353 } 00354 00355 // Indices refer to NGA grid region ..... 00356 00357 if (( i >= LIMIT1 ) && ( i <= LIMIT2 )) 00358 { 00359 latitude = 00360 ( _baseLatitude + _dLat * double( i )); 00361 00362 longitude = 00363 ( _baseLongitude + _dLon * double( j )); 00364 } 00365 00366 // Indices refer to Northern pad region ..... 00367 00368 if ( i > LIMIT2 ) 00369 { 00370 latitude = 00371 ( PIDIV2 - _dLat * double( i - LIMIT2 )); 00372 00373 k = j + ( _nOrigCols / 2 ); 00374 00375 if ( k >= _nGridCols ) k -= _nOrigCols; 00376 00377 longitude = 00378 ( _baseLongitude + _dLon * double( k )); 00379 } 00380 00381 } // End of exceptions' try block 00382 00383 catch ( ... ) { return( 1 ); } 00384 00385 return( 0 ); // Normal-return flag 00386 00387 } // End of function Egm2008GeoidGrid::loadGridCoords 00388 00389 00390 // ********************************************** 00391 // * Compute a one-dimensional spline's moments * 00392 // ********************************************** 00393 00394 int 00395 Egm2008GeoidGrid::initSpline( 00396 int n, // input 00397 const double posts[], // input 00398 double moments[] ) // output 00399 { 00400 // November 4, 2010: Version 1.00 00401 // February 12, 2011: Version 2.00: 00402 // new variable names following code review 00403 00404 // Using equally-spaced abscissae, this 00405 // function computes a one-dimensional spline's moments. 00406 00407 // This function assumes that the endmost 00408 // support points' second derivatives are zero. 00409 00410 // Definitions: 00411 00412 // moments: Spline moments (2nd derivatives) 00413 // that are evaluated at the spline's support 00414 // points; this array contains at least n elements. 00415 // n: The spline's number of support points; the 00416 // calling function ensures that 3 <= n <= 20. 00417 // posts: An array containing the support points' 00418 // ordinates; this array contains at least n elements. 00419 00420 // w: An array that, on completion of the Gauss 00421 // elimination step, contains the triangularized 00422 // coefficient matrix' superdiagonal. Element w[0] 00423 // DOES NOT contain a superdiagonal matrix element. 00424 00425 // Thread locks are not needed here, because this function 00426 // can only be invoked from a bicubic spline geoidHeight function; 00427 // the thread locks reside in the bicubic spline's geoidHeight function. 00428 00429 // This C++ function 00430 // refines FORTRAN subroutine INITSP, 00431 // which was originally written (July 1983) 00432 // by Rene Forsberg, Institut Fur Erdmessung, 00433 // Universitat Hannover, Federal Republic of Germany. 00434 00435 // NGA used this FORTRAN function in 00436 // its internal EGM 2008 interpolation software. 00437 00438 // References: 00439 00440 // Einfuhrung in die Numerische Mathematik, Band I, 00441 // Josef Stoer, 00442 // Springer-Verlag, Heidelberg, 1972, pp 82 und 86 00443 00444 // Introduction to Numerical Analysis, 00445 // Josef Stoer and Roland Bulirsch, 00446 // Springer-Verlag, New York 1980, 00447 // Sections 2.4.1 - 2.4.2 (pp 93 - 102) 00448 00449 // Algorithm: 00450 00451 // For the special case of equally- 00452 // spaced abscissae, where the first and 00453 // last moments are zero, the moments are found by 00454 // solving the linear tridiagonal system of equations 00455 00456 // .. .. .. .. .. .. 00457 // : : : : : : 00458 // : 2.0 0.0 : : m : : 0 : 00459 // : : : 0 : : : 00460 // : 0.5 2.0 0.5 : : m : : rhs : 00461 // : : : 1 : : 1 : 00462 // : 0.5 2.0 0.5 : : m : : rhs : 00463 // : : : 2 : = : 2 :, 00464 // : . : : : : : 00465 // : . : : : : . : 00466 // : . : : : : . : 00467 // : 0.5 2.0 0.5 : : : : . : 00468 // : : : : : : 00469 // : 0.0 2.0 : : m : : 0 : 00470 // :. .: :. n-1.: :. .: 00471 00472 // where 00473 00474 // 1) n is the number of support data points, 00475 00476 // 2) m is the j-th support point's unknown second derivative, 00477 // j 00478 00479 // .. 00480 // : 0.0, j = 0, 00481 // : 00482 // : 00483 // 3) rhs = : 3 * (y - 2 y + y ), 1 <= j <= (n-2), 00484 // j : j+1 j j-1 00485 // : 00486 // : 0.0, j = n - 1. 00487 // :. 00488 00489 // This system of equations is never singular, 00490 // so tests for singularity are not necessary. 00491 00492 // The coefficient matrix is tridiagonal, 00493 // so system solution time is proportional 00494 // to the number of moments being computed. 00495 00496 try { 00497 00498 const int TWENTY = 20; 00499 00500 int k; 00501 00502 double p; 00503 00504 double w[ TWENTY ]; 00505 00506 // EDIT THE INPUT AND INITIALIZE ..... 00507 00508 if ( TWENTY != MAX_WSIZE ) return( 1 ); 00509 00510 if ( n < 3 ) return( 1 ); 00511 00512 w[ 0 ] = 0.0; 00513 moments[ 0 ] = 0.0; 00514 00515 // APPLY GAUSS ELIMINATION TO THE 00516 // MOMENTS' TRI-DIAGONAL SYSTEM OF EQUATIONS ..... 00517 00518 for (k = 2; k < n; k++) 00519 { 00520 p = ( w[ k-2 ] / 2.0 ) + 2.0; 00521 w[ k-1 ] = -0.5 / p; 00522 moments[ k-1 ] = 00523 (3.0 * 00524 ( posts[ k ] - ( 2.0 * posts[ k-1 ]) + 00525 posts[ k-2 ]) - ( moments[ k-2 ] / 2.0 )) / p; 00526 } 00527 00528 // SOLVE THE RESULTING BI-DIAGONAL SYSTEM ..... 00529 00530 moments[ n-1 ] = 0.0; 00531 00532 for (k = n - 1; k >= 2; k--) 00533 { 00534 moments[ k-1 ] += ( w[ k-1 ] * moments[ k ]); 00535 } 00536 00537 } // End of exceptions' try block 00538 00539 catch ( ... ) { return( 1 ); } 00540 00541 return( 0 ); // Normal-return flag 00542 00543 } // End of function Egm2008GeoidGrid::initSpline 00544 00545 00546 // *************************************************** 00547 // * Specialized one-dimensional spline interpolator * 00548 // *************************************************** 00549 00550 double 00551 Egm2008GeoidGrid::spline( 00552 int n, // input 00553 double x, // input 00554 const double posts[], // input 00555 const double moments[] ) // input 00556 { 00557 // November 4, 2010: Version 1.00 00558 // February 12, 2011: Version 2.00: 00559 // new variable names following code review 00560 00561 // Using support points having 00562 // equally-spaced abscissae, this function quickly 00563 // evaluates a cubic spline at an abscissa of interest. 00564 00565 // Definitions: 00566 00567 // moments: Spline moments (2nd derivatives) 00568 // that were evaluated at the spline's support 00569 // points; this array contains at least n elements. 00570 // n: The spline's number of support points; the 00571 // calling function ensures that 3 <= n <= 20. 00572 // posts: An array containing the support points' 00573 // ordinates; this array has at least n elements. 00574 // x: The interpolation argument; 00575 // (x == 0.0) ==> x equals the first support abscissa, 00576 // (x == double(n-1) ==> x equals the final support abscissa. 00577 00578 // Thread locks are not needed here, because this function 00579 // can only be invoked from a bicubic spline geoidHeight function; 00580 // the thread locks reside in the bicubic spline's geoidHeight function. 00581 00582 // This C++ function 00583 // refines FORTRAN function SPLINE, 00584 // which was originally written (June 1983) 00585 // by Rene Forsberg, Institut Fur Erdmessung, 00586 // Universitat Hannover, Federal Republic of Germany. 00587 00588 // NGA used this FORTRAN function in 00589 // its internal EGM 2008 interpolation software. 00590 00591 // References: 00592 00593 // Einfuhrung in die Numerische Mathematik, Band I, 00594 // Josef Stoer, 00595 // Springer-Verlag, Heidelberg, 1972, p 81 00596 00597 // Introduction to Numerical Analysis, 00598 // Josef Stoer and Roland Bulirsch, 00599 // Springer-Verlag, New York 1980, 00600 // Sections 2.4.1 and 2.4.2 (pp 93 - 102) 00601 00602 // When (x < 0) or (x > n-1), this function linearly 00603 // extrapolates the dependent variable from the nearest 00604 // support points. Recall that the endpoints' moments are zero. 00605 00606 double height; 00607 00608 const double MIN_ABSCISSA = 0.0; 00609 00610 int j; 00611 00612 double dx; 00613 double maxAbscissa; 00614 double mJ; 00615 double mJp1; 00616 00617 // Edit the input and initialize ..... 00618 00619 if ( n < 3 ) 00620 { 00621 throw MSP::CCS::CoordinateConversionException( 00622 "Error: Egm2008GeoidGrid::spline: not enough points"); 00623 } 00624 00625 try { 00626 00627 maxAbscissa = double( n - 1 ); 00628 00629 // Use one-dimensional extrapolation or 00630 // interpolation to determine geoid height ..... 00631 00632 if ( x <= MIN_ABSCISSA ) // linear extrapolation 00633 { 00634 height = 00635 posts[ 0 ] + 00636 ( x - MIN_ABSCISSA ) * 00637 ( posts[ 1 ] - posts[ 0 ] - 00638 ( moments[ 1 ] / 6.0 )); 00639 } 00640 else if ( x >= maxAbscissa ) // linear extrapolation 00641 { 00642 height = 00643 posts[ n-1 ] + 00644 ( x - maxAbscissa ) * 00645 ( posts[ n-1 ] - posts[ n-2 ] + 00646 ( moments[ n-2 ] / 6.0 )); 00647 } 00648 else // cubic spline interpolation 00649 { 00650 j = int( floor( x )); 00651 00652 dx = x - double( j ); 00653 mJ = moments[ j ]; 00654 mJp1 = moments[ j+1 ]; 00655 00656 height = // compute geoid height using Horner's rule 00657 posts[ j ] + 00658 dx * (( posts[ j+1 ] - posts[ j ] - 00659 ( mJ / 3.0 ) - ( mJp1 / 6.0 )) + 00660 dx * (( mJ / 2.0 ) + 00661 dx * ( mJp1 - mJ ) / 6.0 )); 00662 } 00663 00664 } // End of exceptions' try block 00665 00666 catch ( ... ) 00667 { 00668 throw MSP::CCS::CoordinateConversionException( 00669 "Error: Egm2008GeoidGrid::spline"); 00670 } 00671 00672 return ( height ); 00673 00674 } // End of function Egm2008GeoidGrid::spline 00675 00676 00677 // ************** 00678 // * Swap bytes * 00679 // ************** 00680 00681 void 00682 Egm2008GeoidGrid::swapBytes( 00683 void* buffer, // input & output 00684 size_t size, // input 00685 size_t count ) // input 00686 { 00687 // November 8, 2010: Version 1.00 00688 00689 // This function swaps 00690 // bytes when transforming between 00691 // BIG-ENDIAN and LITTLE-ENDIAN binary formats. 00692 00693 // Definitions: 00694 00695 // buffer: Buffer containing binary numbers. 00696 // count: Number of data items in the buffer. 00697 // size: Size (bytes) of each buffer data item. 00698 00699 // Thread locks are not needed here, because this function 00700 // can only be invoked from a bicubic spline geoidHeight function; 00701 // the thread locks reside in the bicubic spline's geoidHeight function. 00702 00703 char temp; 00704 00705 char* b = (char *) buffer; 00706 00707 for ( size_t c = 0; c < count; c ++ ) 00708 { 00709 for ( size_t s = 0; s < ( size / 2 ); s++ ) 00710 { 00711 temp = b[ c * size + s ]; 00712 b[ c * size + s] = b[ c * size + size - s - 1 ]; 00713 b[ c * size + size - s - 1] = temp; 00714 } 00715 } 00716 00717 } // End of function Egm2008GeoidGrid::swapBytes 00718 00719 00720 // ************************************* 00721 // * Find Southwest grid-point indices * 00722 // ************************************* 00723 00724 int 00725 Egm2008GeoidGrid::swGridIndices( 00726 double latitude, // input 00727 double longitude, // input 00728 int& i, // output 00729 int& j ) // output 00730 { 00731 // November 19, 2010: Version 1.00 00732 00733 // Given geodetic latitude and longitude, 00734 // this function finds indices to the nearest 00735 // grid point TO THE SOUTHWEST of the input point. 00736 // The indices refer to the padded version 00737 // of NGA's worldwide worldwide geoid-height grid. 00738 00739 // Definitions: 00740 00741 // i: Worldwide grid's latitude index. 00742 // j: Worldwide grid's longitude index. 00743 // latitude: Input geodetic latitude (radians). 00744 // longitude: Input geodetic longitude (radians). 00745 00746 // Thread locks are not needed here, because this function 00747 // can only be invoked from a bicubic spline geoidHeight function; 00748 // the thread locks reside in the bicubic spline's geoidHeight function. 00749 00750 // Edit the input and initialize ..... 00751 00752 try { 00753 00754 double maxLatitude; 00755 00756 maxLatitude = 00757 ( _baseLatitude + double( _nGridRows - 1 ) * _dLat ); 00758 00759 if (( latitude < _baseLatitude ) || 00760 ( latitude > maxLatitude )) return( 1 ); 00761 00762 while ( longitude < 0.0 ) longitude += TWOPI; 00763 while ( longitude > TWOPI ) longitude -= TWOPI; 00764 00765 // Compute the indices ..... 00766 00767 i = _nGridPad + int(( latitude + PIDIV2 ) / _dLat ); 00768 00769 j = _nGridPad + int( longitude / _dLon ); 00770 00771 } // End of exceptions' try block 00772 00773 catch ( ... ) { return( 1 ); } 00774 00775 return( 0 ); // Normal-return flag 00776 00777 } // End of function Egm2008GeoidGrid::swGridIndices 00778 00780 // UNCLASSIFIED UNCLASSIFIED UNCLASSIFIED UNCLASSIFIED // 00782