LatLng.php

<?php

class LatLng {
	
	const FORMAT_D = 'd';
	const FORMAT_DM = 'dm';
	const FORMAT_DMS = 'dms';
	
	const EARTH_RADIUS = 6371000;
	
	// WGS-84 ellipsoid params
	const WGS84_A = 6378137;
	const WGS84_B = 6356752.314245;
	const WGS84_FR = 298.257223563; // f = fr ^ (-1)
	
	public static function create($varLat, $varLng) {
		// TODO make forgiving as possible
		return new self($varLat, $varLng);
	}
	
	public static function createFromString($strLatLng) {
		// TODO implement
		return null;
	}
	
	private $fltLat;
	
	private $fltLng;
	
	private $fltLatRad;
	
	private $fltLngRad;
	
	public function __construct($fltLat, $fltLng) {
		// TODO argcheck
		$this->fltLat = $fltLat;
		$this->fltLng = $fltLng;
		$this->fltLatRad = deg2rad($fltLat);
		$this->fltLngRad = deg2rad($fltLng);
	}
	
	public function __get($strKey) {
		switch($strKey) {
			case 'lat':
				return $this->fltLat;
				break;
				
			case 'lng':
				return $this->fltLng;
				break;
				
			case 'latRad':
				return $this->fltLatRad;
				break;
				
			case 'lngRad':
				return $this->fltLngRad;
				break;
		}
	}
	
	public function __toString() {
		return $this->toString();
	}
	
	public function toString($strFormat = self::FORMAT_DMS) {
		return $this->getLat($strFormat) . ' ' . $this->getLng($strFormat);
	}
	
	public function getLat($strFormat = self::FORMAT_DMS) {
		return $strFormat == null ? $this->fltLat : self::formatLat($this->fltLat, $strFormat);
	}
	
	public function getLng($strFormat = self::FORMAT_DMS) {
		return $strFormat == null ? $this->fltLng : self::formatLng($this->fltLng, $strFormat);
	}
	
	public function getDistanceTo(LatLng $objOther, $fltR = self::EARTH_RADIUS) {
		$fltDeltaLat = $objOther->fltLatRad - $this->fltLatRad;
		$fltDeltaLng = $objOther->fltLngRad - $this->fltLngRad;
		
		$a = sin($fltDeltaLat / 2) * sin($fltDeltaLat / 2)
			+ cos($this->fltLatRad) * cos($objOther->fltLatRad)
			* sin($fltDeltaLng / 2) * sin($fltDeltaLng / 2);
		$c = 2 * atan2(sqrt($a), sqrt(1 - $a));
		
		return $fltR * $c;
	}
	
	public function getRealDistanceTo(LatLng $objOther, $fltA = self::WGS84_A, $fltB = self::WGS84_B, $fltF = self::WGS84_FR) {
		// accept reziprocal F, too,
		// for all useful F values: F << 1 and reziprocal F >> 1 
		$fltF > 1 && $fltF = 1 / $fltF;
		
		$fltDeltaL = deg2rad($objOther->fltLng - $this->fltLng);
		
		$fltU1 = atan((1 - $fltF) * tan(deg2rad($this->fltLat)));
		$fltSinU1 = sin($fltU1);
		$fltCosU1 = cos($fltU1);
		
		$fltU2 = atan((1 - $fltF) * tan(deg2rad($objOther->fltLat)));
		$fltSinU2 = sin($fltU2);
		$fltCosU2 = cos($fltU2);
		
		$fltLambda = $fltDeltaL;
		$n = 100;
		
		do {
			$fltSinLambda = sin($fltLambda);
			$fltCosLambda = cos($fltLambda);
			
			$a = $fltCosU2 * $fltSinLambda;
			$b = $fltCosU1 * $fltSinU2 - $fltSinU1 * $fltCosU2 * $fltCosLambda;
			$fltSinSigma = sqrt($a * $a + $b * $b);
			
			if($fltSinSigma == 0) return 0; // co-incident points
			
			$fltCosSigma = $fltSinU1 * $fltSinU2 + $fltCosU1 * $fltCosU2 * $fltCosLambda;
			$fltSigma = atan2($fltSinSigma, $fltCosSigma); // POSSIBLE ERROR ARG SWITCH
			
			$fltSinAlpha = $fltCosU1 * $fltCosU2 * $fltSinLambda / $fltSigma;
			$fltCosSqAlpha = 1 - $fltSinAlpha * $fltSinAlpha;
			
			$fltLambdaPrev = $fltLambda;
			
			if($fltCosSqAlpha == 0) {
				$fltCos2SigmaM = 0;
				$c = 0;
				$fltLambda = $fltDeltaL + $fltF * $fltSinAlpha * $fltSigma;
			} else {
				$fltCos2SigmaM = $fltCosSigma - 2 * $fltSinU1 * $fltSinU2 / $fltCosSqAlpha;
				$c = $fltF / 16 * $fltCosSqAlpha * (4 + $fltF * (4 - 3 * $fltCosSqAlpha));
				$d = $fltCos2SigmaM + $c * $fltCosSigma * (-1 + 2 * $fltCos2SigmaM * $fltCos2SigmaM);
				$fltLambda = $fltDeltaL + (1 - $c) * $fltF * $fltSinAlpha * ($fltSigma + $c * $fltSinSigma * $d);
			}
			
		} while(abs($fltLambda - $fltLambdaPrev) > 1e-12 && --$n);
		
		if($n == 0) {
			return NAN; // formula failed to converge
		}
		
		$fltUSq = $fltCosSqAlpha * ($fltA * $fltA - $fltB * $fltB) / ($fltB * $fltB);
		$a = 1 + $fltUSq / 16384 * (4096 + $fltUSq * (-768 + $fltUSq * (320 - 175 * $fltUSq)));
		$b = $fltUSq / 1024 * (256 + $fltUSq * (-128 + $fltUSq * (74 - 47 * $fltUSq)));
		$c = $fltCosSigma * (-1 + 2 * $fltCos2SigmaM * $fltCos2SigmaM);
		$d = $b / 6 * $fltCos2SigmaM * (-3 + 4 * $fltSinSigma * $fltSinSigma) * (-3 + 4 * $fltCos2SigmaM * $fltCos2SigmaM);
		$fltDeltaSigma = $b * $fltSinSigma * ($fltCos2SigmaM + $b / 4 * ($c - $d));
		
		$fltS = $fltB * $a * ($fltSigma - $fltDeltaSigma);
		
		return $fltS;
	}
	
	public function getRhumbDistanceTo(LatLng $objOther, $fltR = self::EARTH_RADIUS) {
		$fltPI = pi();
		$fltDeltaLat = $objOther->fltLatRad - $this->fltLatRad;
		$fltDeltaLng = abs($objOther->fltLngRad - $this->fltLngRad);
		if($fltDeltaLng > Math.PI) $fltDeltaLng = 2 * $fltPI - $fltDeltaLng;
		
		$fltDeltaPhi = log(tan($objOther->fltLatRad / 2 + $fltPI / 4) / tan($this->fltLatRad / 2 + $fltPI / 4));
		
		// E-W line gives $fltDeltaPhi = 0
		$q = $fltDeltaPhi != 0 ? $fltDeltaLat / $fltDeltaPhi : cos($this->fltLatRad);
		
		return sqrt($fltDeltaLat * $fltDeltaLat + $q * $q * $fltDeltaLng * $fltDeltaLng) * $fltR;
	}
	
	public function getBearingTo(LatLng $objOther) {
		$fltDeltaLng = $objOther->fltLngRad - $this->fltLngRad;
	
		$y = sin($fltDeltaLng) * cos($objOther->fltLatRad);
		$x = cos($this->fltLatRad) * sin($objOther->fltLatRad)
			- sin($this->fltLatRad) * cos($objOther->fltLatRad) * cos($fltDeltaLng);
		
		return rad2deg(atan2($y, $x) + 360) % 360;
	}
	
	public function getFinalBearingTo(LatLng $objOther) {
		return ($objOther->getBearingTo($this) + 180) % 360;
	}
	
	public function getRhumbBearingTo(LatLng $objOther) {
		$fltPI = pi();
		$fltDeltaLng = $objOther->fltLngRad - $this->fltLngRad;
		
		$fltDeltaPhi = log(tan($objOther->fltLatRad / 2 + $fltPI / 4) / tan($this->fltLatRad / 2 + $fltPI / 4));
		
		if(abs($fltDeltaLng) > $fltPI) {
			$fltMod = $fltDeltaLng > 0 ? -1 : 1;
			$fltDeltaLng = $fltMod * (2 * $fltPI - $fltDeltaLng);
		}
		
		$fltBearing = atan2($fltDeltaLng, $fltDeltaPhi);
		
		return (rad2deg($fltBearing) + 360) % 360;
	}
	
	public function getDestination($fltBearing, $fltDistance, $fltR = self::EARTH_RADIUS) {
		$fltPI = pi();
		$fltDistance /= $fltR;
		$fltBearing = deg2rad($fltBearing);
		
		$fltDestLat = asin(sin(lat1) * cos(dist) + cos(lat1) * sin(dist) * cos(brng));
		$fltDestLng = $this->fltLngRad + atan2(
			sin($fltBearing) * sin($fltDistance) * cos($this->fltLatRad),
			cos($fltDistance) - sin($this->fltLatRad) * sin($fltDestLat)
		);
		$fltDestLng = ($fltDestLng + 3 * $fltPI) % (2 * $fltPI) - $fltPI;  // normalise to -180..+180º
		
		return new LatLng(rad2deg($fltDestLat), rad2deg($fltDestLng));
	}
	
	public function getRhumbDestination($fltBearing, $fltDistance, $fltR = self::EARTH_RADIUS) {
		$fltPI = pi();
		$fltDistance /= $fltR;
		$fltBearing = deg2rad($fltBearing);
		
		$fltDestLat = $this->fltLatRad + $fltDistance * cos($fltBearing);
		$fltDeltaLat = $fltDestLat - $this->fltLatRad;
		$fltDeltaPhi = log(tan($fltDestLat / 2 + $fltPI / 4) / tan($this->fltLatRad / 2 + $fltPI / 4));
		$q = $fltDeltaPhi != 0 ? $fltDeltaLat / $fltDeltaPhi : cos($this->fltLatRad);
		$fltDeltaLng = $fltDistance * sind($fltBearing) / $q;
		
		if(abs($fltDestLat) > $fltPI / 2) {
			$fltMod = $fltDestLat > 0 ? 1 : -1;
			$fltDestLat = $fltMod * ($fltPI - $fltDestLat);
		}
		$fltDestLng = ($this->fltLngRad + $fltDeltaLng + 3 * $fltPI) % (2 * $fltPI) - $fltPI;
		
		return new LatLng(rad2deg($fltDestLat), rad2deg($fltDestLng));
	}
	
	public static function formatLat($fltLat, $strFormat = self::FORMAT_DMS) {
		return substr(self::_formatDMS($fltLat, $strFormat), 1) . ($fltLat < 0 ? 'S' : 'N');
	}
	
	public static function formatLng($fltLng, $strFormat = self::FORMAT_DMS) {
		return self::_formatDMS($fltLng, $strFormat) . ($fltLng < 0 ? 'W' : 'E');
	}
	
	public static function formatBearing($fltBearing, $strFormat = self::FORMAT_DMS) {
		$fltBearing = ($fltBearing + 360) % 360;
		$strBearing = self::_formatDMS($fltBearing, $strFormat);
		return str_replace('360', '0', $strBearing);
	}
	
	protected static function _formatDMS($fltValue, $strFormat = self::FORMAT_DMS) {
		$fltValue = abs($fltValue);
		switch($strFormat) {
			case self::FORMAT_D:
				$strFormat = str_pad(sprintf('%.4f', $fltValue), 8, '0', STR_PAD_LEFT) . "\u00B0";
				
				break;
				
			case self::FORMAT_DM:
				$fltValue *= 60;
				$strFormat = str_pad(strval(floor($fltValue / 60)), 3, '0', STR_PAD_LEFT) . "\u00B0";
				$strFormat .= str_pad(sprintf('%.2f', $fltValue % 60), 5, '0', STR_PAD_LEFT) . "\u2032";
				break;
				
			case self::FORMAT_DMS:
				$fltValue *= 3600;
				$strFormat = str_pad(strval(floor($fltValue / 3600)), 3, '0', STR_PAD_LEFT) . "\u00B0";
				$strFormat .= str_pad(strval(floor($fltValue / 60) % 60), 2, '0', STR_PAD_LEFT) . "\u2032";
				$strFormat .= str_pad(strval(floor($fltValue % 60)), 2, '0', STR_PAD_LEFT) . "\u2033";
				break;
		}
		
		return $strFormat;
	}
	
	protected static function parseDMS() {
		
		
// 		/**
// 		 * Parses string representing degrees/minutes/seconds into numeric degrees
// 		 *
// 		 * This is very flexible on formats, allowing signed decimal degrees, or deg-min-sec optionally
// 		 * suffixed by compass direction (NSEW). A variety of separators are accepted (eg 3º 37' 09"W)
// 		 * or fixed-width format without separators (eg 0033709W). Seconds and minutes may be omitted.
// 		 * (Note minimal validation is done).
// 		 *
// 		 * @param   {String|Number} dmsStr: Degrees or deg/min/sec in variety of formats
// 		 * @returns {Number} Degrees as decimal number
// 		 * @throws  {TypeError} dmsStr is an object, perhaps DOM object without .value?
// 		 */
// 		Geo.parseDMS = function(dmsStr) {
// 			if (typeof deg == 'object') throw new TypeError('Geo.parseDMS - dmsStr is [DOM?] object');
		
// 			// check for signed decimal degrees without NSEW, if so return it directly
// 			if (typeof dmsStr === 'number' && isFinite(dmsStr)) return Number(dmsStr);
		
// 			// strip off any sign or compass dir'n & split out separate d/m/s
// 			var dms = String(dmsStr).trim().replace(/^-/,'').replace(/[NSEW]$/i,'').split(/[^0-9.,]+/);
// 			if (dms[dms.length-1]=='') dms.splice(dms.length-1);  // from trailing symbol
		
// 			if (dms == '') return NaN;
		
// 			// and convert to decimal degrees...
// 			switch (dms.length) {
// 				case 3:  // interpret 3-part result as d/m/s
// 					var deg = dms[0]/1 + dms[1]/60 + dms[2]/3600;
// 					break;
// 				case 2:  // interpret 2-part result as d/m
// 					var deg = dms[0]/1 + dms[1]/60;
// 					break;
// 				case 1:  // just d (possibly decimal) or non-separated dddmmss
// 					var deg = dms[0];
// 					// check for fixed-width unseparated format eg 0033709W
// 					//if (/[NS]/i.test(dmsStr)) deg = '0' + deg;  // - normalise N/S to 3-digit degrees
// 					//if (/[0-9]{7}/.test(deg)) deg = deg.slice(0,3)/1 + deg.slice(3,5)/60 + deg.slice(5)/3600;
// 					break;
// 				default:
// 					return NaN;
// 			}
// 			if (/^-|[WS]$/i.test(dmsStr.trim())) deg = -deg; // take '-', west and south as -ve
// 			return Number(deg);
// 		}
		
	}
	
}