PHPで、ある緯度経度からある方角にある距離だけ移動したときの緯度経度を求める

# 東京から北西に10000km進むとどこに着くか？
`Vincenty法の逆解法`を使うと、二地点間の`方位角`を求めるのに使えますが、
`Vincenty法の順解法`をを使うことにより、ある地点からある方角にある距離移動したときの**緯度経度**を求められます。
PHPで組んでみます。

# Vincenty法の順解法

```php
/**
 * Vincenty法の順解法により、指定の地点に指定の方位角へ指定の測地線長を移動したときの緯度経度を求める。
 *
 * @param  float $lat1     緯度1 (十進法度単位)
 * @param  float $lng1     経度1 (十進法度単位)
 * @param  float $distance 加算する測地線長 (m)
 * @param  float $azimuth  加算する方位角 (ラジアン)
 * @return array           加算後の緯度経度 (十進法度単位) [0 => 緯度, 1 => 経度]
 */
function dest_vincenty($lat1, $lng1, $distance, $azimuth) {
	$earth_minor_radius = 6356752.0;
	$earth_major_radius = 6378137.0;
	$earth_flattening   = 0.0033528106811823;
	$acceptable_error   = 1E-12;

	$lat1        = deg2rad($lat1);
	$lng1        = deg2rad($lng1);
	$sin_azimuth = sin(-$azimuth);
	$cos_azimuth = cos(-$azimuth);

	// 測地緯度から化成緯度を求める
	$rlat1     = atan($earth_minor_radius / $earth_major_radius * tan($lat1));
	$sin_rlat1 = sin($rlat1);
	$cos_rlat1 = cos($rlat1);
	$tan_rlat1  = (1 - $earth_flattening) * tan($lat1);

	// 必要な値を事前に計算
	$sigma1     = atan2($tan_rlat1, $cos_azimuth);
	$sin_alpha  = $cos_rlat1 * $sin_azimuth;
	$cos2_alpha = 1 - $sin_alpha * $sin_alpha;
	$u2         = $cos2_alpha * (pow($earth_major_radius, 2) - pow($earth_minor_radius, 2)) / pow($earth_minor_radius, 2);
	$a          = 1 + $u2 / 16384 * (4096 + $u2 * (-768 + $u2 * (320 - 175 * $u2)));
	$b          = $u2 / 1024 * (256 + $u2 * (-128 + $u2 * (74 - 47 * $u2)));

	$sigma      = $distance / $earth_minor_radius / $a;
	$prev_sigma = $sigma + 1;
	$sin_sigma  = sin($sigma);
	$cos_sigma  = cos($sigma);
	$sigmam     = null;
	$cos_sigmam = null;
	// シグマ値を計算する。 (値が収束するまで繰り返す)
	for ($i = 0; $i<100 && abs($sigma - $prev_sigma) > $acceptable_error; $i++) {
		$sigmam     = 2 * $sigma1 + $sigma;
		$cos_sigmam = cos($sigmam);
		$tmp        = 1.0 / 6.0 * $b * $cos_sigmam * (-3 + 4 * pow($sin_sigma, 2)) * (-3 + 4 * pow($cos_sigmam, 2));
		$delta      = $b * $sin_sigma * ($cos_sigmam + 0.25 * $b * ($cos_sigma * (-1 + 2 * pow($cos_sigmam, 2)) - $tmp));
		$prev_sigma = $sigma;
		$sigma      = $distance / $earth_minor_radius / $a + $delta;
		$sin_sigma  = sin($sigma);
		$cos_sigma  = cos($sigma);
	}

	// 緯度を算出
	$tmp    = sqrt($sin_alpha * $sin_alpha + pow($sin_rlat1 * $sin_sigma - $cos_rlat1 * $cos_sigma * $cos_azimuth, 2));
	$lat2   = rad2deg(atan2($sin_rlat1 * $cos_sigma + $cos_rlat1 * $sin_sigma * $cos_azimuth, (1 - $earth_flattening) * $tmp));

	// 経度を算出
	$lambda = atan2($sin_sigma * $sin_azimuth, $cos_rlat1 * $cos_sigma - $sin_rlat1 * $sin_sigma * $cos_azimuth);
	$c      = $earth_flattening / 16.0 * $cos2_alpha * (4 + $earth_flattening * (4 - 3 * $cos2_alpha));
	$tmp    = ($sigma + $c * $sin_sigma * ($cos_sigmam + $c * $cos_sigma * (-1 + 2 * pow($cos_sigmam, 2))));
	$lng2   = rad2deg($lng1 - $lambda + (1 - $c) * $earth_flattening * $sin_alpha * $tmp);
	return [round($lat2 > 180 ? $lat2 - 360 : $lat2, 6), round($lng2 > 180 ? $lng2 - 360 : $lng2, 6)];
}
```

# 結果

```php
$tokyo = [35.681236, 139.767125];
$destination = dest_vincenty($tokyo[0], $tokyo[1], 10000000, -M_PI / 4);
echo "東京から北西に10000km移動したときの緯度経度: {$destination[0]}, {$destination[1]}\n";
```

```sh
東京から北西に10000km移動したときの緯度経度: 35.214532, 19.76941
```

結果は地中海のど真ん中でした。

※引数の方位角には`方位角のラジアン値`を指定します。
　これは、**真北**を`0`とし、**時計回りの方角**を`正の数`としたものです。`3.14...`が**真南**ということになります。
　度数からラジアン値に変換する場合は`deg2rad`関数を使ってください。