<?php
// Lat Long functions in PHP
//
// found by RWH at http://www.megalithia.com/search/llfuncshighlight.php
//
// With thanks to Andy, G4JNT for inspiration in GEOG, and to the OSGB for their white paper on coordinate transformation
// describing the iterative method used
// thanks to the Ordnance survey of Ireland for details of the true and false origins of the Irish grid
//
// You may use and redistribute this code under the terms of the GPL see http://www.gnu.org/copyleft/gpl.html
//
//
// Written by Richard
// www.megalithia.com
// v0.something 27/2/2000
// v 1.01 28 June 2004
// v 1.02 6 Aug 2004 line 89 add "0" to chars in $ngr=stripcharsnotinbag Thx Andy
// v 1.03 9 Mar 2005 Jan (Klingon) added conversion to WGS84 map datum and removed limitation of digits of the grid ref
// v 1.04 10 Aug 2005 Richard correct error trapping (only manifest on malformed ngrs
//
// This code is predicated on the assumption that your are ONLY feeding it Irish or UK Grid references.
// It uses the single char prefix of Irish grid refs to tell the difference, UK grid refs have a two letter prefix.
// We would like an even number of digits for the rest of the grid ref.
// Anything in the NGR other than 0-9, A-Z, a-z is eliminated.
// WARNING this assumes there are no decimal points in your NGR components.
// (before v 1.03) you must have at least 4 numerical digits eg TM1234 and no more than eight eg TM12345678 otherwise
// you and your data are viewed with suspicion.
// NEW (Jan): at least the letter prefix of the Grid ref, no upper limit
//
// The transformation from OSGB36/Ireland 1965 to WGS84 is more precise than 5 m.
//
// The function is case insensitive
//
//call_user_func("NGR2LL",$ngr);
// Return value is array($country,$error,$lat,$long);
// you may test error = OK if blank, otherwise you will usually get a "malformed NGR"

// Removes all characters which do NOT appear in string $bag
// from string $s.
function stripCharsNotInBag($s,$bag)
{
$returnString="";
// Search through string's characters one by one.
// If character is in bag, append to returnString.
for($i=0;$i<strlen($s);$i++) {
      // Check that current character isn't whitespace.
      $c=substr($s,$i,1);
      if(strstr($bag,$c)) { $returnString.=$c; }
}
return $returnString;
}

// haversine distance computation.
// Haversines are good for small diffs in lat,long which is usually the case in terrestrial GPS calcs
// you can KO this function if not wanted - it is not part of the NGR2LL calc
function GCdist($lat1,$lon1,$lat2,$lon2)
{
// expects input in radians, give result in km
// mean radius of the earth in km
$earth_rad=6366.71;
// haversinehttp://mathforum.org/library/drmath/view/51879.html
$dlon=$lon2-$lon1;
$dlat=$lat2-$lat1;
// above all in degrees mind you
$a=pow(sin($dlat/2.0),2.0)+cos($lat1)*cos($lat2)*pow(sin($dlon/2.0),2.0);
$c=2*atan2(sqrt($a),sqrt(1.0-$a));
$d=$earth_rad*$c;
return $d;
// old non-haversine return acos(sin($lat1)*sin($lat2)+cos($lat1)*cos($lat2)*cos($lon1-$lon2))*$earth_rad;
}

function bearing($lat1,$lon1,$lat2,$lon2)
{
// Inspiration thanks to http://www.movable-type.co.uk/scripts/LatLong.html
$dlong=($lon2-$lon1)*M_PI/180.0;
$lat1=$lat1*M_PI/180.0;
$lat2=$lat2*M_PI/180.0;
return( 180/M_PI*atan2(sin($dlong)*cos($lat2),cos($lat1)*sin($lat2) - sin($lat1)*cos($lat2)*cos($dlong)) );
}

// Inner function of NGR2LL() for GB accepting eastings and northings in metres
function GBgr_num2ll($e,$n) 
{
     // deg2rad
     $dr=180.0/M_PI;
     // True Origin is 2 deg W
     $phi0uk=-2.0;
     // True Origin is 49 deg N
     $lambda0uk=49.0;
     // scale factor @ central meridian
     $F0uk=0.9996012717;
     // True origin in 400 km E of false origin
     $E0uk=400000.0;
     //True origin is 100 km S of false origin
     $N0uk=-100000.0;
     // semi-major axis (in line to equator) 0.996012717 is yer scale @ central meridian
     $auk=6377563.396*$F0uk;
     //semi-minor axis (in line to poles)
     $buk=6356256.91*$F0uk;
     // flatness=a1-b1/(a1+b1)
     $n1uk=0.00167322025032508731869331280635710896296;
     // first eccentricity squared=2*f-f^2where f=(a1-b1)/a1
     $e2uk=0.006670539761597529073698869358812557054558;
     // radius of earth
     //$re=6371.29;
     $k=($n-$N0uk)/$auk+$lambda0uk/$dr;
     $nextcounter=0;
     do {
         $nextcounter=$nextcounter+1;
         $k3=$k-$lambda0uk/$dr;
         $k4=$k+$lambda0uk/$dr;
         $j3=(1.0+$n1uk+1.25*pow($n1uk,2.0)+1.25*pow($n1uk,3.0))*$k3;
         $j4=(3.0*$n1uk+3.0*pow($n1uk,2.0)+2.625*pow($n1uk,3.0))*sin($k3)*cos($k4);
         $j5=(1.875*pow($n1uk,2.0)+1.875*pow($n1uk,3.0))*sin(2.0*$k3)*cos(2.0*$k4);
         $j6=35.0/24.0*pow($n1uk,3.0)*sin(3.0*$k3)*cos(3.0*$k4);
         $m=$buk*($j3-$j4+$j5-$j6);
         $k=$k+($n-$N0uk-$m)/$auk;
     } while((abs($n-$N0uk-$m)>0.000000001)&&($nextcounter<100));
   # echo "loop 1 ". abs($n-$N0uk-$m) ." $nextcounter llll<br/>";
     $v=$auk/sqrt(1.0-$e2uk*pow(sin($k),2.0));
     $r=$v*(1.0-$e2uk)/(1.0-$e2uk*pow(sin($k),2.0));
     $h2=$v/$r-1.0;
     $y1=$e-$E0uk;
     $j3=tan($k)/(2.0*$r*$v);
     $j4=tan($k)/(24.0*$r*pow($v,3.0))*(5.0+3.0*pow(tan($k),2.0)+$h2-9.0*pow(tan($k),2.0)*$h2);
     $j5=tan($k)/(720.0*$r*pow($v,5.0))*(61.0+90.0*pow(tan($k),2.0)+45.0*pow(tan($k),4.0));
     $k9=$k-$y1*$y1*$j3+pow($y1,4.0)*$j4-pow($y1,6.0)*$j5;
     $j6=1.0/(cos($k)*$v);
     $j7=1.0/(cos($k)*6.0*pow($v,3.0))*($v/$r+2.0*pow(tan($k),2.0));
     $j8=1.0/(cos($k)*120.0*pow($v,5.0))*(5.0+28.0*pow(tan($k),2.0)+24.0*pow(tan($k),4.0));
     $j9=1.0/(cos($k)*5040.0*pow($v,7.0));
     $j9=$j9*(61.0+662.0*pow(tan($k),2.0)+1320.0*pow(tan($k),4.0)+720.0*pow(tan($k),6.0));
     $long=$phi0uk+$dr*($y1*$j6-$y1*$y1*$y1*$j7+pow($y1,5.0)*$j8-pow($y1,7.0)*$j9);
     $lat=$k9*$dr;
     // v1.04 this bracket moved to just before elsif // }
    // convert long/lat to Cartesian coordinates
    $v=6377563.396/sqrt(1.0-$e2uk*pow(sin($k),2.0));
    $cartxa=$v*cos($k9)*cos($long/$dr);
    $cartya=$v*cos($k9)*sin($long/$dr);
    $cartza=(1.0-$e2uk)*$v*sin($k9);
    // Helmert-Transformation from OSGB36 to WGS84 map date
    $rotx=-0.1502/3600.0*M_PI/180.0;
    $roty=-0.2470/3600.0*M_PI/180.0;
    $rotz=-0.8421/3600.0*M_PI/180.0;
    $scale=-20.4894/1000000.0;
    $cartxb=446.448+(1.0+$scale)*$cartxa+$rotz*$cartya-$roty*$cartza;
    $cartyb=-125.157-$rotz*$cartxa+(1.0+$scale)*$cartya+$rotx*$cartza;
    $cartzb=542.06+$roty*$cartxa-$rotx*$cartya+(1.0+$scale)*$cartza;
    // convert Cartesian to long/lat
    $awgs84=6378137.0;
    $bwgs84=6356752.3141;
    $e2wgs84=0.00669438003551279089034150031998869922791;
    $lambdaradwgs84=atan($cartyb/$cartxb);
    $long=$lambdaradwgs84*180.0/M_PI;
    $pxy=sqrt(pow($cartxb,2.0)+pow($cartyb,2.0));
    $phiradwgs84=atan($cartzb/$pxy/(1.0-$e2wgs84));
    $nextcounter=0;
    do {
        $nextcounter=$nextcounter+1;
        $v=$awgs84/sqrt(1.0-$e2wgs84*pow(sin($phiradwgs84),2.0));
        $phinewwgs84=atan(($cartzb+$e2wgs84*$v*sin($phiradwgs84))/$pxy);
        $phiradwgs84=$phinewwgs84;
    } while((abs($phinewwgs84-$phiradwgs84)>0.000000000001)&&($nextcounter<100));
  # echo "l2l2 $nextcounter l2l2 ";
    $lat=$phiradwgs84*180.0/M_PI;
    // } v 1.04 mod
    return array($lat,$long);
}

function NGR2LL($ngr)
{
// returns a error,lat,long list
$country="";
$ngr=stripCharsNotInBag(strtoupper($ngr),"0123456789ABCDEFGHJKLMNOPQRSTUVWXYZ");
$lett=stripCharsNotInBag($ngr,"ABCDEFGHJKLMNOPQRSTUVWXYZ");

if(strlen($lett)==1)
{
$error="";
$lat=0.0;
$long=0.0;
$num=stripCharsNotInBag($ngr,"01232456789");
$le=strlen($num);
$country="Irish";
if($le%2==1) {
    // bust odd numerical parts
    $error="Malformed numerical part of NGR";
    $lat=0.0;
    $long=0.0;
} else {
    $pr=$le/2;
    // split into northings
    $n=substr($num,$pr);
    // and eastings
    $e=substr($num,0,$pr);
    $pr=pow(10.0,(5.0-$pr));
    $T1=ord(substr($lett,0,1))-65;
    if($T1>8)
    {
        $T1=$T1-1;
    }
    $e=100000.0*($T1%5.0)+$e*$pr;
    $n=$n*$pr+100000.0*(4.0-floor($T1/5.0));
    // USE IRISH values
    // okay up to here
    // now for the heavy stuff
    // deg2rad
    $dr=180.0/M_PI;
    // True Origin is 8 deg W
    $phi0ir=-8.0;
    // True Origin is 53.5 deg N
    $lambda0ir=53.5;
    // scale factor @ central meridian
    $F0ir=1.000035;
    // True origin in 200 km E of false origin
    $E0ir=200000.0;
    //True origin is 250km N of false origin
    $N0ir=250000.0;
    // semi-major axis (in line to equator) 1.000035 is yer scale @ central meridian
    $air=6377340.189*$F0ir;
    //semi-minor axis (in line to poles)
    $bir=6356034.447*$F0ir;
    // flatness=a1-b1/(a1 + b1)
    $n1ir=0.001673220384152058651484728058385228837777;
    // first eccentricity squared=2*f-f^2 where f=(a1-b1)/a1
    $e2ir=0.006670540293336110419293763349975612794125;
    // radius of earth
    //$re=6371.29;
    $k=($n-$N0ir)/$air+$lambda0ir/$dr;
    $nextcounter=0;
    do
    {
    $nextcounter=$nextcounter+1;
    $k3=$k-$lambda0ir/$dr;
    $k4=$k+$lambda0ir/$dr;
    $j3=(1.0+$n1ir+1.25*pow($n1ir,2.0)+1.25*pow($n1ir,3.0))*$k3;
    $j4=(3.0*$n1ir+3.0*pow($n1ir,2.0)+2.625*pow($n1ir,3.0))*sin($k3)*cos($k4);
    $j5=(1.875*pow($n1ir,2.0)+1.875*pow($n1ir,3.0))*sin(2.0*$k3)*cos(2.0*$k4);
    $j6=35.0/24.0*pow($n1ir,3.0)*sin(3.0*$k3)*cos(3.0*$k4);
    $m=$bir*($j3-$j4+$j5-$j6);
    $k=$k+($n-$N0ir-$m)/$air;
    }
    // Loop
    while((abs($n-$N0ir-$m)>0.000000000001)&&($nextcounter<10000));
    $v=$air/sqrt(1.0-$e2ir*pow(sin($k),2.0));
    $r=$v*(1.0-$e2ir)/(1.0-$e2ir*pow(sin($k),2.0));
    $h2=$v/$r-1.0;
    $y1=$e-$E0ir;
    $j3=tan($k)/(2.0*$r*$v);
    $j4=tan($k)/(24.0*$r*pow($v,3.0))*(5.0+3.0*pow(tan($k),2.0)+$h2-9.0*pow(tan($k),2.0)*$h2);
    $j5=tan($k)/(720.0*$r*pow($v,5.0))*(61.0+90.0*pow(tan($k),2.0)+45.0*pow(tan($k),4.0));
    $k9=$k-$y1*$y1*$j3+pow($y1,4.0)*$j4-pow($y1,6.0)*$j5;
    $j6=1.0/(cos($k)*$v);
    $j7=1.0/(cos($k)*6.0*pow($v,3.0))*($v/$r+2.0*pow(tan($k),2.0));
    $j8=1.0/(cos($k)*120.0*pow($v,5.0))*(5.0+28.0*pow(tan($k),2.0)+24.0*pow(tan($k),4.0));
    $j9=1.0/(cos($k)*5040.0*pow($v,7.0));
    $j9=$j9*(61.0+662.0*pow(tan($k),2.0)+1320.0*pow(tan($k),4.0)+720.0*pow(tan($k),6.0));
    $long=$phi0ir+$dr*($y1*$j6-$y1*$y1*$y1*$j7+pow($y1,5.0)*$j8-pow($y1,7.0)*$j9);
    $lat=$k9*$dr;

    // v1.04 this bracket moved to just before elsif // }
    // convert long/lat to Cartesian coordinates
    $v=6377340.189/sqrt(1.0-$e2ir*pow(sin($k),2.0));
    $cartxa=$v*cos($k9)*cos($long/$dr);
    $cartya=$v*cos($k9)*sin($long/$dr);
    $cartza=(1.0-$e2ir)*$v*sin($k9);
    // Helmert-Transformation from Ireland 1965 to WGS84 map date
    $rotx=1.042/3600.0*M_PI/180.0;
    $roty=0.214/3600.0*M_PI/180.0;
    $rotz=0.631/3600.0*M_PI/180.0;
    $scale=8.15/1000000.0;
    $cartxb=482.53+(1.0+$scale)*$cartxa+$rotz*$cartya-$roty*$cartza;
    $cartyb=-130.596-$rotz*$cartxa+(1.0+$scale)*$cartya+$rotx*$cartza;
    $cartzb=564.557+$roty*$cartxa-$rotx*$cartya+(1.0+$scale)*$cartza;
    // convert Cartesian to long/lat
    $awgs84=6378137.0;
    $bwgs84=6356752.3141;
    $e2wgs84=0.00669438003551279089034150031998869922791;
    $lambdaradwgs84=atan($cartyb/$cartxb);
    $long=$lambdaradwgs84*180.0/M_PI;
    $pxy=sqrt(pow($cartxb,2.0)+pow($cartyb,2.0));
    $phiradwgs84=atan($cartzb/$pxy/(1.0-$e2wgs84));
    $nextcounter=0;
    do
    {
    $nextcounter=$nextcounter+1;
    $v=$awgs84/sqrt(1.0-$e2wgs84*pow(sin($phiradwgs84),2.0));
    $phinewwgs84=atan(($cartzb+$e2wgs84*$v*sin($phiradwgs84))/$pxy);
    $phiradwgs84=$phinewwgs84;
    }
    // Loop
    while((abs($phinewwgs84-$phiradwgs84)>0.000000000001)&&($nextcounter<10000));
    $lat=$phiradwgs84*180.0/M_PI;
    }
        } // v 1.04 mod
elseif(strlen($lett)==2)
{
    // British
    // first caclulate e,n
    $country="UK";
    $num=stripCharsNotInBag($ngr,"0123456789");
    $le=strlen($num);

    if($le%2==1) {
    // bust odd numerical parts
        $error="Malformed numerical part of NGR";
        $lat=0.0;
        $long=0.0;
    } else {
        $pr=$le/2;
        // split into northings
        $n=substr($num,$pr);
        // and eastings
        $e=substr($num,0,$pr);
        $pr=pow(10.0,(5.0-$pr));
        $T1=ord(substr($lett,0,1))-65;
        if($T1>8)
            {
            $T1=$T1-1;
            }
        $T2=ord(substr($lett,1,1))-65;
        if($T2>8)
            {
            $T2=$T2-1;
            }
        $e=500000.0*($T1%5.0)+100000.0*($T2%5.0)-1000000.0+$e*$pr;
        $n=1900000.0-500000.0*floor($T1/5.0)-100000.0*floor($T2/5.0)+$n*$pr;
        // USE BRitish values
        // okay up to here
        // now for the heavy stuff
        list($lat,$long)=GBgr_num2ll($e,$n);
    }
}
else
{
    $error="Malformed NGR";
    $lat=0.0;
    $long=0.0;
}
return array($country,$error,$lat,$long);
}
?>