緯度差1°あたりの経線弧の長さを求めるプログラム

• 高校地学の教科書を読んでいたら計算したくなったので。久しぶりのプログラミング楽しかったです。
• 計算結果は教科書に載っているものと微妙に異なる数値が出てきます。
• 楕円の弧長の近似値を求めるのに区分求積法を使いました。もっと良い方法がありそうです。
• 区分求積法を実装したついでにπの近似値も求めてみたり
  
• 調べてみると、緯度差1°あたりの経線弧の長さを求めるという問題からは、どうやら「楕円積分」という話題に発展できるみたいですね。
  • 楕円積分 〜 振り子の周期を求める [物理のかぎしっぽ]
  • 楕円積分 - Wikipedia

# 緯度差1°あたりの経線弧の長さを求める

include Math

A = 6378.137 # 赤道半径
B = 6356.752 # 極半径
N = 50000 # 近似の細かさ
EPSILON = 1.0/180*PI # 1°

def main
  p lng(deg2rad(90)) # 極
  p lng(deg2rad(66.0 + 20.0/60)) # 北フィンランド
  p lng(deg2rad(45.0)) # フランス
  p lng(deg2rad(1.0 + 31.0/60)) # エクアドル
  p lng(deg2rad(0)) # 赤道
end

def deg2rad(deg)
  deg/180.0*PI
end

def lng(alpha)
  theta_1 = alpha_to_theta(alpha-EPSILON/2)
  theta_2 = alpha_to_theta(alpha+EPSILON/2)
  arc(theta_1, theta_2, N)
end

# 楕円の二つの媒介変数の値からその間の弧の長さを求める
def arc(theta_1, theta_2, n)
  integrate(->theta { sqrt((A*sin(theta))**2 + (B*cos(theta))**2) }, theta_1, theta_2, n)
end

# 緯度αから楕円の媒介変数θの値を得る
def alpha_to_theta(alpha)
  #atan(B/A * tan(alpha))
  atan2(B*sin(alpha), A*cos(alpha))
end

# 区分求積法により定積分の近似値を計算する
def integrate(f, a, b, n)
  delta_x = (b-a).to_f / n
  sum = 0.0
  n.times do |i|
    x = a + delta_x * i
    sum += f.(x) * delta_x
  end
  sum
end

main() if $0 == __FILE__