Ruby script to calculate a central angle betweetn 2 points on earth.

calc_central_angle.rb

#! /usr/local/bin/ruby
#---------------------------------------------------------------------------------
# 2地点の緯度・経度から中心角を計算
# * 2地点の緯度(Beta)／経度(Lambda)／楕円体高(Height)から
#   地球中心から2地点への線分のなす中心角を求める。
# * CalcCentralAngle クラスが実行クラス
# * CentralAngle クラスが計算の本質部分
#
#   Date          Author          Version
#   2019.08.29    mk-mode.com     1.00 新規作成
#
# Copyright(C) 2019 mk-mode.com All Rights Reserved.
#---------------------------------------------------------------------------------
# コマンドライン引数：
#   第1: 緯度-1
#   第2: 経度-1
#   第3: 緯度-2
#   第4: 経度-2
#   第5: 地球楕円体("GRS80"|"WGS84"|"BESSEL") (optional)
#   * 緯度は、東経が +, 西経が - で指定。
#   * 経度は、北緯が +, 南緯が - で指定。
#   * 2地点の高度(標高)は 0m に限定。
#   * 地球楕円体省略時は "GRS80" とみなす。
#---------------------------------------------------------------------------------
#
class CalcCentralAngle
  USAGE = "[USAGE] ./calc_central_angle.rb LAT_1 LNG_1 LAT_2 LNG_2 [ELLIPSOID]"

  # Initialization
  # * コマンドライン引数の取得
  def initialize
    if ARGV.size < 4
      puts USAGE
      exit
    end
    args = ARGV[0, 4].map(&:to_f)
    @pt_1 = args[0, 2]
    @pt_2 = args[2, 2]
    @el   = "GRS80"
    @el   = ARGV[4] if ARGV[4]
  end

  # Execution
  def exec
    puts "      POINT-1: %8.4f° %9.4f°" % @pt_1
    puts "      POINT-2: %8.4f° %9.4f°" % @pt_2
    puts "    ELLIPSOID: #{@el}"
    obj = CentralAngle.new(*((@pt_1 + @pt_2) << @el))
    ca = obj.calc_central_angle
    puts "CENTRAL-ANGLE: %8.4f°" % ca
  rescue => e
    msg = "[#{e.class}] #{e.message}\n"
    e.backtrace.each { |tr| msg << "\t#{tr}\n" }
    $stderr.puts msg
    exit 1
  end
end

class CentralAngle
  # 地球楕円体(長半径, 扁平率の逆数)
  ELLIPSOID = {
    "GRS80"  => [6_378_137.0  , 298.257_222_101],
    "WGS84"  => [6_378_137.0  , 298.257_223_563],
    "BESSEL" => [6_377_397.155, 299.152_813_000]
  }
  PI_180 = Math::PI / 180.0

  def initialize(lat_1, lng_1, lat_2, lng_2, ellipsoid="GRS80")
    @lat_1, @lng_1 = lat_1, lng_1
    @lat_2, @lng_2 = lat_2, lng_2
    el  = ELLIPSOID[ellipsoid]
    @a  = el[0]
    f   = 1 / el[1]
    @e2 = f * (2 - f)
  end

  # 中心角の計算
  # * 2ベクトルのなす角とみなして計算
  #
  # @return: 中心角 (Unit: °)
  def calc_central_angle
    x_1, y_1, z_1 = blh2ecef(@lat_1, @lng_1)
    x_2, y_2, z_2 = blh2ecef(@lat_2, @lng_2)
    return Math.acos(
      (x_1 * x_2 + y_1 * y_2 + z_1 * z_2) /
      (Math.sqrt(x_1 * x_1 + y_1 * y_1 + z_1 * z_1) *
       Math.sqrt(x_2 * x_2 + y_2 * y_2 + z_2 * z_2))
    ) / PI_180
  rescue => e
    raise
  end

private
  # BLH座標 -> ECEF直交座標 変換
  #
  # @param        lat: BLH  座標 緯度 (Unit: °)
  # @param        lng:           経度 (Unit: °)
  # @param         ht:           高度 (Unit: m)
  # @return [x, y, z]: ECEF 座標      (Unit: m)
def blh2ecef(lat, lng, ht=0.0)
    s = Math.sin(lat * PI_180)
    n = @a / Math.sqrt(1.0 - @e2 * s * s)
    x = (n + ht) \
      * Math.cos(lat * PI_180) \
      * Math.cos(lng * PI_180)
    y = (n + ht) \
      * Math.cos(lat * PI_180) \
      * Math.sin(lng * PI_180)
    z = (n * (1.0 - @e2) + ht) \
      * Math.sin(lat * PI_180)
    return [x, y, z]
  rescue => e
    raise
  end
end

CalcCentralAngle.new.exec if __FILE__ == $0