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

東京から北西に10000km進むとどこに着くか？

Vincenty法の逆解法を使うと、二地点間の方位角を求めるのに使えますが、
Vincenty法の順解法をを使うことにより、ある地点からある方角にある距離移動したときの緯度経度を求められます。
PHPで組んでみます。

Vincenty法の順解法

/**
 * 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)];
}

結果

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

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

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

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