[PHP] Convert EPSG2100 coordinates to WGS84 coordinates

functions.php

function transformEPSG2100toWGS84($x, $y)
{
    // EPSG:2100 parameters
    $a = 6378137.0; // Semi-major axis
    $f = 1 / 298.257222101; // Flattening
    $e2 = 2 * $f - $f ** 2; // Square of eccentricity
    $k0 = 0.9996; // Scale factor
    $lambda0 = deg2rad(24); // Central meridian in radians
    $E0 = 500000; // False easting
    $N0 = 0; // False northing

    // Remove false easting and northing
    $E = $x - $E0;
    $N = $y - $N0;

    // Calculate meridional arc
    $M = $N / $k0;

    // Calculate the latitude using the iterative method
    $mu = $M / ($a * (1 - $e2 / 4 - 3 * $e2 ** 2 / 64 - 5 * $e2 ** 3 / 256));
    $e1 = (1 - sqrt(1 - $e2)) / (1 + sqrt(1 - $e2));

    $phi1 = $mu +
        (3 * $e1 / 2 - 27 * $e1 ** 3 / 32) * sin(2 * $mu) +
        (21 * $e1 ** 2 / 16 - 55 * $e1 ** 4 / 32) * sin(4 * $mu) +
        (151 * $e1 ** 3 / 96) * sin(6 * $mu) +
        (1097 * $e1 ** 4 / 512) * sin(8 * $mu);

    // Calculate the transverse radius of curvature
    $N1 = $a / sqrt(1 - $e2 * sin($phi1) ** 2);
    $T1 = tan($phi1) ** 2;
    $C1 = $e2 / (1 - $e2) * cos($phi1) ** 2;
    $R1 = $a * (1 - $e2) / (1 - $e2 * sin($phi1) ** 2) ** 1.5;
    $D = $E / ($N1 * $k0);

    // Calculate latitude and longitude
    $lat = $phi1 - ($N1 * tan($phi1) / $R1) *
        ($D ** 2 / 2 -
            (5 + 3 * $T1 + 10 * $C1 - 4 * $C1 ** 2 - 9 * $e2) * $D ** 4 / 24 +
            (61 + 90 * $T1 + 298 * $C1 + 45 * $T1 ** 2 - 252 * $e2 - 3 * $C1 ** 2) * $D ** 6 / 720);

    $lon = $lambda0 + ($D -
        (1 + 2 * $T1 + $C1) * $D ** 3 / 6 +
        (5 - 2 * $C1 + 28 * $T1 - 3 * $C1 ** 2 + 8 * $e2 + 24 * $T1 ** 2) * $D ** 5 / 120) / cos($phi1);

    // Convert radians to degrees
    $lat = rad2deg($lat);
    $lon = rad2deg($lon);

    // Apply the full Helmert transformation (TOWGS84 parameters)
    $helmertParams = [
        'dx' => -199.87, // Translation in X (meters)
        'dy' => 74.79,   // Translation in Y (meters)
        'dz' => 246.62,  // Translation in Z (meters)
        'rx' => 0.0,     // Rotation about X-axis (arc-seconds)
        'ry' => 0.0,     // Rotation about Y-axis (arc-seconds)
        'rz' => 0.0,     // Rotation about Z-axis (arc-seconds)
        'scale' => 0.0   // Scale factor (ppm)
    ];

    // Convert lat/lon to Cartesian coordinates
    $phi = deg2rad($lat);
    $lambda = deg2rad($lon);
    $cartesian = geodeticToCartesian($phi, $lambda, $a, $f);

    // Apply Helmert transformation
    $transformedCartesian = helmertTransform($cartesian, $helmertParams);

    // Convert back to geodetic coordinates (lat/lon)
    $result = cartesianToGeodetic($transformedCartesian, $a, $f);

    return ['lat' => rad2deg($result['lat']), 'lon' => rad2deg($result['lon'])];
}

function geodeticToCartesian($phi, $lambda, $a, $f)
{
    $e2 = 2 * $f - $f ** 2; // Square of eccentricity
    $N = $a / sqrt(1 - $e2 * sin($phi) ** 2); // Radius of curvature in the prime vertical
    $x = $N * cos($phi) * cos($lambda);
    $y = $N * cos($phi) * sin($lambda);
    $z = $N * (1 - $e2) * sin($phi);
    return ['x' => $x, 'y' => $y, 'z' => $z];
}

function helmertTransform($cartesian, $params)
{
    $x = $cartesian['x'] + $params['dx'];
    $y = $cartesian['y'] + $params['dy'];
    $z = $cartesian['z'] + $params['dz'];

    // Apply rotations and scale if needed
    return ['x' => $x, 'y' => $y, 'z' => $z];
}

function cartesianToGeodetic($cartesian, $a, $f)
{
    $x = $cartesian['x'];
    $y = $cartesian['y'];
    $z = $cartesian['z'];

    $e2 = 2 * $f - $f ** 2; // Square of eccentricity
    $p = sqrt($x ** 2 + $y ** 2);
    $phi = atan2($z, $p * (1 - $e2));
    $lambda = atan2($y, $x);

    // Iteratively refine latitude
    $N = $a / sqrt(1 - $e2 * sin($phi) ** 2);
    $phiPrev = 0;

    while (abs($phi - $phiPrev) > 1e-12) {
        $phiPrev = $phi;
        $phi = atan2($z + $e2 * $N * sin($phi), $p);
        $N = $a / sqrt(1 - $e2 * sin($phi) ** 2);
    }

    return ['lat' => $phi, 'lon' => $lambda];
}