Implementation for WGS84 (lat,lon) conversion to Greek Grid (GGRS87)

In my last  post I said that I had found a site with an algorithm that would convert WGS84 (lat/lon) coordinates to the Greek Grid system (GGRS87).  I successfully converted that algorithm into a .php program and I present it here.  The lines commented out are the lines from the algorithm.  The immediately following line is my php interpretation of this.  As it is it rounds down to the whole meter.  But its accuracy appears to be of the order of one one-thousandth of a meter (millimeter).  Plenty good enough for the likes of us.


<?php

// from algorithm at: http://studentguru.gr/b/kokyri/archive/2012/02/21/convert-coordinates-from-greek-grid-ggrs87-87-espg-2100-to-latitude-longitude-wsg84

//public static double[] ToGreekGrid(double lat, double lng)
function ToGreekGrid($lat, $lng)
//        {
{
//            lat = lat * Math.PI / 180;
$lat = $lat * M_PI / 180.0;
//            lng = lng * Math.PI / 180;
$lng = $lng * M_PI / 180.0;
//            var a = 6378137;
$a = 6378137;
//            var b = 6356752.31424518;
$b = 6356752.31424518;
//            var e2 = (Math.Pow(a,2) - Math.Pow(b,2)) / Math.Pow(a, 2);
$e2 = ($a * $a - $b * $b) / ($a * $a);
//            var v = a / Math.Sqrt(1 - e2 * (Math.Sin(lat) * Math.Sin(lat)));
$v = $a / sqrt(1 - $e2 * (sin($lat) * sin($lat)));
//            var X = v * Math.Cos(lat) * Math.Cos(lng);
$X = $v * cos($lat) * cos($lng);
//            var Y = v * Math.Cos(lat) * Math.Sin(lng);
$Y = $v * cos($lat) * sin($lng);
//            var Z = (v * (1 - e2)) * Math.Sin(lat);
$Z = ($v * (1.0 - $e2)) * sin($lat);
//            var px = 199.723;
$px = 199.723;
//            var py = -74.03;
$py = -74.03;
//            var pz = -246.018;
$pz = -246.018;
//            var rw = 0;
$rw = 0;
//            var rf = 0;
$rf = 0;
//            var rk = 0;
$rk = 0;
//            var ks = 1;
$ks = 1;
//            var c1 = Math.Cos(rw);
$c1 = cos($rw);
//            var c2 = Math.Cos(rf);
$c2 = cos($rf);
//            var c3 = Math.Cos(rk);
$c3 = cos($rk);
//            var s1 = Math.Sin(rw);
$s1 = sin($rw);
//            var s2 = Math.Sin(rf);
$s2 = sin($rf);
//            var s3 = Math.Sin(rk);
$s3 = sin($rk);
//            var D11 = c2 * c3;
$D11 = $c2 * $c3;

//            var D21 = -c2 * s3;
$D21 = $c2*(-1) * $s3;

//            var D31 = Math.Sin(rf);
$D31 = sin($rf);

//            var D12 = s1 * s2 * c3 + c1 * s3;
$D12 = $s1 * $s2 * $c3 + $c1 * $s3;

//            var D22 = -s1 * s2 * s3 + c1 * c3;
$D22 = $s1*(-1) * $s2 *$s3 + $c1 * $c3;

//            var D32 = -s1 * c2;
$D32 = $s1*(-1) * $c2;

//            var D13 = -c1 * s2 * c3 + s1 * s3;
$D13 = $c1*(-1) * $s2 * $c3 + $s1 * $s3;

//            var D23 = c1 * s2 * s3 + s1 * c3;
$D23 = $c1 * $s2 * $s3 + $s1 * $c3;

//            var D33 = c1 * c2;
$D33 = $c1 * $c2;

//            var X1 = px + ks * (D11 * X + D12 * Y + D13 * Z);
$X1 = $px + $ks * ($D11 * $X + $D12 * $Y + $D13 * $Z);

//            var Y1 = py + ks * (D21 * X + D22 * Y + D23 * Z);
$Y1 = $py + $ks * ($D21 * $X + $D22 * $Y + $D23 * $Z);

//            var Z1 = pz + ks * (D31 * X + D32 * Y + D33 * Z);
$Z1 = $pz + $ks * ($D31 * $X + $D32 * $Y + $D33 * $Z);

//            lng = Math.Atan2(Y1, X1);
$lng = atan2($Y1, $X1);

//            var lat0 = Math.Atan2(Z1, Math.Sqrt(X1 * X1 + Y1 * Y1));
$lat0 = atan2($Z1, sqrt($X1 * $X1 + $Y1 * $Y1));

//            b = 6356752.31414036;
$b = 6356752.31414036;

//            e2 = (Math.Pow(a, 2) - Math.Pow(b,2)) / Math.Pow(a, 2);
$e2 = ($a * $a - $b * $b) / ($a * $a);

//            while (Math.Abs(lat - lat0) > 0.0000000001)
while (abs($lat - $lat0) > 0.0000000001)
//            {
{

//                var No = a / Math.Sqrt(1 - e2 * Math.Sin(lat0) * Math.Sin(lat0));
$No = $a / sqrt(1 - $e2 * sin($lat0) * sin($lat0));


//                var h = Math.Sqrt(X1 * X1 + Y1 * Y1) / Math.Cos(lat0) - No;
$h = sqrt($X1 * $X1 + $Y1 * $Y1) / cos($lat0) - $No;

//                lat = lat0;
$lat = $lat0;

//                lat0 = Math.Atan(Z1 / Math.Sqrt(X1 * X1 + Y1 * Y1) * (1 / (1 - e2 * No / (No+ h))));
$lat0 = atan($Z1 /sqrt($X1 * $X1 + $Y1 * $Y1) * (1 / (1 - $e2 * $No / ($No + $h))));
//            }
}
//            lng = lng - 24 * Math.PI / 180;
$lng = $lng - 24 * M_PI / 180;

//            var m0 = 0.9996;
$m0 = 0.9996;

//            var es2 = (Math.Pow(a,2) - Math.Pow(b,2)) / (Math.Pow(b,2));
$es2 = ($a * $a - $b * $b) / ($b * $b);

//            var V = Math.Sqrt(1 + es2 * Math.Cos(lat) * Math.Cos(lat));
$V = sqrt(1 + $es2 * cos($lat) * cos($lat));

//            var eta = Math.Sqrt(es2 * Math.Cos(lat) * Math.Cos(lat));
$eta = sqrt($es2 * cos($lat) * cos($lat));

//            var Bf = Math.Atan(Math.Tan(lat) / Math.Cos(V * lng) * (1 + eta * eta / 6 * (1 - 3 * Math.Sin(lat) * Math.Sin(lat)) * lng * lng * lng * lng));
                     $Bf =  atan(tan($lat) / cos($V * $lng) * (1 + $eta * $eta / 6 * (1 - 3 * sin($lat) * sin($lat)) * $lng * $lng * $lng * $lng));

//            var Vf = Math.Sqrt(1 + es2 * Math.Cos(Bf) * Math.Cos(Bf));
$Vf = sqrt(1 + $es2 * cos($Bf) * cos($Bf));

//            var etaf = Math.Sqrt(es2 * Math.Cos(Bf) * Math.Cos(Bf));
$etaf = sqrt($es2 * cos($Bf) * cos($Bf));

//            var n = (a - b) / (a + b);
$n = ($a - $b) / ($a + $b);

//            var r1 = (1 + n * n / 4 + n * n * n * n / 64) * Bf;
$r1 = (1 + $n * $n / 4 + $n * $n * $n * $n / 64) * $Bf;

//            var r2 = 3.0 / 2.0 * n * (1 - n * n / 8) * Math.Sin(2 * Bf);
$r2 = 3.0 / 2.0 * $n * (1 - $n * $n / 8) * sin(2 * $Bf);

//            var r3 = 15.0 / 16.0 * n * n * (1 - n * n / 4) * Math.Sin(4 * Bf);
$r3 = 15.0 / 16.0 * $n * $n * (1 - $n * $n / 4) * sin (4 * $Bf);

//            var r4 = 35.0 / 48.0 * n * n * n * Math.Sin(6 * Bf);
$r4 = 35.0 / 48.0 * $n * $n * $n * sin(6 * $Bf);

//            var r5 = 315.0 / 512.0 * n * n * n * n * Math.Sin(8 * Bf);
$r5 = 315.0 / 512.0 * $n * $n * $n * $n * sin(8 * $Bf);

//            var Northing = a / (1 + n) * (r1 - r2 + r3 - r4 + r5) * m0 - 0.001 + 4202812 - 4207988.1206046063;
$Northing = $a / (1 + $n) * ($r1 - $r2 + $r3 - $r4 + $r5) * $m0 - 0.001 + 4202812 - 4207988.1206046063;

//            var ys = Math.Tan(lng) * Math.Cos(Bf) / Vf * (1 + etaf * etaf * lng * lng * Math.Cos(Bf) * Math.Cos(Bf) * (etaf * etaf / 6 + lng * lng / 10));

$ys = tan($lng) * cos($Bf) / $Vf * (1 + $etaf * $etaf * $lng * $lng * cos($Bf) * cos($Bf) * ($etaf * $etaf / 6 + $lng * $lng / 10));

//            ys = Math.Log(ys + Math.Sqrt(ys * ys + 1));
$ys = log($ys + sqrt($ys * $ys + 1));

//            var Easting = m0 * Math.Pow(a,2) / b * ys + 500000;
$Easting = $m0 * pow($a, 2) / $b * $ys + 500000;

$Northing = $Northing + 5176.122; // Added correction per comments

// echo $Easting.PHP_EOL;
// echo $Northing.PHP_EOL.PHP_EOL;

return array(0 => round($Easting), 1 => round($Northing));

//        }
}

/// main

//ToGreekGrid (37.97412198000, 22.59898503000);

list ($Easting, $Northing) = ToGreekGrid(37.97412198000, 22.59898503000);

echo PHP_EOL.PHP_EOL;
echo $Easting.PHP_EOL;
echo $Northing;

?>


Commentary

You'll need a php interpreter to run this, of course.  If you want to copy this stand-alone program from Google Drive then click on this.

The main routine does nothing but call the routine ToGreekGrid and print out the Easting and Northing results.  You'll modify that, of  course.  I see no barrier now to my adding the GGRS87 coordinates to every entry  in the Mycenaean Atlas Project.

If you try this then let me know what you think.