<?php

// By Douglas Myers-Turnbull
/*
 * This file is part of VincentyCalculator. VincentyCalculator is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version. VincentyCalculator is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See
 * the GNU General Public License for more details. You should have received a copy of the GNU General Public License
 * along with VincentyCalculator. If not, see <http://www.gnu.org/licenses/>.
 */

// Set the internal encoding
mb_language('uni');
mb_internal_encoding('UTF-8');

// The values for ELLIPSOID_A, ELLIPSOID_B, and ELLIPSOID_F are given by WGS84
// The code should also fit any other ellipsoidal model

// Does not affect the calculation
// Only used for display
define("LENGTH_UNIT", "meters");

// A is the ellipsoid’s major semiaxis
define("ELLIPSOID_A", 6378137);

// B is the ellipsoid’s minor semiaxis
define("ELLIPSOID_B", 6356752.314245);

// F is flattening, given by (a − b) / a
// It is here to avoid computation time
define("ELLIPSOID_F", 1/298.257223563);

// The maximum number of iterations independent of Δλ
define("MAX_ITERATIONS", 20);

// The value considered to be negligible change (Δλ)
// If Δλ < this value, the iterations cease
define("NEGLIGIBLE_CHANGE", .00000000001);

findDistance();

/**
 * Finds the distance between two points and sends a location header to display the results.
 */
function findDistance() {

	// Get each parameter in a PHP‐recognized format
	$lon1 = formatParam(replaceChars(trim(getGet('lon1'))), 'lon', $lon1err);
	$lat1 = formatParam(replaceChars(trim(getGet('lat1'))), 'lat', $lat1err);
	$lon2 = formatParam(replaceChars(trim(getGet('lon2'))), 'lon', $lon2err);
	$lat2 = formatParam(replaceChars(trim(getGet('lat2'))), 'lat', $lat2err);

	// formatParam returns false on error
	if (!$lon1 || !$lat1 || !$lon2 || !$lat2) {
		$error_str = 'bad&' . ($lon1err ? 'lon1err&' : '') . ($lon2err ? 'lon2err&' : '') . ($lat1err ? 'lat1err&' : '') . ($lat2err ? 'lat2err&' : '');
		header("Location: index.html?" . $error_str . "lon1=" . $_GET['lon1'] . "&lat1=" . $_GET['lat1'] . "&lon2=" . $_GET['lon2'] . "&lat2=" . $_GET['lat2']);
		die();
	}

	// Find the distance and format it
	$distance = calculateDistance($lon1, $lat1, $lon2, $lat2);
	$result = formatDistance($distance);

	header('Location: index.html?lon1=' . $_GET['lon1'] . '&lat1=' . $_GET['lat1'] . '&lon2=' . $_GET['lon2'] . '&lat2=' . $_GET['lat2'] . '&result=' . $result . '&units=' . LENGTH_UNIT);
	//header("Location: index.html?lon1=$lon1&lat1=$lat1&lon2=$lon2&lat2=$lat2&result=$result&units=" . LENGTH_UNIT);
	die();

}

/**
 * Replaces certain characters with characters with equivalent meanings.
 * @param $param The string in which to search
 * @return The new string with replaced characters
 */
function replaceChars($param) {

	/*
	Equivalent characters:
		′ (U+2032) ⟷ ' (U+0027)
		″ (U+2033) ⟷ " (U+0022)
		− (U+2212) ⟷ - (U+002D)
	*/

	$equiv = array('′' => '\'', '″' => '"', '−' => '-');

	$result = $param;

	foreach ($equiv as $index => $value) {
		$result = str_replace($index, $value, $result);
	}

	return $result;

}

/**
 * Formats a parameter for use in Vincenty’s formula. Takes a direction (N, E, S, W) into account and combines degrees, minutes, and seconds.
 * @param $param The parameter to format
 * @param $type The type of the parameter; can be either 'lon' for a longitude or 'lat' for a latitude
 * @param $error Passed by reference; a comma-separated list of problems
 * @return The formatted parameter in a PHP‐recognized numeric format or false on error
 */
function formatParam($param, $type, &$error) {

/*
Correct:
	34.5S       ⟶ −34.5
Both negative sign and direction:
	−34.5S      ⟶ −34.5, warning
	−34.5N      ⟶   34.5, warning
Invalid direction (assume a latitude):
	34.5E       ⟶   error
	34.5W       ⟶ error
Correct:
	34 30       ⟶  34.5
	34 30 30    ⟶  34.55
	34 30′30″    ⟶  34.55
	34 30′ 30″   ⟶  34.55
Negative minutes or seconds:
	34 −30      ⟶  33.5, warning
	34 −30 −15  ⟶ 33.575, warning
Misuse of prime marks:
	34 30′ 15′    ⟶  34.525, warning
	34 30″ 15″   ⟶  34.525, warning
	34 30″ 15′   ⟶  34.525, warning
*/

	// If something is wrong, false must always be returned
	if ($param === false) {
		return false;
	}

	$param = convertDirection($param, $type, $error);
	$param = toDecimalFormat($param, $error);

	return $param;

}

/**
 * Accounts for a direction specified by N, E, S, and W.
 * @param $value The value to format
 * @param $type The type of the parameter; can be either 'lon' for a longitude or 'lat' for a latitude
 * @param $error Passed by reference; a comma-separated list of problems
 * @return The resulting value
 */
function convertDirection($value, $type, &$error) {

	$dir = strtolower($value[strlen($value) - 1]);

	// Only convert for values with a direction
	if (!is_numeric($dir) && $dir != "'" && $dir != '"') {

		// Check longitudes
		if ($type == 'lon') {
			if ($dir == 'w') {
				return -substr($value, 0, strlen($value) - 1);
			} else if ($dir == 'e') {
				return substr($value, 0, strlen($value) - 1);
			} else {
				$error .= "dir$dir,";
				return false;
			}
		}

		// Check latitudes
		if ($type == 'lat') {
			if ($dir == 's') {
				return -substr($value, 0, strlen($value) - 1);
			} else if ($dir == 'n') {
				return substr($value, 0, strlen($value) - 1);
			} else {
				$error .= "dir$dir,";
				return false;
			}
		}

	}

	return $value;

}

/**
 * Convert the format [degrees minutes seconds] to [degrees].
 * @param $value The value to format
 * @param $error Passed by reference; a comma-separated list of problems
 * @return The converted value
 */
function toDecimalFormat($value, &$error) {

	$result = 0;

	// Single quote must be before double quote
	if (stripos($value, '\'') !== false && stripos($value, '"') !== false && stripos($value, '\'') > stripos($value, '"'))
	{
		$err .= 'qte,';
		return false;
	}

	// Because dd dd'dd" is accepted, convert quote and double quote to a space
	$piece = explode(' ', str_replace("\\'", ' ', str_replace("\\\"", ' ', $value)));

	if (count($piece) > 0) {

		// The number of components (degrees, minutes, seconds) included into $result
		$a = 0;

		foreach ($piece as $i => $v) {

			// $v can be blank with formats such as dd dd' dd"
			if ($v != '') {

				if (is_numeric($v)) {
					// Add the component to $result
					$result += $v / pow(60, $a);
				} else {
					$error .= 'chr,';
					return false;
				}

				$a++;

			}

		}

	} else {

		// If minutes and seconds were not specified, return $value as‐is
		$result = $value;

	}

	return $result;

}

/**
 * @return The formatted distance, in the form: nn,nnn.ddd ddd d units.
 * @param The distance as a number
 */
function formatDistance($distance) {

	// Format the number
	$index = (strlen($distance) - strpos($distance, '.'));
	$result = number_format($distance, $index);

	// Split the decimal portion by thousands with spaces
	// number_format will NOT do this :(
	$index = strpos($result, '.');
	$buffer = substr($result, 0, $index);
	for ($i = 0; $i < strlen($distance); $i++) {
		$buffer .= $result[$i + $index];
		if ($i % 3 == 0 && $i != 0) {
			$buffer .= ' ';
		}
	}

	return trim(trim($buffer), '.');

}

/**
 * @return The value of a GET parameter specified by $name or false on error
 * @param $name The name of the GET parameter
 */
function getGet($name) {
	if (isset($_GET[$name]) && $_GET[$name] != '') {
		return $_GET[$name];
	}
	return false;
}

/**
 * Finds the distance between two points using Vincenty’s formula.
 * @param $lon1 The first point’s longitude, in radians
 * @param $lat1 The first point’s latitude, in radians
 * @param $lon2 The second point’s longitude, in radians
 * @param $lat2 The first point’s latitude, in radians
 * @return The geodesic’s length or false on failure
 */
function calculateDistance($lon1, $lat1, $lon2, $lat2) {

	// L is the difference in latitude (positive east)
	// The longitudes are used nowhere else
	$L = toRadians($lon2 - $lon1);

	// U is reduced latitude, (1−f)⋅tan(φ), where φ₁ ≍ $lat₁ and φ₂ ≍ $lat₂
	// Here, phi is $lat1 and $lat2
	$U1 = atan((1 - ELLIPSOID_F) * tan(toRadians($lat1)));
	$U2 = atan((1 - ELLIPSOID_F) * tan(toRadians($lat2)));

	// Compute the sine and cosine of U1 and U2 beforehand because they are used often
	$cosU1 = cos($U1);
	$cosU2 = cos($U2);
	$sinU1 = sin($U1);
	$sinU2 = sin($U2);

	// Described in §4, step 13
	// Lambda (λ) is the difference in longitude on an auxillary sphere
	// Set λ to L for a first approximation
	// Lambda will change—and become a closer approximation—each iteration
	$lambda = $L;

	// Iterate until either:
	// A) Δλ is negligible (here, Δλ < NEGLIGIBLE_CHANGE), which indicates that minimal error exists
	// B) $iter > MAX_ITERATIONS:
	// The needed number of iterations before Δλ < NEGLIGIBLE_CHANGE is loosely
	// proportional to how close to antipodal the two points are
	$iter = 0;
	do {

		// Calculate the sine and cosine of λ beforehand because they are used often
		// Must calculate them each iteration because λ changes
		$sinLambda = sin($lambda);
		$cosLambda = cos($lambda);

		// Described in §4, steps 14, 15, and 16, respectively
		// Sigma (σ) is angular distance on the axuillary sphere
		// Vincenty’s paper uses sin²σ instead of sinσ
		$sinSigma = sqrt(sqr($cosU2 * $sinLambda) + sqr($cosU1 * $sinU2 - $sinU1 * $cosU2 * $cosLambda));

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

		$cosSigma = $sinU1 * $sinU2 + $cosU1 * $cosU2 * $cosLambda;
		$tanSigma = $sinSigma / $cosSigma;

		// Described in §4, step 17
		// Alpha (α) is the angular distance between the two points on the auxillary sphere
		$sinAlpha = $cosU1 * $cosU2 * $sinLambda / $sinSigma;

		// Needed in step 18
		$cosSqrAlpha = 1 - sqr($sinAlpha);

		// Described in §4, step 18
		// σm is the angular distance on the auxillary sphere from the equator to the line’s midpoint
		// On the equator, cos²α = 0: in such case, cos2σm should be zero (mentioned in §6)
		$cos2SigmaM = 0;
		if ($cosSqrAlpha != 0) {
			$cos2SigmaM = $cosSigma - 2 * $sinU1 * $sinU2 / $cosSqrAlpha;
		}

		// Described in §3, step 10
		$C = (ELLIPSOID_F / 16) * $cosSqrAlpha * (4 + ELLIPSOID_F * (4 - 3 * $cosSqrAlpha));

		// Remember λ to find Δλ, which will be |λ − λ₀|
		$lambdaLast = $lambda;

		// This is needed in step 11
		$sigma = atan2($sinSigma, $cosSigma);

		// Described in §3, step 11
		$lambda = $L + (1 - $C) * ELLIPSOID_F * $sinAlpha *
		($sigma + $C * $sinSigma * ($cos2SigmaM + $C * $cosSigma * (-1 + 2 * sqr($cos2SigmaM))));
		$iter++;

		// Let the calculation fail at MAX_ITERATIONS iterations
		// Will occur when |λ| > π
		if ($iter > MAX_ITERATIONS) {
			return false;
		}

		$iter++;

	} while (abs($lambda - $lambdaLast) > NEGLIGIBLE_CHANGE);

	// Needed in §3, step 3
	// As explained in the definitions, u² = cos²α(a² − b²) / b²
	$uSqr = $cosSqrAlpha * (sqr(ELLIPSOID_A) - sqr(ELLIPSOID_B)) / sqr(ELLIPSOID_B);

	// Described in §3, step 3
	$A = 1 + $uSqr / 16384 * (4096 + $uSqr * (-768 + $uSqr * (320 - 175 * $uSqr)));

	// Described in §3, step 4
	$B = $uSqr / 1024 * (256 + $uSqr * (-128 + $uSqr * (74 - 47 * $uSqr)));

	// Described in §3, step 6
	$deltaSigma = $B * $sinSigma * ($cos2SigmaM + ($B / 4) * ($cosSigma * (-1 + 2 * sqr($cos2SigmaM)) -
	$B / 6 * $cos2SigmaM * (-3 + 4 * sqr($sinSigma)) * (-3 + 4 * sqr($cos2SigmaM))));

	// Described in §4, step 19
	// Ess (s) is the length of the geodesic ☺
	$ess = ELLIPSOID_B * $A * ($sigma - $deltaSigma);

	// Steps 20 and 21 in §4 are not needed
	return $ess;

}

/**
 * Converts degrees to radians.
 * @param $degrees The angle in degrees
 * @return The angle in radians
 */
function toRadians($degrees) {
	return $degrees * M_PI / 180;
}

/**
 * Squares a number.
 * @param $number The number to square
 * @return The squared number
 */
function sqr($number) {
	return $number * $number;
}

?>