115 using namespace MSP::CCS;
 
  122 #define EPSILON 1.75e-7    
  152 const double PI = 3.14159265358979323e0;  
 
  167 const double _6 = (6.0 * (
PI / 180.0));
 
  168 const double _8 = (8.0 * (
PI / 180.0));
 
  169 const double _72 = (72.0 * (
PI / 180.0));
 
  170 const double _80 = (80.0 * (
PI / 180.0));
 
  174 #define _500000  500000.0 
  182   double northing_offset; 
 
  186   {{
LETTER_C, 1100000.0, -72.0, -80.5, 0.0},
 
  187   {
LETTER_D, 2000000.0, -64.0, -72.0, 2000000.0},
 
  188   {
LETTER_E, 2800000.0, -56.0, -64.0, 2000000.0},
 
  189   {
LETTER_F, 3700000.0, -48.0, -56.0, 2000000.0},
 
  190   {
LETTER_G, 4600000.0, -40.0, -48.0, 4000000.0},
 
  191   {
LETTER_H, 5500000.0, -32.0, -40.0, 4000000.0},
 
  192   {
LETTER_J, 6400000.0, -24.0, -32.0, 6000000.0},
 
  193   {
LETTER_K, 7300000.0, -16.0, -24.0, 6000000.0},
 
  194   {
LETTER_L, 8200000.0, -8.0, -16.0, 8000000.0},
 
  195   {
LETTER_M, 9100000.0, 0.0, -8.0, 8000000.0},
 
  197   {
LETTER_P, 800000.0, 16.0, 8.0, 0.0},
 
  198   {
LETTER_Q, 1700000.0, 24.0, 16.0, 0.0},
 
  199   {
LETTER_R, 2600000.0, 32.0, 24.0, 2000000.0},
 
  200   {
LETTER_S, 3500000.0, 40.0, 32.0, 2000000.0},
 
  201   {
LETTER_T, 4400000.0, 48.0, 40.0, 4000000.0},
 
  202   {
LETTER_U, 5300000.0, 56.0, 48.0, 4000000.0},
 
  203   {
LETTER_V, 6200000.0, 64.0, 56.0, 6000000.0},
 
  204   {
LETTER_W, 7000000.0, 72.0, 64.0, 6000000.0},
 
  205   {
LETTER_X, 7900000.0, 84.5, 72.0, 6000000.0}};
 
  212   long ltr2_high_value;   
 
  213   long ltr3_high_value;   
 
  214   double false_easting;   
 
  215   double false_northing;  
 
  255   char alphabet[] = 
"ABCDEFGHIJKLMNOPQRSTUVWXYZ";
 
  259     i = sprintf (USNGString+i,
"%2.2ld",zone);
 
  261     strncpy(USNGString, 
"  ", 2);  
 
  264     USNGString[i++] = alphabet[letters[j]];
 
  265   divisor = pow (10.0, (5.0 - precision));
 
  266   easting = fmod (easting, 100000.0);
 
  267   if (easting >= 99999.5)
 
  269   east = (long)(easting/divisor);
 
  270   i += sprintf (USNGString+i, 
"%*.*ld", precision, precision, east);
 
  271   northing = fmod (northing, 100000.0);
 
  272   if (northing >= 99999.5)
 
  274   north = (long)(northing/divisor);
 
  275   i += sprintf (USNGString+i, 
"%*.*ld", precision, precision, north);
 
  304   while (USNGString[i] == 
' ')
 
  307   while (isdigit(USNGString[i]))
 
  315       strncpy (zone_string, USNGString+j, 2);
 
  317       sscanf (zone_string, 
"%ld", zone);
 
  318       if ((*zone < 1) || (*zone > 60))
 
  327   while (isalpha(USNGString[i]))
 
  330   if (num_letters == 3)
 
  333     letters[0] = (toupper(USNGString[j]) - (long)
'A');
 
  336     letters[1] = (toupper(USNGString[j+1]) - (long)
'A');
 
  339     letters[2] = (toupper(USNGString[j+2]) - (long)
'A');
 
  346   while (isdigit(USNGString[i]))
 
  349   if ((num_digits <= 10) && (num_digits%2 == 0))
 
  353     char north_string[6];
 
  362       strncpy (east_string, USNGString+j, n);
 
  364       sscanf (east_string, 
"%ld", &east);
 
  365       strncpy (north_string, USNGString+j+n, n);
 
  367       sscanf (north_string, 
"%ld", &north);
 
  368       multiplier = pow (10.0, 5.0 - n);
 
  369       *easting = east * multiplier;
 
  370       *northing = north * multiplier;
 
  388 USNG::USNG( 
double ellipsoidSemiMajorAxis, 
double ellipsoidFlattening, 
char* ellipsoidCode ) :
 
  403   double inv_f = 1 / ellipsoidFlattening;
 
  405    if (ellipsoidSemiMajorAxis <= 0.0)
 
  409   if ((inv_f < 250) || (inv_f > 350))
 
  417   strncpy (USNGEllipsoidCode, ellipsoidCode, 2);
 
  418    USNGEllipsoidCode[2] = 
'\0';
 
  428   ups = 
new UPS( *( u.ups ) );
 
  429   utm = 
new UTM( *( u.utm ) );
 
  433   strcpy( USNGEllipsoidCode, u.USNGEllipsoidCode );
 
  451     ups->operator=( *u.ups );
 
  452     utm->operator=( *u.utm );
 
  456     strcpy( USNGEllipsoidCode, u.USNGEllipsoidCode );
 
  500   double latitude  = geodeticCoordinates->
latitude();
 
  501   double longitude = geodeticCoordinates->
longitude();
 
  507   if ((longitude < -
PI) || (longitude > (2*
PI)))
 
  523         mgrsorUSNGCoordinates = fromUTM(
 
  524            utmCoordinates, longitude, latitude, precision );
 
  525         delete utmCoordinates;
 
  530         mgrsorUSNGCoordinates = fromUPS( upsCoordinates, precision );
 
  531         delete upsCoordinates;
 
  535      delete utmCoordinates;
 
  536      delete upsCoordinates;
 
  540   return mgrsorUSNGCoordinates;
 
  562   double usng_northing;
 
  570      &zone, letters, &usng_easting, &usng_northing, &precision );
 
  576         utmCoordinates = toUTM(
 
  577            zone, letters, usng_easting, usng_northing, precision );
 
  579         delete utmCoordinates;
 
  583         upsCoordinates = toUPS( letters, usng_easting, usng_northing );
 
  585         delete upsCoordinates;
 
  590      delete utmCoordinates;
 
  591      delete upsCoordinates;
 
  595   return geodeticCoordinates;
 
  616   long zone       = utmCoordinates->
zone();
 
  617   char hemisphere = utmCoordinates->
hemisphere();
 
  618   double easting  = utmCoordinates->
easting();
 
  619   double northing = utmCoordinates->
northing();
 
  621   if ((zone < 1) || (zone > 60))
 
  623   if ((hemisphere != 
'S') && (hemisphere != 
'N'))
 
  640   double latitude = geodeticCoordinates->
latitude();
 
  646         mgrsorUSNGCoordinates = fromUTM(
 
  647            utmCoordinates, geodeticCoordinates->
longitude(),
 
  648            latitude, precision );
 
  652         mgrsorUSNGCoordinates = fromUPS( upsCoordinates, precision );
 
  657      delete geodeticCoordinates;
 
  658      delete upsCoordinates;
 
  662   delete geodeticCoordinates;
 
  663   delete upsCoordinates;
 
  665   return mgrsorUSNGCoordinates;
 
  687   double usng_easting, usng_northing;
 
  694      &zone, letters, &usng_easting, &usng_northing, &precision );
 
  699         utmCoordinates = toUTM(
 
  700            zone, letters, usng_easting, usng_northing, precision );
 
  706         upsCoordinates = toUPS( letters, usng_easting, usng_northing );
 
  713      delete utmCoordinates;
 
  714      delete geodeticCoordinates;
 
  715      delete upsCoordinates;
 
  719   delete geodeticCoordinates;
 
  720   delete upsCoordinates;
 
  722   return utmCoordinates;
 
  745   char hemisphere = upsCoordinates->
hemisphere();
 
  746   double easting  = upsCoordinates->
easting();
 
  747   double northing = upsCoordinates->
northing();
 
  749   if ((hemisphere != 
'N') && (hemisphere != 
'S'))
 
  766   double latitude = geodeticCoordinates->
latitude();
 
  772         mgrsorUSNGCoordinates = fromUPS( upsCoordinates, precision );
 
  776         double longitude = geodeticCoordinates->
longitude();
 
  777         mgrsorUSNGCoordinates = fromUTM(
 
  778            utmCoordinates, longitude, latitude, precision );
 
  783      delete geodeticCoordinates;
 
  784      delete utmCoordinates;
 
  788   delete geodeticCoordinates;
 
  789   delete utmCoordinates;
 
  791   return mgrsorUSNGCoordinates;
 
  814   double usng_northing;
 
  821      &zone, letters, &usng_easting, &usng_northing, &precision );
 
  826         upsCoordinates = toUPS( letters, usng_easting, usng_northing );
 
  832         utmCoordinates = toUTM(
 
  833            zone, letters, usng_easting, usng_northing, precision );
 
  840      delete geodeticCoordinates;
 
  841      delete utmCoordinates;
 
  842      delete upsCoordinates;
 
  846   delete geodeticCoordinates;
 
  847   delete utmCoordinates;
 
  849   return upsCoordinates;
 
  871   double pattern_offset;      
 
  872   double grid_northing;       
 
  874   long ltr2_high_value;       
 
  880   long zone = utmCoordinates->
zone();
 
  881   char hemisphere = utmCoordinates->
hemisphere();
 
  882   double easting = utmCoordinates->
easting();
 
  883   double northing = utmCoordinates->
northing();
 
  885   getLatitudeLetter( latitude, &letters[0] );
 
  890     natural_zone = (long)(31 + ((longitude) / 
_6));
 
  892     natural_zone = (long)(((longitude) / 
_6) - 29);
 
  894   if (natural_zone > 60)
 
  897   if (zone != natural_zone) 
 
  902     UTMCoordinates* utmCoordinatesOverride = utmOverride.convertFromGeodetic(
 
  903        &geodeticCoordinates );
 
  905     zone       = utmCoordinatesOverride->
zone();
 
  906     hemisphere = utmCoordinatesOverride->
hemisphere();
 
  907     easting    = utmCoordinatesOverride->
easting();
 
  908     northing   = utmCoordinatesOverride->
northing();
 
  910     delete utmCoordinatesOverride;
 
  911     utmCoordinatesOverride = 0;
 
  917     if ((zone == 31) && (easting >= 
_500000))
 
  922      if ((zone == 32) && (easting < 
_500000)) 
 
  924      else if (((zone == 32) && (easting >= 
_500000)) || 
 
  925         ((zone == 34) && (easting < 
_500000))) 
 
  927      else if (((zone == 34) && (easting >= 
_500000)) || 
 
  928         ((zone == 36) && (easting < 
_500000))) 
 
  930      else if ((zone == 36) && (easting >= 
_500000)) 
 
  940         utmOverride.convertFromGeodetic( &geodeticCoordinates );
 
  942      zone       = utmCoordinatesOverride->
zone();
 
  943      hemisphere = utmCoordinatesOverride->
hemisphere();
 
  944      easting    = utmCoordinatesOverride->
easting();
 
  945      northing   = utmCoordinatesOverride->
northing();
 
  947      delete utmCoordinatesOverride;
 
  948      utmCoordinatesOverride = 0;
 
  952   double divisor = pow (10.0, (5.0 - precision));
 
  953   easting        = (long)(easting/divisor) * divisor;
 
  954   northing       = (long)(northing/divisor) * divisor;
 
  956   if( latitude <= 0.0 && northing == 1.0e7)
 
  962   getGridValues( zone, <r2_low_value, <r2_high_value, &pattern_offset );
 
  964   grid_northing = northing;
 
  966   while (grid_northing >= 
TWOMIL)
 
  968     grid_northing = grid_northing - 
TWOMIL;
 
  970   grid_northing = grid_northing + pattern_offset;
 
  971   if(grid_northing >= 
TWOMIL)
 
  972     grid_northing = grid_northing - 
TWOMIL;
 
  974   letters[2] = (long)(grid_northing / 
ONEHT);
 
  976     letters[2] = letters[2] + 1;
 
  979     letters[2] = letters[2] + 1;
 
  981   letters[1] = ltr2_low_value + ((long)(easting / 
ONEHT) -1);
 
  983     letters[1] = letters[1] + 1;
 
  985   makeUSNGString( USNGString, zone, letters, easting, northing, precision );
 
 1012   double min_northing;
 
 1013   double northing_offset;
 
 1014   long ltr2_low_value;
 
 1015   long ltr2_high_value;
 
 1016   double pattern_offset;
 
 1017   double upper_lat_limit;     
 
 1018   double lower_lat_limit;     
 
 1019   double grid_easting;        
 
 1020   double grid_northing;       
 
 1021   double temp_grid_northing = 0.0;
 
 1022   double fabs_grid_northing = 0.0;
 
 1023   double latitude  = 0.0;
 
 1024   double longitude = 0.0;
 
 1025   double divisor   = 1.0;
 
 1028   if((letters[0] == 
LETTER_X) && ((zone == 32) || (zone == 34) || (zone == 36)))
 
 1030   else if ((letters[0] == 
LETTER_V) && (zone == 31) && (letters[1] > 
LETTER_D))
 
 1039      getGridValues(zone, <r2_low_value, <r2_high_value, &pattern_offset);
 
 1044      if((letters[1] < ltr2_low_value) ||
 
 1045         (letters[1] > ltr2_high_value) ||
 
 1049      grid_easting = (double)((letters[1]) - ltr2_low_value + 1) * 
ONEHT;
 
 1051         grid_easting = grid_easting - 
ONEHT;
 
 1053      double row_letter_northing = (double)(letters[2]) * 
ONEHT;
 
 1055         row_letter_northing = row_letter_northing - 
ONEHT;
 
 1058         row_letter_northing = row_letter_northing - 
ONEHT;
 
 1060      if (row_letter_northing >= TWOMIL)
 
 1061         row_letter_northing = row_letter_northing - 
TWOMIL;
 
 1063      getLatitudeBandMinNorthing(letters[0], &min_northing, &northing_offset);
 
 1065      grid_northing = row_letter_northing - pattern_offset;
 
 1066      if(grid_northing < 0)
 
 1069      grid_northing += northing_offset;
 
 1071      if(grid_northing < min_northing)
 
 1074      easting = grid_easting + easting;
 
 1075      northing = grid_northing + northing;
 
 1079         zone, hemisphere, easting, northing );
 
 1089         delete utmCoordinates;
 
 1093      divisor = pow (10.0, (
double)precision);
 
 1094      getLatitudeRange(letters[0], &upper_lat_limit, &lower_lat_limit);
 
 1096      double latitude = geodeticCoordinates->
latitude();
 
 1098      delete geodeticCoordinates;
 
 1100      if(!(((lower_lat_limit - 
PI_OVER_180/divisor) <= latitude) &&
 
 1101          (latitude <= (upper_lat_limit + 
PI_OVER_180/divisor))))
 
 1105   return utmCoordinates;
 
 1125   double false_easting;      
 
 1126   double false_northing;     
 
 1127   double grid_easting;       
 
 1128   double grid_northing;      
 
 1129   long ltr2_low_value;       
 
 1133   char USNGString[21];
 
 1135   char hemisphere = upsCoordinates->
hemisphere();
 
 1136   double easting  = upsCoordinates->
easting();
 
 1137   double northing = upsCoordinates->
northing();
 
 1139   divisor  = pow (10.0, (5.0 - precision));
 
 1140   easting  = (long)(easting/divisor) * divisor;
 
 1141   northing = (long)(northing/divisor) * divisor;
 
 1143   if (hemisphere == 
'N')
 
 1145     if (easting >= TWOMIL)
 
 1150     index = letters[0] - 22;
 
 1157     if (easting >= TWOMIL)
 
 1167   grid_northing = northing;
 
 1168   grid_northing = grid_northing - false_northing;
 
 1169   letters[2] = (long)(grid_northing / ONEHT);
 
 1172     letters[2] = letters[2] + 1;
 
 1175     letters[2] = letters[2] + 1;
 
 1177   grid_easting = easting;
 
 1178   grid_easting = grid_easting - false_easting;
 
 1179   letters[1] = ltr2_low_value + ((long)(grid_easting / ONEHT));
 
 1181   if (easting < TWOMIL)
 
 1184       letters[1] = letters[1] + 3;
 
 1187       letters[1] = letters[1] + 2;
 
 1192       letters[1] = letters[1] + 2;
 
 1195       letters[1] = letters[1] + 1;
 
 1198       letters[1] = letters[1] + 3;
 
 1201   makeUSNGString( USNGString, 0, letters, easting, northing, precision );
 
 1208    long letters[USNG_LETTERS],
 
 1224   long ltr2_high_value;       
 
 1225   long ltr3_high_value;       
 
 1226   long ltr2_low_value;        
 
 1227   double false_easting;       
 
 1228   double false_northing;      
 
 1229   double grid_easting;        
 
 1230   double grid_northing;       
 
 1238     index = letters[0] - 22;
 
 1261   if ((letters[1] < ltr2_low_value) || (letters[1] > ltr2_high_value) ||
 
 1265       (letters[2] > ltr3_high_value))
 
 1268   grid_northing = (double)letters[2] * ONEHT + false_northing;
 
 1270     grid_northing = grid_northing - 
ONEHT;
 
 1273     grid_northing = grid_northing - 
ONEHT;
 
 1275   grid_easting = (double)((letters[1]) - ltr2_low_value) * ONEHT +false_easting;
 
 1279       grid_easting = grid_easting - 300000.0;
 
 1282       grid_easting = grid_easting - 200000.0;
 
 1287       grid_easting = grid_easting - 200000.0;
 
 1290       grid_easting = grid_easting - 
ONEHT;
 
 1293       grid_easting = grid_easting - 300000.0;
 
 1296   easting = grid_easting + easting;
 
 1297   northing = grid_northing + northing;
 
 1301      hemisphere, easting, northing);
 
 1305 void USNG::getGridValues(
 
 1307    long*   ltr2_low_value,
 
 1308    long*   ltr2_high_value,
 
 1309    double *pattern_offset )
 
 1326   set_number = zone % 6;
 
 1331   if ((set_number == 1) || (set_number == 4))
 
 1336   else if ((set_number == 2) || (set_number == 5))
 
 1341   else if ((set_number == 3) || (set_number == 6))
 
 1348   if ((set_number % 2) ==  0)
 
 1349     *pattern_offset = 500000.0;
 
 1351     *pattern_offset = 0.0;
 
 1355 void USNG::getLatitudeBandMinNorthing(
 
 1357    double* min_northing,
 
 1358    double* northing_offset )
 
 1390 void USNG::getLatitudeRange( 
long letter, 
double* north, 
double* south )
 
 1422 void USNG::getLatitudeLetter( 
double latitude, 
int* letter )
 
 1435   if (latitude >= 
_72 && latitude < 
_84_5)
 
 1437   else if (latitude > -
_80_5 && latitude < 
_72)
 
 1439     band = (long)(((latitude + 
_80) / 
_8) + 1.0e-12);