se/SE/se-lib/EpsgConversion.php

<?php

class EpsgConversion {

	private static $fe = 500000;

	private static function CalculateESquared($a, $b) {
		return (pow($a, 2) - pow($b, 2)) / pow($a, 2);
	}

	private static function CalculateE2Squared($a, $b) {
		return (pow($a, 2) - pow($b, 2)) / pow($b, 2);
	}

	private static function denom ($es, $sphi) {
		$sinSphi = sin($sphi);
		return sqrt(1 - $es * pow($sinSphi, 2));
	}

	private static function sphsr ($a, $es, $sphi) {
		$dn = self::denom($es, $sphi);
		return $a * (1 - $es) / pow($dn, 3);
	}

	private static function sphsn ($a, $es, $sphi) {
		$sinSphi = sin($sphi);
		return $a / sqrt(1 - $es * pow($sinSphi, 2));
	}

	private static function sphtmd ($ap, $bp, $cp, $dp, $ep, $sphi) {
		return ($ap * $sphi) - ($bp * sin(2 * $sphi)) + ($cp * sin(4 * $sphi)) - ($dp * sin(6 * $sphi)) + ($ep * sin(8 * $sphi));
	}

	private static function LonLatToPUWG($a, $f, $lon, $lat, $proj) {
		if ($lon < 13.5 || $lon >= 25.5) throw new Exception("Wrong longitude - {$lon}");

		if ($proj == 1) {
			$alam = 19 * M_PI / 180;
			$strf = 0;
			$nfn = -5300000.0;
			$ok = 0.9993;
		} else {
			$nfn = 0;
			$olam = floor(($lon + 1.5)/3) * M_PI / 60;
			$strf = floor(($lon + 1.5)/3) * 1000000;
			$ok = 0.999923;
		}
		
		$result = new StdClass();
		$latRad = $lat * M_PI / 180;
		$lonRad = $lon * M_PI / 180;
		$recf = 1 / $f;
		$b = $a * ($recf - 1) / $recf;
		$eSquared = self::CalculateESquared($a, $b);
		$e2Squared = self::CalculateE2Squared ($a, $b);
		$tn = ($a - $b) / ($a + $b);
		$ap = $a * (1 - $tn + 5 * (pow($tn, 2) - pow($tn, 3)) / 4 + 81 * (pow($tn, 4) - pow($tn, 5)) / 64);
		$bp = 3 * $a * ($tn - pow($tn, 2) + 7 * (pow($tn, 3) - pow($tn, 4)) / 8 + 55 * pow($tn, 5) / 64) / 2;
		$cp = 15 * $a * (pow($tn, 2) - pow($tn, 3) + 3 * (pow($tn, 4) - pow($tn, 5)) / 4) / 16;
		$dp = 35 * $a * (pow($tn, 3) - pow($tn, 4) + 11 * pow($tn, 5) / 16) / 48;
		$ep = 315 * $a * (pow($tn, 4) - pow($tn, 5)) / 512;
		$dlam = $lonRad - $olam;
		$s = sin($latRad);
		$c = cos($latRad);
		$t = $s / $c;
		$eta = $e2Squared * pow($c, 2);
		$sn = self::sphsn($a, $eSquared, $latRad);
		$tmd = self::sphtmd($ap, $bp, $cp, $dp, $ep, $latRad);
		$t1 = $tmd * $ok;
		$t2 = $sn * $s * $c * $ok / 2;
		$t3 = $sn * $s * pow($c, 3) * $ok * (5 - pow($t, 2) + 9 * $eta + 4 * pow($eta, 2)) / 24;
		$t4 = $sn * $s * pow($c, 5) * $ok * (61 - 58 * pow($t, 2) + pow($t, 4) + 270 * $eta - 330 * pow($t, 2) * $eta + 445 * pow($eta, 2) + 324 * pow($eta, 3) - 680 * pow($t, 2) * pow($eta, 2) + 88 * pow($eta, 4) - 600 * pow($t, 2) * pow($eta, 3) - 192 * pow($t, 2) * pow($eta, 4)) / 720;
		$t5 = $sn * $s * pow($c, 7) * $ok * (1385 - 3111 * pow($t, 2) + 543 * pow($t, 4) - pow($t, 6)) / 40320;
		$t6 = $sn * $c * $ok;
		$t7 = $sn * pow($c, 3) * $ok * (1 - pow($t, 2) + $eta) / 6;
		$t8 = $sn * pow($c, 5) * $ok * (5 - 18 * pow($t, 2) + pow($t, 4) + 14 * $eta - 58 * pow($t, 2) * $eta + 13 * pow($eta, 2) + 4 * pow($eta, 3) - 64 * pow($t, 2) * pow($eta, 2) - 24 * pow($t, 2) * pow($eta, 3)) / 120;
		$t9 = $sn * pow($c, 7) * $ok * (61 - 479 * pow($t, 2) + 179 * pow($t, 4) - pow($t, 6)) / 5040;
		$result->x = self::$fe + $strf + $dlam * $t6 + pow($dlam, 3) * $t7  + pow($dlam, 5) * $t8 + pow($dlam, 7) * $t9;
		$result->y = $nfn + $t1 + pow($dlam, 2) * $t2 + pow($dlam, 4) * $t3 + pow($dlam, 6) * $t4 + pow($dlam, 8) * $t5;
		return $result;
	}

	private static function PUWGToLonLat ($a, $f, $lon, $lat, $proj) {
		$result = new StdClass();
		$ok =  0.999923;
		if ($proj == 1) $ok = 0.9993;
		$recf = 1 / $f;
		$b = $a * ($recf - 1) / $recf;
		$eSquared = self::CalculateESquared($a, $b);
		$e2Squared = self::CalculateE2Squared($a, $b);
		$tn = ($a - $b) / ($a + $b);
		$ap = $a * (1 - $tn + 5 * (pow($tn, 2) - pow($tn, 3)) / 4 + 81 * (pow($tn, 4) - pow($tn, 5)) / 64);
		$bp = 3 * $a * ($tn - pow($tn, 2) + 7 * (pow($tn, 3) - pow($tn, 4)) / 8 + 55 * pow($tn, 5) / 64) / 2;
		$cp = 15 * $a * (pow($tn, 2) - pow($tn, 3) + 3 * (pow($tn, 4) - pow($tn, 5)) / 4) / 16;
		$dp = 35 * $a * (pow($tn, 3) - pow($tn, 4) + 11 * pow($tn, 5) / 16) / 48;
		$ep = 315 * $a * (pow($tn, 4) - pow($tn, 5)) / 512;
		if ($proj == 1) {
			$olam = 19 * M_PI / 180;
			$strf = 0;
			$nfn = -5300000;
          
		} else {
			$nfn = 0;
			$strf = floor($lon / 1000000) * 1000000;
			$olam = $strf / 60000000 * M_PI;
		}
		$tmd = ($lat - $nfn) / $ok;
		$sr = self::sphsr($a, $eSquared, 0);
		$ftphi = $tmd / $sr;
		for ($i = 0; $i < 5; $i++) {
			$t10 = self::sphtmd($ap, $bp, $cp, $dp, $ep, $ftphi);
			$sr = self::sphsr($a, $eSquared, $ftphi);
			$ftphi = $ftphi + ($tmd - $t10) / $sr;
		}
		$sr = self::sphsr($a, $eSquared, $ftphi);
		$sn = self::sphsn($a, $eSquared, $ftphi);
		$s = sin($ftphi);
		$c = cos($ftphi);
		$t = $s / $c;
		$eta = $e2Squared * pow($c, 2);
		$de = $lon - self::$fe - $strf;
		$t10 = $t / (2 * $sr * $sn * pow($ok, 2));
		$t11 = $t * (5 + 3 * pow($t, 2) + $eta - 4 * pow($eta, 2) - 9 * pow($t, 2) * $eta) / (24 * $sr * pow($sn, 3) * pow($ok, 4));
        
		$t12 = $t * (61 + 90 * pow($t, 2) + 46 * $eta + 45 * pow($t, 4) - 252 * pow($t, 2) * $eta - 3 * pow($eta, 2) + 100 * pow($eta, 3) - 66 * pow($t, 2) * pow($eta, 2) - 90 * pow($t, 4) * $eta + 88 * pow($eta, 4) + 225 * pow($t, 4) * pow($eta, 2) + 84 * pow($t, 2) * pow($eta, 3) - 192 * pow($t, 2) * pow($eta, 4)) / (720 * $sr * pow($sn, 5) * pow($ok, 6));
		$t13 = $t * (1385 + 3633 * pow($t, 2) + 4095 * pow($t, 4) + 1575  * pow($t, 6)) / (40320 * $sr * pow($sn, 7) * pow($ok, 8));
		$t14 = 1 / ($sn * $c * $ok);
		$t15 = (1 + 2 * pow($t, 2) + $eta) / (6 * pow($sn, 3) * $c * pow($ok, 3));
		$t16 = 1 * (5 + 6 * $eta + 28 * pow($t, 2) - 3 * pow($eta, 2) + 8 * pow($t, 2) * $eta + 24 * pow($t, 4) - 4 * pow($eta, 3) + 4 * pow($t, 2) * pow($eta, 2) + 24 * pow($t, 2) * pow($eta, 3)) / (120 * pow($sn, 5) * $c * pow($ok, 5));
		$t17 = 1 * (61 + 662 * pow($t, 2) + 1320 * pow($t, 4) + 720 * pow($t, 6)) / (5040 * pow($sn, 7) * $c * pow($ok, 7));
		$dlam = $de * $t14 - pow($de, 3) * $t15 + pow($de, 5) * $t16 - pow($de, 7) * $t17;
		$result->x = ($olam + $dlam) * 180 / M_PI;
		$result->y = ($ftphi - pow($de, 2) * $t10 + pow($de, 4) * $t11 - pow($de, 6) * $t12 + pow($de, 8) * $t13) * 180 / M_PI;
		return $result;
	}

	public static function getEpsg($lon, $lat = 0) {
		if ($lon < 13.5 || $lon >= 25.5) throw new Exception("Wrong longitude - {$lon}");
		return floor(($lon + 1.5) / 3) + 2171;
	}

	public static function LonLatToPUWGWGS84($lon, $lat) {
		$result = self::LonLatToPUWG(6378137, 1 / 298.257223563, $lon, $lat, 4);
		$result->epsg = self::getEpsg($lon);
		return $result;
	}

	public static function PUWGToLonLatWGS84($lon, $lat) {
		$result = self::PUWGToLonLat(6378137, 1 / 298.257223563, $lon, $lat, 4);
		$result->epsg = self::getEpsg($result->x);
		return $result;
	}

}