<?php

/**
 * datum shift using Helmert transformation
 * parameters for transformation WGS84 --> OSGB36 hard coded
 * possibility to implement other parameters should be evident
 *
 * height above the ellipsoid is ignored
 * latitude and longitude must be given in decimal degrees
 *
 * no transformation if abs(latitude)>89 or abs(longitude)>89
 * (lifting this restriction requires some more lines of code)
 */
function HelmertDatumShift ( $latitude , $longitude )
{
    /* return if latitude or longitude out of range */

    if ( abs($latitude) > 89. OR abs($longitude) > 89. )
    {
	$result[0] = $latitude;
	$result[1] = $longitude;
	return $result;
    }

    /* parameters for transformation WGS84 --> OSGB36 */

    $a1 = 6378137.0;
    $b1 = 6356752.3141;

    $a2 = 6377563.396;
    $b2 = 6356256.910;

    $tx = -446.448;
    $ty = +125.157;
    $tz = -542.060;
    $s0 = 1.0000204894;
    $rx = deg2rad ( -.1502/3600. );
    $ry = deg2rad ( -.2470/3600. );
    $rz = deg2rad ( -.8241/3600. );

    /* convert latitude and longitude to cartesian coordinates */

    $phi = deg2rad ( $latitude );
    $lambda = deg2rad ( $longitude );

    $cosphi = cos ( $phi );
    $sinphi = sin ( $phi );
    $coslambda = cos ( $lambda );
    $sinlambda = sin ( $lambda );

    $aq = $a1 * $a1;
    $e1 = ( $aq - $b1 * $b1 ) / $aq;

    $ny = $a1 / sqrt ( 1. - $e1*$sinphi*$sinphi );

    $x1 = $ny * $cosphi * $coslambda;
    $y1 = $ny * $cosphi * $sinlambda;
    $z1 = $ny * ( 1. - $e1 ) * $sinphi;

    /* the helmert transformation proper */

    $x2 = $tx + ( $s0 * $x1 - $rz * $y1 + $ry * $z1 );
    $y2 = $ty + ( $rz * $x1 + $s0 * $y1 - $rx * $z1 );
    $z2 = $tz + ( -$ry *$x1 + $rx * $y1 + $s0 * $z1 );

    /* convert cartesian coordinates to latitude and longitude */

    $aq = $a2 * $a2;
    $e2 = ( $aq - $b2 * $b2 ) / $aq;

    $lambda = atan2 ( $y2 , $x2 );

    $p2 = sqrt ( $x2*$x2 + $y2*$y2 );
    $phi = atan2 ( $z2 , $p2*(1.-$e2) );

    $grenze = 1. / $b2;
    $n0 = 11;

    do
    {
	$phi_alt = $phi;
	$ny = $a2 / sqrt ( 1. - $e2 * sin($phi) * sin($phi) );
	$phi = atan2 ( $z2 + $e2 * $ny * sin($phi) , $p2 );
	$n0 -= 1;
    } while ( abs ( $phi - $phi_alt ) > $grenze AND $n0 > 0 );

    $result[0] = rad2deg ( $phi );
    $result[1] = rad2deg ( $lambda );

    return $result;
}

?>