<?php
// datum converter
// Original: Nowral 01/3/9
//---------------------
//public
//Settings
$pi	= 4 * atan2(1,1);	// pi
$rd	= $pi / 180;		// radian

//---------------------
//private
function geoconv($b, $l, $h, $dat, $dat_){
	//datum
	$a_WGS84 = 6378137.0; //WGS-84 Equatorial Radius (a)
	$f_WGS84 = 1 / 298.257223563; //WGS-84 Flattening (f)
	
	$lines = load_datumlist();
	
	if(!$lines==false){
		foreach ($lines as $r => $line) {
			list($dat1, $ellp, $da, $df, $dx0, $dy0, $dz0, $datname) = split(',', $line);
			$a{$datname} = $a_WGS84 - $da;
			$f = $f_WGS84 - $df/10000;
			$e2{$datname} = 2*$f - $f*$f; # First Eccentricity
			$dx{$datname} = $dx0;
			$dy{$datname} = $dy0;
			$dz{$datname} = $dz0;
		}
	}
	/*$file_name = 'datum.csv';
	if(is_file($file_name)){
		$text = file_get_contents($file_name);
		$lines = explode("\n", $text);
		foreach ($lines as $r => $line) {
			list($dat1, $ellp, $da, $df, $dx0, $dy0, $dz0, $datname) = split(',', $line);
			$a{$datname} = $a_WGS84 - $da;
			$f = $f_WGS84 - $df/10000;
			$e2{$datname} = 2*$f - $f*$f; # First Eccentricity
			$dx{$datname} = $dx0;
			$dy{$datname} = $dy0;
			$dz{$datname} = $dz0;
		}
	}*/
	// target
	if(isset($a{$dat_})&&isset($a{$dat})){
		$a_ = $a{$dat_};
		$e2_ = $e2{$dat_};
		$dx_ = $dx{$dat_};
		$dy_ = $dy{$dat_};
		$dz_ = $dz{$dat_};
		
		$dx{$dat} = $dx{$dat} - $dx_;
		$dy{$dat} = $dy{$dat} - $dy_;
		$dz{$dat} = $dz{$dat} - $dz_;
	}else{
		echo 'Error: Illegal parameter.';
		exit();
	}
	//convert
	list($x, $y, $z) = llh2xyz($b, $l, $h, $a{$dat}, $e2{$dat});
	return xyz2llh($x+$dx{$dat}, $y+$dy{$dat}, $z+$dz{$dat}, $a_, $e2_);
	#list($b, $l, $h) = xyz2llh($x+$dx{$dat}, $y+$dy{$dat}, $z+$dz{$dat}, $a_, $e2_);
	#return array($b, $l, $h);
}
//---------------------
function load_datumlist(){
	$file_name = 'datum.csv';
	if(is_file($file_name)){
		$text = file_get_contents($file_name);
		$lines = explode("\n", $text);
		return $lines;
	}else{
		return false;
	}
}
//---------------------
function get_datum_name(){
	$lines = load_datumlist();
	if(!$lines==false){
		foreach ($lines as $r => $line) {
			list($dat1, $ellp, $da, $df, $dx0, $dy0, $dz0, $datname) = split(',', $line);
			$datum_name{$datname} = $datname;
		}
		return $datum_name;
	}
}
//---------------------
function llh2xyz($b, $l, $h, $a, $e2) { # 楕円体座標 -> 直交座標
	global $rd;
	
	$b *= $rd;
	$l *= $rd;
	$sb = sin($b);
	$cb = cos($b);
	$rn = $a / sqrt(1-$e2*$sb*$sb);
	
	$x = ($rn+$h) * $cb * cos($l);
	$y = ($rn+$h) * $cb * sin($l);
	$z = ($rn*(1-$e2)+$h) * $sb;
	
	return array($x, $y, $z);
}
//---------------------
function xyz2llh($x, $y, $z, $a, $e2) { # 直交座標 -> 楕円体座標
	global $rd;
	$bda = sqrt(1-$e2); # b/a
	
	$p = sqrt($x*$x+$y*$y);
	$t = atan2($z, $p*$bda);
	$st = sin($t);
	$ct = cos($t);
	$b = atan2($z+$e2*$a/$bda*$st*$st*$st, $p-$e2*$a*$ct*$ct*$ct);
	$l = atan2($y, $x);
	
	$sb = sin($b);
	$rn = $a / sqrt(1-$e2*$sb*$sb);
	$h = $p/cos($b) - $rn;
	
	return array($b/$rd, $l/$rd, $h);
}
//---------------------
function deg2gdms($d) { # GPSy format
	$d = ($d>180)? $d-360: $d; # Latitude, Longitude
	if($d<0) {
		$sgn = "-";
		$d *= -1;
	}else{
		$sgn = "";
	}
	$sf = round($d*36000);
	$s = round($sf/10) % 60;
	$m = round($sf/600) % 60;
	$d = round($sf/36000);
	$sf %= 10;
	return sprintf("%s%d.%02d'%02d.%d\"", $sgn, $d, $m, $s, $sf);
}
//---------------------
function deg2dms($d) {
	$sf = round($d*360000);
	$s = $sf / 100 % 60;
	$m = $sf / 6000 % 60;
	$d = floor($sf/360000);
	$sf %= 100;
	return sprintf("%d.%02d'%02d.%02d\"", $d, $m, $s, $sf);
}
//---------------------
/*$debug=0;
if($debug){
	// ex.latlng of Tokyo : Lat 35°20′39.984328", Lon 138°35′08.086122", Height 697.681000 [m]
	#$b = 36.84006462037767;			// Latitude(decimal)
	#$l = 140.0042724609375;			// Longitude(decimal)
	#$b = 35 + 20/60 + 39.984328/3600;	// Latitude
	#$l = 138 + 35/60 +	8.086122/3600;	// Longitude
	$h = 697.681000;					// Height [m]
	$h = (isset($_GET['h']))? $_GET['h']:$h;
	$latlng = "-36.84006462037767,-140.0042724609375";
	list($b,$l) = (isset($_GET['latlng']))? explode(",", htmlspecialchars($_GET['latlng'])):explode(",", $latlng);

	//datum
	$dat  = (isset($_GET['datum'])) ? $_GET['datum']: "Tokyo";	// from default:Tokyo
	$dat_ = (isset($_GET['datum_']))? $_GET['datum_']:"WGS 84";	// to	 default:WGS 84

	echo "<p>Input \nLat： ".htmlspecialchars(deg2dms($b))."\nLong： ".htmlspecialchars(deg2dms($l))."\nHeight： ".htmlspecialchars($h)."</p>\n\n";
	list($b, $l, $h) = geoconv($b, $l, $h, $dat, $dat_);
	echo "<p>Output \nLat： ".htmlspecialchars(deg2gdms($b))."\nLong： ".htmlspecialchars(deg2gdms($l))."\nHeight： ".htmlspecialchars($h)."</p>\n\n";
}*/