「地球上の２点間の最短距離」を求めるPHPコード：【ヒュベニの公式】緯度＋経度で座標指定。地球儀の球面上で，最短航路の距離[m]＋方角[θ]を算出するアルゴリズムのサンプルプログラム

地球上で「２点間の最短距離とその方角(東西南北)」を，緯度・経度から求める方法は？

曲面上で「三平方の定理」は使えない！

「ヒュベニの公式」サンプルコードと動作デモ

地球上（地球儀の球面上）で
実際に最短距離と方角を計算する
PHPのサンプルプログラムです。

ここでは
「ハワイ語たん」と
「パプア諸語たん」の
球面上での最短距離を算出しています。

<?php


// ヒュベニの公式で，地球上の２点間の最短距離とその方角を求めるコード
// ※日付変更線をまたぐ場合も正確に動作する
//
// @author:rwanda_go_tan




// ---------- 最短航路を計算する関数 ---------- 


// ２点間の最短距離をヒュベニの公式で求める。（メートル単位）
// 引数の緯度・経度は -180 ～ +180 の範囲内で渡すこと
// @author:rwanda_go_tan
function get_hubeny_distance( 

  $latitude1_deg, $longitude1_deg, // 座標1

  $latitude2_deg, $longitude2_deg  // 座標2

)
{
    // ---------- 緯度・経度の加工 ---------- 

    // 緯度・経度をラジアンに変換
    $latitude1_rad  = deg2rad( $latitude1_deg );
    $longitude1_rad = deg2rad( $longitude1_deg );
    $latitude2_rad  = deg2rad( $latitude2_deg );
    $longitude2_rad = deg2rad( $longitude2_deg );

    // 緯度の差分（ラジアン）
    $diff_latitude_rad = abs( $latitude1_rad - $latitude2_rad );

    // 経度の差算（ラジアン）：
    // 日付変更線をまたぐ場合にも
    // 必ず短いほうの弧長を計算するように
    $diff_longitude_rad = min([

      // 単純な差分で済む場合
      abs( $longitude1_rad - $longitude2_rad ),

      // 単純な差分だけでは長くなってしまう場合，
      // 残りの円弧が最短経路となる
      2*pi() - abs( $longitude1_rad - $longitude2_rad )

    ]);

    // 平均緯度（ラジアン）
    $average_latitude_rad = ( $latitude1_rad + $latitude2_rad ) / 2.0;


    // ---------- 各種定数の計算 ---------- 

    // 地球の長半径，赤道半径(Semi-major Axis of the Earth)
    $semi_major_axis = 6378137.0; 

    // 地球の短半径，極半径(Semi-minor Axis of the Earth)
    $semi_minor_axis = 6356752.314245; 

      // 定数の引用元：
      // https://www.trail-note.net/tech/calc_distance/
      // 世界的に使用されている WGS84 (1984年) を採用

    // 離心率を求める
    $E = sqrt(
      ( 
        // 長半径の2乗
        $semi_major_axis * $semi_major_axis 
        - 
        // 短半径の2乗
        $semi_minor_axis * $semi_minor_axis 
      ) 
      / 
      // 長半径の2乗
      ( $semi_major_axis * $semi_major_axis )
    );

    // 子午線(しごせん)曲率半径 M を求める
    $sin_lat = sin( $average_latitude_rad );
    $W = sqrt(
      1.0 
      - 
      $E * $E * ( $sin_lat * $sin_lat )
    );
    $M = $semi_major_axis * ( 1 - $E * $E ) / ( $W * $W * $W );

    // 卯酉線(ぼうゆうせん)曲率半径 N 
    $N = $semi_major_axis / $W; 


    // ---------- 距離算出 ---------- 

    $t1 = $M * $diff_latitude_rad;
    $t2 = $N * cos( $average_latitude_rad ) * $diff_longitude_rad;
    $result = sqrt(
      ( $t1 * $t1 ) 
      + 
      ( $t2 * $t2 )
    );


    // メートル単位で返却
    return $result;
}
// NOTE: Qiitaに載っているサンプルコードは
// ２点の最短航路が日付変更線をまたぐ場合に計算結果がおかしくなる。
// 上記コードではそのバグを修正してある。
// https://qiita.com/chiyoyo/items/b10bd3864f3ce5c56291

// ２点の位置情報から方位角を算出
// 引数の緯度・経度は -180 ～ +180 の範囲で渡す。
// 起点から見た終点へのラジアンを0～2πで返す。北が0で時計回りに角度が増える。
// @author:rwanda_go_tan
function get_azimuth_of_two_points(

  $latitude_from_deg, $longitude_from_deg, // 座標1，起点

  $latitude_to_deg, $longitude_to_deg  // 座標2，終点

){

  // Azimuth : 方位角
  // https://ja.wikipedia.org/wiki/%E6%96%B9%E4%BD%8D%E8%A7%92

  // xの差分
  $diff_longitude_deg = $longitude_to_deg - $longitude_from_deg;

  // NOTE:
  // 緯度・経度は360°単位であり，三角関数での計算時にはラジアンに変換

  // 方位角
  $theta = pi() / 2.0 
    - atan2(

      cos( $latitude_from_deg  / 180.0 * pi() )
      *
      tan( $latitude_to_deg    / 180.0 * pi() )

      -

      sin( $latitude_from_deg  / 180.0 * pi() )
      *
      cos( $diff_longitude_deg / 180.0 * pi() )
      ,

      sin( $diff_longitude_deg / 180.0 * pi() )

    )
  ;

  // 計算式：
  // https://keisan.casio.jp/exec/system/1257670779

  // PHP, atan2( y, x )
  // https://www.php.net/manual/ja/function.atan2.php

  // 角度を0から2πまでに正規化
  while( $theta >= 2*pi() ){ $theta -= 2*pi(); }
  while( $theta < 0       ){ $theta += 2*pi(); }

  return $theta;
}



// 2点の方位角を方角の名称に変換。
// 引数はラジアンで，0～2π。北が0で時計回りに角度が増える。
// @author:rwanda_go_tan
function translate_azimuth_to_direction_name( $azimuth_rad ) {

  // 方角リスト
  $direction_names = [
    "北", "北東", "東", "南東", 
    "南", "南西", "西", "北西", 
    "北" // π/16 だけずらす分を，この末尾要素で回収する
  ];

  // NOTE: 方位を8分割している。
  // つまり，- π/16 ～ + π/16 を「北」とみなす。

  // 方位角を配列のインデックスに変換
  $direction_index = floor(
    (
      // - π/16 を北とみなすために
      ( $azimuth_rad + ( pi() ) / 16.0 )
      /
      ( 2 * pi() )
    )
    *
    8.0 // 2πの回転を，全方位の個数のスキャンに対応させる
  );

  // 方角の名称
  $direction_name = $direction_names[ $direction_index ];

  return $direction_name;
}

// ---------- サンプル座標データ ---------- 


// ハワイ語たん
// https://gogakutan.jishotter.info/detail/hawaii_tan
$hawaigotan = [
  "lat" => 19.8967662,
  "lng" => -155.5827818 
];


// パプア諸語たん
// https://gogakutan.jishotter.info/detail/papuan_tan
$paputan = [
  "lat" => -6.314993,
  "lng" => 143.95555 
];




// ---------- 動作デモ ---------- 



// ヒュベニの公式で最短距離を出す
$distance_hybany = get_hubeny_distance(
  $hawaigotan["lat"], $hawaigotan["lng"], 
  $paputan["lat"],    $paputan["lng"]
);
echo "最短距離 = " . ( $distance_hybany / 1000.0 ) . " km.<br><br><br>";





// ２点の方位角を出す
$azimuth_rad = get_azimuth_of_two_points(
  $hawaigotan["lat"], $hawaigotan["lng"], 
  $paputan["lat"],    $paputan["lng"]
);
echo "方位角 = " . $azimuth_rad . " rad.<br><br><br>";





// ２点の方角名称を出す
$direction_name = translate_azimuth_to_direction_name( $azimuth_rad );
echo "方位名称 = " . $direction_name . " .<br><br><br>";



/*

上記サンプルの出力：

最短距離 = 7285.1919373465 km.
方位角 = 4.4095462906492 rad.
方位名称 = 南西 .

*/


?>

距離計算の用語：大円距離，測地線，方位角

地球上の「方角」は，直感とは違った答えになる！

「ヒュベニの式」世間のサンプルコードには，バグが潜んでいる場合もある！

正しく計算できるように，バグを除去しました。