Little helper: Geodistance

Today I present you a helper class int PHP offering methods for the most famous variants in calculating geographical distances by latitude and longitude (Circle, Cosines, Haversine and Vincenty included) as well as cartesian coordinates (a great alternative for use in SQL-queries).

During development of some PHP and JS applications using geolocations and distances I finally found the time to refactor the knowledge I gained into a tiny PHP-helper-class offering fast (but not that accurate) to slow (but accurate to the max) static methods for your daily geodistance calculations.

<?php
	class Geographical {
		const METHOD_CIRCLE        = 1;
		const METHOD_COSINES       = 2;
		const METHOD_HAVERSINE     = 3;
		const METHOD_VINCENTY      = 4;
		const RADIUS_METER         = 6371000.785;
		const SEMIAX_MAJOR         = 6378137;
		const SEMIAX_MINOR         = 6356752.3141;

			public static function getCartesian($longitude, $latitude) {
			$lambda    = $longitude * pi() / 180;
			$phi       = $latitude * pi() / 180;
			$return    = new stdClass();
			$return->x = self::RADIUS_METER * cos($phi) * cos($lambda);
			$return->y = self::RADIUS_METER * cos($phi) * sin($lambda);
			$return->z = self::RADIUS_METER * sin($phi);

				unset($lambda);
				unset($phi);

				return $return;
			}

			public static function getDistanceByCartesian($x1, $y1, $z1, $x2, $y2, $z2) {
				return 2 * self::RADIUS_METER *
					asin(
						sqrt(
							pow($x1 - $x2, 2)
							+ pow($y1 - $y2, 2)
							+ pow($z1 - $z2, 2)
						) / (2 * self::RADIUS_METER)
					);
			}

			public static function getDistanceByLatLng($lat1, $lng1, $lat2, $lng2, $method = self::METHOD_COSINES) {
				$return = false;
					switch($method) {
						case self::METHOD_CIRCLE:
							$factor  = deg2rad(1);
							$lat1   *= $factor;
							$lng1   *= $factor;
							$lat2   *= $factor;
							$lng2   *= $factor;

							$return = rad2deg(acos(sin($lat1) * sin($lat2) + cos($lat1) * cos($lat2))) * 60 * 1852;

							unset($factor);
							unset($lat1);
							unset($lng1);
							unset($lat2);
							unset($lng2);

							break;
						case self::METHOD_COSINES:
							$factor  = deg2rad(1);
							$lat1   *= $factor;
							$lng1   *= $factor;
							$lat2   *= $factor;
							$lng2   *= $factor;

							$return = self::RADIUS_METER * acos(sin($lat1) * sin($lat2) + cos($lat1) * cos($lat2) * cos($lng1 - $lng2));

							unset($factor);
							unset($lat1);
							unset($lng1);
							unset($lat2);
							unset($lng2);

							break;
						case self::METHOD_HAVERSINE:
							$factor  = deg2rad(1);
							$lat1   *= $factor;
							$lng1   *= $factor;
							$lat2   *= $factor;
							$lng2   *= $factor;
							$dlat    = $lat2 - $lat1;
							$dlng    = $lng2 - $lng1;

							$a = sin($dlat / 2) * sin($dlat / 2) + cos($lat1) * cos($lat2) * sin($dlng / 2) * sin($dlng / 2);

							$return = self::RADIUS_METER * (2 * atan2(sqrt($a), sqrt(1 - $a)));

							unset($factor);
							unset($lat1);
							unset($lng1);
							unset($lat2);
							unset($lng2);
							unset($dlat);
							unset($dlng);
							unset($a);

							break;
						case self::METHOD_VINCENTY:
							$factor  = deg2rad(1);
							$lat1   *= $factor;
							$lng1   *= $factor;
							$lat2   *= $factor;
							$lng2   *= $factor;
							$f       = (self::SEMIAX_MAJOR - self::SEMIAX_MINOR) / self::SEMIAX_MAJOR;
							$L       = $lng2 - $lng1;
							$U1      = atan((1 - $f) * tan($lat1));
							$U2      = atan((1 - $f) * tan($lat2));
							$sinU1   = sin($U1);
							$sinU2   = sin($U2);
							$cosU1   = cos($U1);
							$cosU2   = cos($U2);
							$lambda  = $L;
							$lambdaP = 2 * pi();
							$i       = 20;

							while(abs($lambda - $lambdaP) > 1e-12 and --$i > 0) {
								$sinLambda = sin($lambda);
								$cosLambda = cos($lambda);
								$sinSigma  = sqrt(($cosU2 * $sinLambda) * ($cosU2 * $sinLambda) + ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda) * ($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda));

								if($sinSigma == 0) {
									return 0;
								}

								$cosSigma   = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosLambda;
								$sigma      = atan2($sinSigma, $cosSigma);
								$sinAlpha   = $cosU1 * $cosU2 * $sinLambda / $sinSigma;
								$cosSqAlpha = 1 - $sinAlpha * $sinAlpha;
								$cos2SigmaM = $cosSigma - 2 * $sinU1 * $sinU2 / $cosSqAlpha;

								if(is_nan($cos2SigmaM)) {
									$cos2SigmaM = 0;
								}

								$c       = $f / 16 * $cosSqAlpha * (4 + $f * (4 - 3 * $cosSqAlpha));
								$lambdaP = $lambda;
								$lambda  = $L + (1 - $c) * $f * $sinAlpha * ($sigma + $c * $sinSigma * ($cos2SigmaM + $c * $cosSigma * (-1 + 2 * $cos2SigmaM * $cos2SigmaM)));
							}

							if($i == 0) {
								return false;
							}

							$uSq        = $cosSqAlpha * (self::SEMIAX_MAJOR * self::SEMIAX_MAJOR - self::SEMIAX_MINOR * self::SEMIAX_MINOR) / (self::SEMIAX_MINOR * self::SEMIAX_MINOR);
							$A          = 1 + $uSq / 16384 * (4096 + $uSq * (-768 + $uSq * (320 - 175 * $uSq)));
							$B          = $uSq / 1024 * (256 + $uSq * (-128 + $uSq * (74 - 47 * $uSq)));
							$deltaSigma = $B * $sinSigma * ($cos2SigmaM + $B / 4 * ($cosSigma * (-1 + 2 * $cos2SigmaM * $cos2SigmaM) - $B / 6 * $cos2SigmaM * (-3 + 4 * $sinSigma * $sinSigma) * (-3 + 4 * $cos2SigmaM * $cos2SigmaM)));

							$return = self::SEMIAX_MINOR * $A * ($sigma - $deltaSigma);

							break;
					}

					return $return;
				}
			}
			?>

I took the code from various sources and refactored it to my needs. For all the calculations taking latitude and longitude directly being quite complex I also added a more simple method using cartesian coordinates and also integrated a method to convert latitude and longitude to cartesian coordinates. The Vincenty-method is by far the most accurate but takes quite some time. I prefer to use the Haversine-formula for it’s good compromise between accuracy and speed.