<?php
/***
 * Class SVY21
 * todo: convert SVY21 Singapore TM to Lat Long
 */
class SVY21
{
    const ELLIPSOID_SEMIMAJORAXIS = 6378137.0;
    const ELLIPSOID_ECCENTRICITY = 0.0818191908426215;
    const ELLIPSOID_FLATTENING = 0.00335281066474746;
    const PROJ_NATURALORIGINLATITUDE = 1.36666666666667;
    const PROJ_NATURALORIGINLONGITUDE = 103.833333333333;
    const PROJ_SCALEFACTOR = 1.0;
    const PROJ_FALSEEASTINGS = 28001.642;
    const PROJ_FALSENORTHINGS = 38744.572;

    function __construct()
    {
    }

    /***
     * @param $N of svy21
     * @param $E of svy21
     * @return float[]|int[] return [lat, long]
     */
    public function CnvEN2LL($N, $E)
    {
        $l = $this->SecondEccentricity();
        $a = $this->calc_e1();
        $b = $this->calc_M1($E);
        $c = $this->calc_u1($b);
        $d = $this->calc_lat1($c, $a);
        $e = $this->calc_T1($d);
        $f = $this->calc_c1($d);
        $g = $this->calc_v1($d);
        $h = $this->calc_p1($d);
        $i = ($N - self::PROJ_FALSEEASTINGS) / ($g * self::PROJ_SCALEFACTOR);
        $j = ($i * $i / 2.0) - (5.0 + 3.0 * $e + 10.0 * $f - 4.0 * $f * $f - 9.0 * pow($l, 2)) * (pow($i, 4) / 24.0) + (61.0 + 90.0 * $e + 298.0 * $f + 45.0 * $e * $e - 252.0 * pow($l, 2) - 3.0 * $f * $f) * (pow($i, 6) / 720.0);
        $j = $d - $g * tan($d) * $j / $h;
        $lat = $this->CnvRadToDeg($j);
        $k = ($i - (1.0 + 2.0 * $e + $f) * pow($i, 3) / 6.0 + (5.0 - 2.0 * $f + 28.0 * $e - 3.0 * $f * $f + 8.0 * pow($l, 2) + 24.0 * $e * $e) * pow($i, 5) / 120.0) / cos($d);
        $k = $this->CnvDegToRad(self::PROJ_NATURALORIGINLONGITUDE) + $k;
        $long = $this->CnvRadToDeg($k);
        $r1 = floor($lat);
        $r2 = floor(($lat - $r1) * 60);
        $r3 = (round(((($lat - $r1) - ($r2 / 60)) * 60 * 60) * 100) / 100);
        $r3 = $this->roundNumber($r3, 4);
        $r4 = floor($long);
        $r4 = floor(($long - $r4) * 60);
        $r5 = (round(((($long - $r4) - ($r4 / 60)) * 60 * 60) * 100) / 100);
        $r5 = $this->roundNumber($r5, 4);
        return [$lat, $long];
    }

    private function SecondEccentricity()
    {
        return (sqrt(self::ELLIPSOID_ECCENTRICITY * self::ELLIPSOID_ECCENTRICITY / (1.0 - self::ELLIPSOID_ECCENTRICITY * self::ELLIPSOID_ECCENTRICITY)));
    }

    private function roundNumber($a, $b)
    {
        return round($a * pow(10, $b)) / pow(10, $b);
    }

    private function CnvRadToDeg($a)
    {
        return (180.0 * $a / pi());
    }

    private function CnvDegToRad($a)
    {
        return (pi() * $a) / 180.0;
    }

    private function calc_M($a)
    {
        $c = self::ELLIPSOID_ECCENTRICITY;
        $d = 1.0 - (pow($c, 2) / 4.0) - (3.0 * pow($c, 4) / 64.0) - (5.0 * pow($c, 6) / 256);
        $e = (3.0 * pow($c, 2) / 8.0) + (3.0 * pow($c, 4) / 32.0) + (45.0 * pow($c, 6) / 1024.0);
        $f = (15.0 * pow($c, 4.0) / 256.0) + (45.0 * pow($c, 6) / 1024.0);
        $g = (35.0 * pow($c, 6) / 3072);
        $b = $d * $a - $e * sin(2.0 * $a) + $f * sin(4.0 * $a) - $g * sin(6 * $a);
        return $b * self::ELLIPSOID_SEMIMAJORAXIS;
    }

    private function calc_T1($a)
    {
        $b = tan($a);
        return $b * $b;
    }

    private function calc_v1($a)
    {
        $b = self::ELLIPSOID_ECCENTRICITY * sin($a);
        $b = sqrt(1.0 - $b * $b);
        return self::ELLIPSOID_SEMIMAJORAXIS / $b;
    }

    private function calc_p1($a)
    {
        $c = self::ELLIPSOID_ECCENTRICITY;
        $b = 1.0 - $c * $c * pow(sin($a), 2);
        $b = pow($b, 3.0 / 2.0);
        return self::ELLIPSOID_SEMIMAJORAXIS * (1.0 - $c * $c) / $b;
    }

    private function calc_c1($a)
    {
        $b = $this->SecondEccentricity();
        $b = $b * cos($a);
        return $b * $b;
    }

    private function calc_e1()
    {
        $a = 1.0 - sqrt(1.0 - pow(self::ELLIPSOID_ECCENTRICITY, 2));
        return $a / (1.0 + sqrt(1.0 - pow(self::ELLIPSOID_ECCENTRICITY, 2)));
    }

    private function calc_M1($a)
    {
        $c = $this->calc_M($this->CnvDegToRad(self::PROJ_NATURALORIGINLATITUDE));
        return ($c + ($a - self::PROJ_FALSENORTHINGS) / self::PROJ_SCALEFACTOR);
    }

    private function calc_u1($a)
    {
        $c = self::ELLIPSOID_ECCENTRICITY;
        $b = self::ELLIPSOID_SEMIMAJORAXIS * (1.0 - ($c * $c / 4.0) - (3.0 * pow($c, 4) / 64.0) - (5.0 * pow($c, 6) / 256.0));
        return $a / $b;
    }

    private function calc_lat1($a, $b)
    {
        $c = $a + ((3.0 * $b / 2.0) - 27.0 * $b * $b / 32.0) * sin(2.0 * $a);
        $c = $c + ((21.0 * $b * $b / 16.0) - 55.0 * pow($b, 4.0) / 32.0) * sin(4.0 * $a);
        $c = $c + (151.0 * pow($b, 3) / 96.0) * sin(6.0 * $a);
        $c = $c + (1097.0 * pow($b, 4) / 512.0) * sin(8.0 * $a);
        return $c;
    }

}