=begin
世界測地系WGS84の定数等をまとめたファイル
=end

class WGS84

    R_E             = 6378137.0        # 赤道半径[m]
    F_E             = (1.0 / 298.257223563)
    OMEGA_EARTH     = 7.292115E-5      # 地球自転速度
    OMEGA_EARTH_IAU = 7.2921151467E-5
    MU_EARTH        = 3.986005E14      # 地球重力定数[m^3/s^2]
    EPSILON_EARTH   = 0.0818191908426  # 第一偏心性
    G_WGS0          = 9.7803267714     # 赤道上重力
    G_WGS1          = 0.00193185138639 # 重力公式定数

    def WGS84.r_meridian(latitude) # 子午線(南北)方向の曲率半径
      return R_E * (1.0 - (EPSILON_EARTH ** 2)) \
                  / ((1.0 - (EPSILON_EARTH ** 2) * (Math.sin(latitude)) ** 2) ** 1.5);
    end
    
    def WGS84.r_normal(latitude) # 卯酉線(東西)方向の曲率半径
      return R_E / Math.sqrt(1.0 - (EPSILON_EARTH ** 2) * (Math.sin(latitude) ** 2));  
    end
    
    def WGS84.gravity(latitude)
      return G_WGS0 * (1.0 + G_WGS1 * (Math.sin(latitude) ** 2)) \
                          / Math.sqrt(1.0 - (EPSILON_EARTH ** 2) * (Math.sin(latitude) ** 2)); 
    end
    
    # 2点間の距離を求める、緯度経度はラジアンで指定のこと
    # 計算にはヒュベニ公式を使う
    # http://yamadarake.web.fc2.com/trdi/2009/report000001.html
    def WGS84.distance(lat1, lng1, lat2, lng2)
      
      # 平均緯度
      lat_avg = (lat1 + lat2) / 2

      # 緯度差
      lat_delta = lat1 - lat2
      
      # 経度差
      lng_delta = lng1 - lng2
      
      lat_s = Math::sin(lat_avg)
      lat_c = Math::cos(lat_avg)
      
      w = Math::sqrt(1.0 - (EPSILON_EARTH * EPSILON_EARTH) * (lat_s * lat_s))
      
      # 子午線曲率半径
      r_M = 6335439.0 / (w ** 3)
      
      # 卯酉線曲率半径
      r_N = R_E / w
      
      t1 = r_M * lat_delta
      t2 = r_N * lat_c * lng_delta
      
      return Math::sqrt((t1 * t1) + (t2 * t2))
    end
end