Calculates distance in meters between 2 GPS coordinates (Android port)

android_distance_port.rb

# Calculates distance in meters between 2 GPS coordinates
# Ported to Ruby from the Android SDK source
# Original paper explaning the process: http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf

# Arguments p1, p2 are instances of the class Point below:
# class Point
#   attr_accessor :lat, :lon
# end

def android_distance_port(p1, p2)
  lat1 = p1.lat
  lon1 = p1.lon
  lat2 = p2.lat
  lon2 = p2.lon

  max_iterations = 20
  lat1 *= Math::PI / 180.0
  lat2 *= Math::PI / 180.0
  lon1 *= Math::PI / 180.0
  lon2 *= Math::PI / 180.0

  a = 6378137.0
  b = 6356752.3142
  f = (a - b) / a
  a_sq_minus_b_sq_over_b_sq = (a * a - b * b) / (b * b)

  cap_l = lon2 - lon1
  cap_a = 0.0
  cap_u1 = Math.atan((1.0 - f) * Math.tan(lat1))
  cap_u2 = Math.atan((1.0 - f) * Math.tan(lat2))

  coscap_u1 = Math.cos(cap_u1)
  coscap_u2 = Math.cos(cap_u2)
  sincap_u1 = Math.sin(cap_u1)
  sincap_u2 = Math.sin(cap_u2)
  coscap_u1coscap_u2 = coscap_u1 * coscap_u2
  sincap_u1sincap_u2 = sincap_u1 * sincap_u2

  sigma = 0.0
  delta_sigma = 0.0
  cos_sqcap_alpha = 0.0
  cos_2_s_m = 0.0
  cos_sigma = 0.0
  sin_sigma = 0.0
  coscap_lambda = 0.0
  sincap_lambda = 0.0

  lambda = cap_l
  (0..max_iterations).each do |iter|
    lambda_orig = lambda
    coscap_lambda = Math.cos(lambda)
    sincap_lambda = Math.sin(lambda)
    t1 = coscap_u2 * sincap_lambda
    t2 = coscap_u1 * sincap_u2 - sincap_u1 * coscap_u2 * coscap_lambda
    sin_sq_sigma = t1 * t1 + t2 * t2
    sin_sigma = Math.sqrt(sin_sq_sigma)
    cos_sigma = sincap_u1sincap_u2 + coscap_u1coscap_u2 * coscap_lambda
    sigma = Math.atan2(sin_sigma, cos_sigma)
    sincap_alpha = 0.0
    unless sin_sigma == 0
      sincap_alpha = coscap_u1coscap_u2 * sincap_lambda / sin_sigma
    end
    cos_sqcap_alpha = 1.0 - sincap_alpha * sincap_alpha
    cos_2_s_m = 0.0
    unless cos_sqcap_alpha == 0
      cos_2_s_m = cos_sigma - 2.0 * sincap_u1sincap_u2 / cos_sqcap_alpha
    end

    u_squared = cos_sqcap_alpha * a_sq_minus_b_sq_over_b_sq
    cap_a = 1 + (u_squared / 16384.0) *
      (4096.0 + u_squared *
      (-768 + u_squared * (320.0 - 175.0 * u_squared)))
    cap_b = (u_squared / 1024.0) *
      (256.0 + u_squared *
      (-128.0 + u_squared * (74.0 - 47.0 * u_squared)))
    cap_c = (f / 16.0) *
      cos_sqcap_alpha *
      (4.0 + f * (4.0 - 3.0 * cos_sqcap_alpha))
    cos_2_s_m_sq = cos_2_s_m * cos_2_s_m
    delta_sigma = cap_b * sin_sigma *
      (cos_2_s_m + (cap_b / 4.0) *
      (cos_sigma * (-1.0 + 2.0 * cos_2_s_m_sq) -
      (cap_b / 6.0) * cos_2_s_m *
      (-3.0 + 4.0 * sin_sigma * sin_sigma) *
      (-3.0 + 4.0 * cos_2_s_m_sq)))

    lambda = cap_l +
      (1.0 - cap_c) * f * sincap_alpha *
      (sigma + cap_c * sin_sigma *
      (cos_2_s_m + cap_c * cos_sigma *
      (-1.0 + 2.0 * cos_2_s_m * cos_2_s_m)))

    delta = (lambda - lambda_orig) / lambda
    if (delta).abs < 1.0e-12
      break
    end
  end

  (b * cap_a * (sigma - delta_sigma)).to_f
end