# Last Modified: January 8, 2003
# Description:
#     This code converts latitude/longitude coordinates into UTM cordinates
#		 using the GRS80 ellipsoid (North American Datum 1983).
#
# History: Based on the National Defense Mapping Agency Program
#				Translated from fortran to TrueBasic by Chuck Connor, 3-3-87
#				Translated to MAC, by Chuck Connor, 5-24-89
#				Translated to C, by Laura Connor, 12-97
#           Translated to Perl by Chuck Connor, 1-03 ...for the love of god, will this never end.....
# Example Usage:                                                                          
# @latlong=(27,-73.5);
# @utm=latlong2utm(@latlong);
# print "@utm\n";
#

sub latlong2utm {
  
  my @a=();
  my @b1= ();
  my ($fo, $r, $z);
  my ($lon, $lat, $lat1, $x4);
  my ($easting, $northing);
  my ($bo, $so, $co, $to, $t9, $eo);
  
  ($lat, $lon) = @_;
 
#$a[14] is the equatorial radius (NAD83)
$a[14] = 6378137.0;

#$a[15] is the first eccentricity^2 (NAD83)
$a[15] = .0066943800229;

#here are the values for a Clarke 1866 spheroid (NAD27 datum)
#$a[14]=6378206.4;
#$a[15]=0.00676866;

#here are the values for WGS84
#$a[14]=6378137.0;
#$a[15]=0.00669437999013;

#start grinding....
  
  $a[9] = ((($a[15] * (7.0/32.0) + 5.0/16.0) * $a[15] + .5) * $a[15] + 1.0 ) * $a[15] * .25;
  $a[0] = -(((($a[9] * (195.0/64.0) + 3.25) * $a[9] + 3.75) * $a[9] + 3.0) * $a[9]);
  $a[1] = ((1455.0/32.0 * $a[9] + 70.0/3.0) * $a[9] + 7.5) * $a[9] * $a[9];
  $a[2] = -((70.0/3.0 + $a[9] * (945.0/8.0)) * $a[9] * $a[9] * $a[9]);
  $a[3] = 315.0/4.0 * $a[9] * $a[9] * $a[9] * $a[9];
  $a[10] = (((7.75 - 657.0/64.0 * $a[9]) * $a[9] - 5.25) * $a[9] + 3.0) * $a[9];
  $a[11] = ((5045.0/32.0 * $a[9] - 151.0/3.0) * $a[9] + 10.5) * $a[9] * $a[9];
  $a[12] = (151.0/3.0 - 3291.0/8.0 * $a[9]) * $a[9] * $a[9] * $a[9];
  $a[13] = 1097.0/4.0 * $a[9] * $a[9] * $a[9] * $a[9];
  
  $fo = $a[9] * $a[9];
  $a[9] = ((225.0/64.0 * $fo + 2.25) * $fo + 1.0) * (1.0 - $fo) * (1.0 - $a[9]) * $a[14];
  $r = 206264.806247;
  $a[4] = 500000.0;
  $a[7] = .9996;
  $a[6] = 0.0;
  $a[5] = 0.0;
  
  $z = 0;
  
  $lon -= 360.0;
  $lon *= 3600.0;
  
  $lat *= 3600.0;
  
  #Compute the latitude 
  $x4 = abs($lat);
  $lat1=int($x4/3600.0);

print "lat1 = $lat1\n";
print "lat = $lat\n";
  if ($lat < 0.0) {$lat1 = -$lat1;}
  
  
  # Compute the longitude 
  $x4 = abs($lon);
  
  if (abs($lat) > 302400.0)
    { 
      $easting = 0.0;
      $northing = 0.0;
      
    }
  else
    { $z = int(30.0 + ($lon/21600.0));
      
      if ($z > 30) {$a[8] = ( ($z - 30.0) * 6.0 - 3.0) * 3600.0;}
      else {$a[8] = ((30.0 - $z) * 6.0 + 3.0) * (-3600.0);}
      
      $b1[9] = ($lon - $a[8])/$r;
      
      $b1[8] = $lat / $r;
      $so = sin($b1[8]);
      $co = cos($b1[8]);
      $ro = $a[14] / (sqrt(1.0 - $a[15] * $so * $so));
      $to = $so / $co;
      $t9 = $to * $to;
      $b1[10] = $co * $co;
      $eo = $a[15] * $b1[10] / (1.0 - $a[15]);
      $b1[0] = $ro * $co;
      $b1[2] = (1.0 - $t9 + $eo) * $b1[0] * $b1[10] / 6.0;
      $b1[4] = (($t9 - 18.0) * $t9 + 5.0 + (14.0 - 58.0 * $t9) * $eo) * $b1[0] * $b1[10] * $b1[10] / 120.0;
      $b1[6] = (((179.0 - $t9) * $t9 - 479.0) * $t9 + 61.0) * $b1[0] * $b1[10] * $b1[10] * $b1[10] / 5040.0;
      $b1[11] = $b1[9] * $b1[9];
      
      # Calculate the easting 
      $easting = ((($b1[6] * $b1[11] + $b1[4]) * $b1[11] + $b1[2]) * $b1[11] + $b1[0]) * $b1[9] * $a[7] + $a[4];
      
      $b1[1] = $ro * $b1[10] * $to / 2.0;
      $b1[3] = ($eo * (9.0 + 4.0 * $eo) + 5.0 - $t9) * $b1[1] * $b1[10] / 12.0;
      $b1[5] = (($t9 - 58.0) * $t9 + 61.0 + (270.0 - 330.0 * $t9) * $eo) * $b1[1] * $b1[10] * $b1[10] / 360.0;
      $b1[7] = (((543.0 - $t9) * $t9 - 3111.0) * $t9 + 1385.0) * $b1[1] * $b1[10] * $b1[10] * $b1[10] / 20160.0;
      
      #Calculate the northing
      $northing = ((($b1[7] * $b1[11] + $b1[5]) * $b1[11] + $b1[3]) * $b1[11] + $b1[1]) * $b1[11];
      $northing += (((($a[3] * $b1[10] + $a[2]) * $b1[10] + $a[1]) * $b1[10] + $a[0]) * $so * $co + $b1[8]) * $a[9];
      $northing = ($northing - $a[6]) * $a[7] + $a[5];
      
      if ($lat < 0) {$northing += 10000000.0;}
    }

return $northing, $easting;
}