#!/usr/bin/perl -W
#  use Image::Magick;
use POSIX;
# Version 1.0 21/12/2006
# Version 1.1 14/3/2007 web form
use CGI;


##vars
#
my $maj1830=6377563.396; # Airy1830
my $min1830=6356256.909;
my $maj84=6378137.0; # GRS80 aka WGS84
my $min84=6356752.3141;
my $E1830=0.00667054007414908; # ref ellipsoid 1830
my $E84=0.00669438003551284; # ref ellipsoid 84
my $pi=3.14159265358979;
my @ns=("S","N");
my @ew=("W","E");
my $mapgen="lmap.htm";
my $mapjoin="lmap.bat";
my $mapinfo="lmap.inf";
my $leadin="montage -geometry +0+0 -tile "; # 13x8  NT07NE08.gif
my $leadout=" lmap.gif";
my ($topleft,$g);
my ($East,$North);
my ($Lat,$Long);
my ($UTMEasting, $UTMNorthing, $UTMZone, $longitudeZone);
my $gifs="";
my $mapname="Map Name";
#
##eovars

my $ostl='NN200300';
my $osbr='NT500300';
my $errortext='';

##check incoming, make ostl and osbr be what came in subject to error checks
$query = new CGI;
$ostlin = $query->param('ostlin');
if ($ostlin) {
  if ($ostlin =~/[A-Z]{2}\d{6}/) {
    $ostl=$ostlin;
  } else {
    $errortext.="Expecting OSGB top left in form NT230450!<br>";
  }
}

$osbrin = $query->param('osbrin');
if ($osbrin) {
  if ($osbrin =~/[A-Z]{2}\d{6}/) {
    $osbr=$osbrin;
  } else {
    $errortext.="Expecting OSGB bottom right in form NT450230!<br>";
  }
}


my $header="<html><head><meta http-equiv='Content-Type' content='text/html; charset=windows-1252'>
<title>OSGB to Lat/Long for UI-View32 Maps</title>
</head><body bgcolor='#3399FF' bgproperties='fixed' background='../../img/back.jpg'>
<table border=0 width=100%>
<tr><td align=left><b><a href='../../'>Back Home</a> -> <a href='../'>Back to Project Index</a> -&gt; <a href='index.html'>Back to APRS Index</a></b></td>
<td align=right><b><a href='http://graham.auld.me.uk'>http://graham.auld.me.uk</a></b></td></tr>
</table>
<h2>Convert OSGB to Decimal Lat & Lon</h2>\n";

my $footer="<hr>Author: <i>osgb2ll.pl v1.1</i><br>Dated: <i>14th March 2007</i></body>\n</html>\n";

print "Content-type: text/html\nCache-Control: max-age=1\n\n";
print "<META HTTP-EQUIV=\"expires\" CONTENT=\"Wed, 26 Feb 2003 08:21:57 GMT\">\n";
print $header;
if ($errortext ne '') {
	print "<font color=red><font size=+2>Error!</font><br>$errortext</font><br>\n";
}
print "<form method=POST>\n";

print "<table border=0><tr>
<td><b>Input OS Grid Corners</b></td><td width=20></td>
<td><b>Output for UI-View32 inf file</b></td>
</tr><tr>
<td background=mapback.png width=200 height=200>";
#input side
print "<table border=0 width=100% height=100%><tr><td valign=top><input name=ostlin type=text value=$ostl size=10><br>OSGB Top Left</td></tr><tr><td valign=bottom align=right>OSGB Bottom Right<br><input value=$osbr name=osbrin type=text size=10></td></tr></table>";
print "</td><td width=20><input type=submit value='Convert -->'</td>";

#output side
print "<td valign=top><br>";



print "<b>mapname.inf</b>:<br>\n";
print "<textarea rows=3 cols=20>";
genmap ($ostl,$osbr);
print "</textarea>";

print "</td></tr></table>";

print "<BR>Want to find the OSGB Grid reference of your map image?<P>Then go to <a href=\"http://www.streetmap.co.uk\">http://www.streetmap.co.uk</a>, find the map area of interest which matches the extent of your own image, clicking on the required position, which will then be marked by an orange coloured arrow. Then look for <B>\"Click <a href=\"http://www.streetmap.co.uk/streetmap.dll?GridConvert\">here</a> to convert/measure coordinates\"</B> The link will take you to a Grid Conversion Results page. The 'LR' row in the table will provide you with the OS Grid reference of your selected position. Enter the OSGB coordinates into the the appropriate corner of the map above. After pressing the convert button, a fully formatted UI-View32 inf file will be displayed on the right hand side. Just copy and paste this into a <B>.inf</B> file in the UI-View32 map area, along with a map image file of the same name.";
print "<P>Note: You are fitting a 2-dimentional map to the 3-dimentional curvature of the earth. On small maps there is little problem but on larger maps, the deviation becomes more obvious. Thus while the identified corners and its major diagonal will be correct, the more distant areas accumulate a greater degree of error. Adjustments may be required to reduce the errors in the primary area of interest.";

print $footer;


###theend###
exit();

############below here is conversion stuff#################
# call like so to generate some text as should appear in .inf file...
# genmap ("NN200300","NT500300");
#





#  my $pos="NT090790";
  #$pos="NT500300";
  #$pos="NT080770";
  #getOSRefFromSixFigureReference($pos);
  #print "Ref $East, $North\n";
  #$ngr=toSixFigureString($East,$North);
  #print "$pos = $ngr\n";
  #exit(0);
  #testpos("NT080770"); # eastings 308000, 677000 - 55.58.632N 3.26.467W for APRS mapping
  #exit(0);



#genmap ("NN200300","NT500300");

#
# Given the 6-figure NGR reference of the top left and bottom right squares, it will generate a 'lmap.htm' page
# which allows all the 'Streetmap.co.uk' 1k x 1k images to be displayed in a large table.
# It will then generate a 'lmap.bat' file which can subsequentially run to stitch all the 1k x 1k images together
# into one large 'lmap.gif' but all the required images must be available in the local folder first, hence check the HTM display.
# A 'lmap.inf' file will be generated which, after name editing, can be used with UI-View32.
##
sub genmap {
  my($fromref,$toref)=@_;

  getOSRefFromSixFigureReference($toref); # bottom right of area
  my $gx1=$East;
  my $gy1=$North;

  getOSRefFromSixFigureReference($fromref); # top left of area
  my $gx=$East;
  my $gy=$North;
#  print "$fromref to $toref - ($gx,$gy) to ($gx1,$gy1)\n";
#  #conversion from...

#  open (H, ">$mapgen") or die "Can't open $mapgen $!"; # gen web file to check map
#  print H "<http>\n<body><table cellspacing=0 cellpadding=0>\n";
#  my $gifs='';

#  for ($j=$gy;$j>=$gy1;$j-=10000) { # top left down
#    print H "<tr>";

#    for ($i=$gx;$i<=$gx1;$i+=10000) { # top left right
      #print "$i,$j ";

      #     $ngr=toSixFigureString($i,$j);
      #my $s1=substr($ngr, 0, 2);
#      my $s2=substr($ngr, 2, 1);
#      my $s2a=$s2;
#      my $s3=substr($ngr, 5, 1);

#      my $smngr=$s1.$s2.$s3;

#      print "$smngr ";
#      print H "<td><img src=\"..\\squares\\$smngr.gif\"</td>";
#      $gifs.=' ..\\squares\\'.$smngr.'.gif';
         
#    }
#    print "\n";
#    print H "</tr>\n";
#  }
#  print H "</table></body>\n</http>\n";
#  close(H);

#  my $hor=($gx1-$gx)/10000+1;
#  my $ver=($gy-$gy1)/10000+1;

#  open (B, ">$mapjoin") or die "Can't open $mapjoin $!"; # batch file to gen full map
#  print B $leadin.$hor.'x'.$ver.$gifs.$leadout."\n";
#  close(B);

  #my $m1=OSreftoDEClatlong($fromref,0,10000-($ver*0.5)); # middle left
  #my $m2=OSreftoDEClatlong($fromref,($hor*0.5),10000); # middle top
  #my $m3=OSreftoDEClatlong($toref,10000,($ver*0.5)); # middle right
  #my $m4=OSreftoDEClatlong($toref,10000-($hor*0.5),0); # middle bottom


  # eastings are ~100m out
  $s4=OSreftoDEClatlong($fromref,0,10000); # move to top left of first square

  toOSRef($Lat,$Long); # Convert to Lat/Long
  $East=int($East+0.5); $North=int($North+0.5);
  $s6=toSixFigureString($East,$North);
  #print "$s6 ($East,$North)";
  $s5=OSreftoDEClatlong($toref,10000,0); # move bottom right of last square

  toOSRef($Lat,$Long); # Convert to Lat/Long
  $East=int($East+0.5); $North=int($North+0.5);
  $s6=toSixFigureString(($East),$North);
  #print " - $s6 ($East,$North)\n";
  #print "$s4 - $s5\n";
  #open (R, ">$mapinfo") or die "Can't open $mapinfo $!"; # APRS map info
#  print R "$s4\n$s5\nNew LMap\n";
  print "$s4\n$s5\n$mapname\n";
  #close(R);
#  print $footer; 
}
#
# Position check routine
#
sub testpos {
  my ($pos)=@_;
  $pos=OSreftoDEClatlong($pos,0,0);
  printf "$pos for APRS mapping\n";
}

# ------------------------------

sub deg2rad { my ($deg)=@_; return $deg*$pi/180 }

sub rad2deg { my ($rad)=@_; return $rad*180/$pi }

sub sinSquared { my ($x)=@_; return sin($x)*sin($x) }

sub cosSquared { my ($x)=@_; return cos($x)*cos($x) }

sub tanSquared { my ($x)=@_; return tan($x)*tan($x) }

sub sec { my ($x)=@_; return 1.0 / cos($x) }
#
# numeric latitide, longditude to decimal minutes for APRS
#
sub OSreftoDEClatlong {
  my ($pos,$xdelta,$ydelta)=@_;
  getOSRefFromSixFigureReference($pos); # Alpha-numeric pos to fully numeric
  $East+=$xdelta;
  $North+=$ydelta;

  gridtoLatLng($East,$North); # Convert to Lat/Long
  #OSGB36ToWGS84($Lat,$Long); # Geodetic Transform GB National grid to World Geodetic system
  my $NS='N';
  my $EW; my $lo2;
  my $la1=int($Lat);
  my $la2=($Lat-$la1)*60;
  if ($Long<0.0) {
    $EW='W'; $lo2=-$Long;
  } else {
    $EW='E'; $lo2=$Long;
  }
  my $lo1=int($lo2);
  $lo2=($lo2-$lo1)*60;
  $lo1= sprintf "%d.%.3f%s,%d.%.3f%s",$la1,$la2,$NS,$lo1,$lo2,$EW;
  return $lo1;
}


# Convert this Latitude Longditude from OSGB36 datum to WGS84 datum.
# input numeric (latitude, longditude)
# result to global $Lat, $Long
#
sub OSGB36ToWGS84 {
  my ($lat,$lng)=@_;
  my $a = $maj1830;
  my $b = $min1830;
  my $eSquared = (($a * $a) - ($b * $b)) / ($a * $a);
     $phi = deg2rad($lat);
     $lambda = deg2rad($lng);
     $v = $a / (sqrt(1 - $eSquared * sinSquared($phi)));
     $H = 0; # height
  my $x = ($v + $H) * cos($phi) * cos($lambda);
  my $y = ($v + $H) * cos($phi) * sin($lambda);
  my $z = ((1 - $eSquared) * $v + $H) * sin($phi);

  my  $tx = 446.448;
#  my  $ty = -124.157;
  my $ty = -125.157;
  my $tz = 542.060;
  my $s  = -0.0000204894;
  my $rx = deg2rad( 0.00004172222);
  my $ry = deg2rad( 0.00006861111);
  my $rz = deg2rad( 0.00023391666);

  my $xB = $tx + ($x * (1 + $s)) + (-$rx * $y) + ($ry * $z);
  my $yB = $ty + ($rz * $x) + ($y * (1 + $s)) + (-$rx * $z);
  my $zB = $tz + (-$ry * $x) + ($rx * $y) + ($z * (1 + $s));

  $a = $maj84;
  $b = $min84;
  $eSquared = (($a * $a) - ($b * $b)) / ($a * $a);

  my $lambdaB = rad2deg(atan($yB / $xB));
  my $p = sqrt(($xB * $xB) + ($yB * $yB));
  my $phiN = atan($zB / ($p * (1 - $eSquared)));
  for ($i = 1; $i < 10; $i++) {
    my $v = $a / (sqrt(1 - $eSquared * sinSquared($phiN)));
    $phiN1 = atan(($zB + ($eSquared * $v * sin($phiN))) / $p);
    $phiN = $phiN1;
  }

  my $phiB = rad2deg($phiN);
      
  $Lat = $phiB;
  $Long = $lambdaB;
}
#
# Convert this grid reference into a latitude and longitude
# input numeric (Eastings, Northings)
# result to global $Lat, $Long
#
sub gridtoLatLng {
  my ($E,$N)=@_;
  my $OSGB_F0 = 0.9996012717; # National grid
  my $N0 = -100000.0;
  my $E0 = 400000.0;
  my $phi0 = deg2rad(49.0);
  my $lambda0 = deg2rad(-2.0);
  my $a = $maj1830;
  my $b = $min1830;
  my $eSquared = (($a * $a) - ($b * $b)) / ($a * $a);
  my $phi = 0.0;
  my $lambda = 0.0;
  my $n = ($a - $b) / ($a + $b);
  my $M = 0.0;
  my $phiPrime = (($N - $N0) / ($a * $OSGB_F0)) + $phi0;
  do {
    $M = ($b * $OSGB_F0) * (((1 + $n + ((5.0 / 4.0) * $n * $n) + ((5.0 / 4.0) * $n * $n * $n)) * ($phiPrime - $phi0))
    - (((3 * $n) + (3 * $n * $n) + ((21.0 / 8.0) * $n * $n * $n)) * sin($phiPrime - $phi0) * cos($phiPrime + $phi0))
    + ((((15.0 / 8.0) * $n * $n) + ((15.0 / 8.0) * $n * $n * $n)) * sin(2.0 * ($phiPrime - $phi0))
    * cos(2.0 * ($phiPrime + $phi0))) - (((35.0 / 24.0) * $n * $n * $n) * sin(3.0 * ($phiPrime - $phi0)) * cos(3.0 * ($phiPrime + $phi0))));
    $phiPrime += ($N - $N0 - $M) / ($a * $OSGB_F0);
  } while (($N - $N0 - $M) >= 0.001);
  my $v = $a * $OSGB_F0 * pow(1.0 - $eSquared * sinSquared($phiPrime), -0.5);
  my $rho = $a * $OSGB_F0 * (1.0 - $eSquared) * pow(1.0 - $eSquared * sinSquared($phiPrime), -1.5);
  my $etaSquared = ($v / $rho) - 1.0;
  my $VII = tan($phiPrime) / (2 * $rho * $v);
  my $VIII = (tan($phiPrime) / (24.0 * $rho * pow($v, 3.0))) * (5.0 + (3.0 * tanSquared($phiPrime)) + $etaSquared - (9.0 * tanSquared($phiPrime) * $etaSquared));
  my $IX = (tan($phiPrime) / (720.0 * $rho * pow($v, 5.0))) * (61.0 + (90.0 * tanSquared($phiPrime)) + (45.0 * tanSquared($phiPrime) * tanSquared($phiPrime)));
  my $X = sec($phiPrime) / $v;
  my $XI = (sec($phiPrime) / (6.0 * $v * $v * $v)) * (($v / $rho) + (2 * tanSquared($phiPrime)));
  my $XII = (sec($phiPrime) / (120.0 * pow($v, 5.0))) * (5.0 + (28.0 * tanSquared($phiPrime)) + (24.0 * tanSquared($phiPrime) * tanSquared($phiPrime)));
  my $XIIA = (sec($phiPrime) / (5040.0 * pow($v, 7.0))) * (61.0 + (662.0 * tanSquared($phiPrime))
    + (1320.0 * tanSquared($phiPrime) * tanSquared($phiPrime)) + (720.0 * tanSquared($phiPrime) * tanSquared($phiPrime) * tanSquared($phiPrime)));
   $phi = $phiPrime - ($VII * pow($E - $E0, 2.0)) + ($VIII * pow($E - $E0, 4.0)) - ($IX * pow($E - $E0, 6.0));
   $lambda = $lambda0 + ($X * ($E - $E0)) - ($XI * pow($E - $E0, 3.0)) + ($XII * pow($E - $E0, 5.0)) - ($XIIA * pow($E - $E0, 7.0));
    $Lat=rad2deg($phi); # global lat return
    $Long=rad2deg($lambda); # global long return
 }
#
# Convert this UTM reference to a latitude and longitude
# input numeric (Eastings, Northings)
# result to global $Lat, $Long
#
sub UTMtoLatLng {
  my ($x,$y,$zoneNumber,$zoneLetter)=@_;
  $x-=500000.0;
  my $UTM_F0   = 0.9996;
  my $a=$maj84;
  my $b=$min84;
  my $eSquared = (($a * $a) - ($b * $b)) / ($a * $a);
  my $ePrimeSquared = $eSquared / (1.0 - $eSquared);
  my $e1 = (1 - sqrt(1 - $eSquared)) / (1 + sqrt(1 - $eSquared));
      #$x = $this->easting - 500000.0;;
      #$y = $this->northing;
      #$zoneNumber = $this->lngZone;
      #$zoneLetter = $this->latZone;

  my $longitudeOrigin = ($zoneNumber - 1.0) * 6.0 - 180.0 + 3.0;

  # Correct y for southern hemisphere
  if ((ord($zoneLetter) - ord("N")) < 0) { $y -= 10000000.0; }

  my $m = $y / $UTM_F0;
  my $mu = $m / ($a * (1.0 - $eSquared / 4.0 - 3.0 * $eSquared * $eSquared / 64.0 - 5.0 * pow($eSquared, 3.0) / 256.0));

  my $phi1Rad = $mu + (3.0 * $e1 / 2.0 - 27.0 * pow($e1, 3.0) / 32.0) * sin(2.0 * $mu)
    + (21.0 * $e1 * $e1 / 16.0 - 55.0 * pow($e1, 4.0) / 32.0) * sin(4.0 * $mu) + (151.0 * pow($e1, 3.0) / 96.0) * sin(6.0 * $mu);

  my $n = $a / sqrt(1.0 - $eSquared * sin($phi1Rad) * sin($phi1Rad));
  my $t = tan($phi1Rad) * tan($phi1Rad);
  my $c = $ePrimeSquared * cos($phi1Rad) * cos($phi1Rad);
  my $r = $a * (1.0 - $eSquared) / pow(1.0 - $eSquared * sin($phi1Rad) * sin($phi1Rad), 1.5);
  my $d = $x / ($n * $UTM_F0);

  my $latitude = ($phi1Rad - ($n * tan($phi1Rad) / $r) * ($d * $d / 2.0 - (5.0 + (3.0 * $t) + (10.0 * $c)
    - (4.0 * $c * $c) - (9.0 * $ePrimeSquared)) * pow($d, 4.0) / 24.0 + (61.0 + (90.0 * $t) + (298.0 * $c)
    + (45.0 * $t * $t) - (252.0 * $ePrimeSquared) - (3.0 * $c * $c)) * pow($d, 6.0) / 720.0)) * (180.0/$pi);

  my $longitude = $longitudeOrigin + (($d - (1.0 + 2.0 * $t + $c) * pow($d, 3.0)/6.0 + (5.0 - (2.0 * $c)
    + (28.0 * $t) - (3.0 * $c * $c) + (8.0 * $ePrimeSquared) + (24.0 * $t * $t)) * pow($d, 5.0) / 120.0) / cos($phi1Rad)) * (180.0/$pi);

  $Lat=$latitude; # global Lat return
  $Long=$longitude; # global Long return
}
#
# Take a string formatted as a six-figure OS grid reference (e.g.
# "TG514131") and return a reference to an OSRef object that represents
# that grid reference. The first character must be H, N, S, O or T.
# The second character can be any uppercase character from A through Z
# excluding I.
# input map reference string (ref)
# result to global $East, $North
#
sub getOSRefFromSixFigureReference {
  my ($ref)=@_;
  my $char1 = substr($ref, 0, 1);
  my $char2 = substr($ref, 1, 1);
  my $east  = int(substr($ref, 2, 3)) * 100;
  my $north = int(substr($ref, 5, 3)) * 100;
  if ($char1 eq 'H') { $north += 1000000;
  } elsif ($char1 eq 'N') { $north += 500000;
  } elsif ($char1 eq 'O') { $north += 500000; $east += 500000;
  } elsif ($char1 eq 'T') { $east += 500000;
  }
  my $char2ord = ord($char2);
  if ($char2ord > 73) {$char2ord--;} # Adjust for no I
  my $nx = (($char2ord - 65) % 5) * 100000;
  my $ny = (4 - floor(($char2ord - 65) / 5)) * 100000;
  $East=$east + $nx; # global OSRef return
  $North=$north + $ny; # global OSRef return
}
#
# Convert this grid reference into a string using a standard six-figure
# grid reference including the two-character designation for the 100km
sub toSixFigureString {
  my ($easting,$northing)=@_;
  my $hundredkmE=floor($easting / 100000);
  my $hundredkmN=floor($northing / 100000);
  my $firstLetter="";
  if ($hundredkmN < 5) {
    if ($hundredkmE < 5) {
      $firstLetter="S";
    } else {
      $firstLetter="T";
    }
  } elsif ($hundredkmN < 10) {
    if ($hundredkmE < 5) {
      $firstLetter="N";
    } else {
      $firstLetter="O";
    }
  } else {
    $firstLetter="H";
  }

  my $secondLetter="";
  my $index=65 + ((4 - ($hundredkmN % 5)) * 5) + ($hundredkmE % 5);
  my $ti=$index;
  if ($index >=73) {$index++};
  $secondLetter=chr($index);

  my $e=floor(($easting - (100000 * $hundredkmE)) / 100);
  my $n=floor(($northing - (100000 * $hundredkmN)) / 100);
  my $es=$e;
  if ($e < 100) {$es="0$es"};
  if ($e < 10)  {$es="0$es"};
  my $ns=$n;
  if ($n < 100) {$ns="0$ns"};
  if ($n < 10)  {$ns="0$ns"};

  return $firstLetter . $secondLetter . $es . $ns;
}

#
# Work out the UTM latitude zone from the latitude
# input latitude (latitude)
# return from subroutine zone letter
#
sub getUTMLatitudeZoneLetter {
  my ($latitude)=@_;
  if ((84 >= $latitude) && ($latitude >= 72)) {return "X"}
  elsif (( 72 > $latitude) && ($latitude >=  64)) {return "W"}
  elsif (( 64 > $latitude) && ($latitude >=  56)) {return "V"}
  elsif (( 56 > $latitude) && ($latitude >=  48)) {return "U"}
  elsif (( 48 > $latitude) && ($latitude >=  40)) {return "T"}
  elsif (( 40 > $latitude) && ($latitude >=  32)) {return "S"}
  elsif (( 32 > $latitude) && ($latitude >=  24)) {return "R"}
  elsif (( 24 > $latitude) && ($latitude >=  16)) {return "Q"}
  elsif (( 16 > $latitude) && ($latitude >=   8)) {return "P"}
  elsif ((  8 > $latitude) && ($latitude >=   0)) {return "N"}
  elsif ((  0 > $latitude) && ($latitude >=  -8)) {return "M"}
  elsif (( -8 > $latitude) && ($latitude >= -16)) {return "L"}
  elsif ((-16 > $latitude) && ($latitude >= -24)) {return "K"}
  elsif ((-24 > $latitude) && ($latitude >= -32)) {return "J"}
  elsif ((-32 > $latitude) && ($latitude >= -40)) {return "H"}
  elsif ((-40 > $latitude) && ($latitude >= -48)) {return "G"}
  elsif ((-48 > $latitude) && ($latitude >= -56)) {return "F"}
  elsif ((-56 > $latitude) && ($latitude >= -64)) {return "E"}
  elsif ((-64 > $latitude) && ($latitude >= -72)) {return "D"}
  elsif ((-72 > $latitude) && ($latitude >= -80)) {return "C"}
  else {return 'Z'}
}
#
# Convert this LatLng object into an OSGB grid reference. Note that this
# function does not take into account the bounds of the OSGB grid -
# beyond the bounds of the OSGB grid, the resulting OSRef object has no
# meaning
#
sub toOSRef {
  my ($lat,$long)=@_;
  my $OSGB_F0=0.9996012717;
  my $N0=-100000.0;
  my $E0=400000.0;
  my $phi0=deg2rad(49.0);
  my $lambda0=deg2rad(-2.0);
  my $a=$maj1830;
  my $b=$min1830;
  my $eSquared=$E1830;
  my $phi=deg2rad($lat);
  my $lambda=deg2rad($long);
  my $E=0.0;
  my $N=0.0;
  my $n=($a - $b) / ($a + $b);
  my $v=$a * $OSGB_F0 * pow(1.0 - $eSquared * sinSquared($phi), -0.5);
  my $rho=$a * $OSGB_F0 * (1.0 - $eSquared) * pow(1.0 - $eSquared * sinSquared($phi), -1.5);
  my $etaSquared=($v / $rho) - 1.0;
  my $M=($b * $OSGB_F0) * (((1 + $n + ((5.0 / 4.0) * $n * $n) + ((5.0 / 4.0) * $n * $n * $n)) * ($phi - $phi0))
       - (((3 * $n) + (3 * $n * $n) + ((21.0 / 8.0) * $n * $n * $n)) * sin($phi - $phi0) * cos($phi + $phi0))
       + ((((15.0 / 8.0) * $n * $n) + ((15.0 / 8.0) * $n * $n * $n)) * sin(2.0 * ($phi - $phi0)) * cos(2.0 * ($phi + $phi0)))
       - (((35.0 / 24.0) * $n * $n * $n) * sin(3.0 * ($phi - $phi0)) * cos(3.0 * ($phi + $phi0))));
  my $I=$M + $N0;
  my $II=($v / 2.0) * sin($phi) * cos($phi);
  my $III=($v / 24.0) * sin($phi) * pow(cos($phi), 3.0) * (5.0 - tanSquared($phi) + (9.0 * $etaSquared));
  my $IIIA=($v / 720.0) * sin($phi) * pow(cos($phi), 5.0) * (61.0 - (58.0 * tanSquared($phi)) + pow(tan($phi), 4.0));
  my $IV=$v * cos($phi);
  my $V=($v / 6.0) * pow(cos($phi), 3.0) * (($v / $rho) - tanSquared($phi));
  my $VI=($v / 120.0) * pow(cos($phi), 5.0) * (5.0 - (18.0 * tanSquared($phi)) + (pow(tan($phi), 4.0)) + (14 * $etaSquared) - (58 * tanSquared($phi) * $etaSquared));
  $North=$I + ($II * pow($lambda - $lambda0, 2.0)) + ($III * pow($lambda - $lambda0, 4.0)) + ($IIIA * pow($lambda - $lambda0, 6.0));
  $East=$E0 + ($IV * ($lambda - $lambda0)) + ($V * pow($lambda - $lambda0, 3.0)) + ($VI * pow($lambda - $lambda0, 5.0));
  # Global North, East return
}
# Convert a latitude and longitude to an UTM reference
#
#
sub toUTMRef {
  my ($latitude,$longitude)=@_;
  my $UTM_F0  =0.9996;
  my $a=$maj84;
  my $eSquared=$E84;

  my $latitudeRad=$latitude * ($pi / 180.0);
  my $longitudeRad=$longitude * ($pi / 180.0);
  $longitudeZone=int(($longitude + 180.0) / 6.0) + 1;

  # Special zone for Norway
  if ($latitude>=56.0 && $latitude<64.0 && $longitude>=3.0 && $longitude<12.0) {
    $longitudeZone=32;
  }

  # Special zones for Svalbard
  if ($latitude>=72.0 && $latitude<84.0) {
    if ($longitude>=0.0 && $longitude<9.0) {
      $longitudeZone=31;
    } elsif ($longitude>=9.0 && $longitude<21.0) {
      $longitudeZone=33;
    } elsif ($longitude>=21.0 && $longitude<33.0) {
      $longitudeZone=35;
    } elsif ($longitude>=33.0 && $longitude<42.0) {
      $longitudeZone=37;
    }
  }

  my $longitudeOrigin=($longitudeZone - 1) * 6 - 180 + 3;
  my $longitudeOriginRad=$longitudeOrigin * (pi() / 180.0);

  $UTMZone=getUTMLatitudeZoneLetter($latitude);
  my $ePrimeSquared=($eSquared) / (1 - $eSquared);

  my $n=$a / sqrt(1 - $eSquared * sin($latitudeRad) * sin($latitudeRad));
  my $t=tan($latitudeRad) * tan($latitudeRad);
  my $c=$ePrimeSquared * cos($latitudeRad) * cos($latitudeRad);
  my $A=cos($latitudeRad) * ($longitudeRad - $longitudeOriginRad);

  my $M=$a * ((1 - $eSquared / 4 - 3 * $eSquared * $eSquared / 64 - 5 * $eSquared * $eSquared * $eSquared / 256) * $latitudeRad
         - (3 * $eSquared / 8 + 3 * $eSquared * $eSquared / 32 + 45 * $eSquared * $eSquared * $eSquared / 1024) * sin(2 * $latitudeRad)
         + (15 * $eSquared * $eSquared / 256 + 45 * $eSquared * $eSquared * $eSquared / 1024) * sin(4 * $latitudeRad)
         - (35 * $eSquared * $eSquared * $eSquared / 3072) * sin(6 * $latitudeRad));

  $UTMEasting=$UTM_F0 * $n * ($A + (1 - $t + $c) * pow($A, 3.0) / 6 + (5 - 18 * $t + $t * $t + 72 * $c - 58 * $ePrimeSquared) * pow($A, 5.0) / 120) + 500000.0;
  $UTMNorthing=$UTM_F0 * ($M + $n * tan($latitudeRad) * ($A * $A / 2 + (5 - $t + (9 * $c) + (4 * $c * $c)) * pow($A, 4.0) / 24 + (61 - (58 * $t) + ($t * $t) + (600 * $c) - (330 * $ePrimeSquared))
        * pow($A, 6.0) / 720));

  # Adjust for the southern hemisphere
  if ($latitude < 0) {
    $UTMNorthing +=10000000.0;
  }
# ($UTMEasting, $UTMNorthing, $UTMZone, $longitudeZone); # Global values returned.
}    

