108                                                                     {
  109  
  110  
  111  
  112  
  113  
  114 
  115  const double PI = 3.14159265;
  116  
  117  const double deg2rad = PI / 180;
  118  
  119 
  122  double k0 = 0.9996;
  123 
  124  double LongOrigin;
  125  double eccPrimeSquared;
  126  double N, T, C, A, M;
  127 
  128  
  129  double LongTemp = (Long + 180) - int((Long + 180) / 360)*360 - 180; 
  130 
  131  double LatRad = Lat*deg2rad;
  132  double LongRad = LongTemp*deg2rad;
  133  double LongOriginRad;
  134  int ZoneNumber;
  135 
  136  ZoneNumber = int((LongTemp + 180) / 6) + 1;
  137 
  138  if (Lat >= 56.0 && Lat < 64.0 && LongTemp >= 3.0 && LongTemp < 12.0)
  139    ZoneNumber = 32;
  140 
  141  
  142  if (Lat >= 72.0 && Lat < 84.0) {
  143    if (LongTemp >= 0.0 && LongTemp < 9.0) ZoneNumber = 31;
  144    else if (LongTemp >= 9.0 && LongTemp < 21.0) ZoneNumber = 33;
  145    else if (LongTemp >= 21.0 && LongTemp < 33.0) ZoneNumber = 35;
  146    else if (LongTemp >= 33.0 && LongTemp < 42.0) ZoneNumber = 37;
  147  }
  148  LongOrigin = (ZoneNumber - 1)*6 - 180 + 3; 
  149  LongOriginRad = LongOrigin * deg2rad;
  150 
  151  
  153 
  154  eccPrimeSquared = (eccSquared) / (1 - eccSquared);
  155 
  156  N = a / sqrt(1 - eccSquared * sin(LatRad) * sin(LatRad));
  157  T = tan(LatRad) * tan(LatRad);
  158  C = eccPrimeSquared * cos(LatRad) * cos(LatRad);
  159  A = cos(LatRad)*(LongRad - LongOriginRad);
  160 
  161  M = a * ((1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256) * LatRad
  162           - (3 * eccSquared / 8 + 3 * eccSquared * eccSquared / 32 + 45 * eccSquared * eccSquared * eccSquared / 1024) * sin(2 * LatRad)
  163           + (15 * eccSquared * eccSquared / 256 + 45 * eccSquared * eccSquared * eccSquared / 1024) * sin(4 * LatRad)
  164           - (35 * eccSquared * eccSquared * eccSquared / 3072) * sin(6 * LatRad));
  165 
  166  UTMEasting = (double) (k0 * N * (A + (1 - T + C) * A * A * A / 6
  167                                   + (5 - 18 * T + T * T + 72 * C - 58 * eccPrimeSquared) * A * A * A * A * A / 120)
  168                         + 500000.0);
  169 
  170  UTMNorthing = (double) (k0 * (M + N * tan(LatRad)*(A * A / 2 + (5 - T + 9 * C + 4 * C * C) * A * A * A * A / 24
  171                                                     + (61 - 58 * T + T * T + 600 * C - 330 * eccPrimeSquared) * A * A * A * A * A * A / 720)));
  172  if (Lat < 0)
  173    UTMNorthing += 10000000.0; 
  174}
static Ellipsoid ellipsoid[]
 
char UTMLetterDesignator(double Lat)
 
double eccentricitySquared