module Geocoder::Calculations

def bearing_between(point1, point2, options = {})
  # set default options
  options[:method] ||= Geocoder.config.distances
  options[:method] = :linear unless options[:method] == :spherical
  # convert to coordinate arrays
  point1 = extract_coordinates(point1)
  point2 = extract_coordinates(point2)
  # convert degrees to radians
  point1 = to_radians(point1)
  point2 = to_radians(point2)
  # compute deltas
  dlat = point2[0] - point1[0]
  dlon = point2[1] - point1[1]
  case options[:method]
  when :linear
    y = dlon
    x = dlat
  when :spherical
    y = Math.sin(dlon) * Math.cos(point2[0])
    x = Math.cos(point1[0]) * Math.sin(point2[0]) -
        Math.sin(point1[0]) * Math.cos(point2[0]) * Math.cos(dlon)
  end
  bearing = Math.atan2(x,y)
  # Answer is in radians counterclockwise from due east.
  # Convert to degrees clockwise from due north:
  (90 - to_degrees(bearing) + 360) % 360
end

def bounding_box(point, radius, options = {})
  lat,lon = extract_coordinates(point)
  radius  = radius.to_f
  units   = options[:units] || Geocoder.config.units
  [
    lat - (radius / latitude_degree_distance(units)),
    lon - (radius / longitude_degree_distance(lat, units)),
    lat + (radius / latitude_degree_distance(units)),
    lon + (radius / longitude_degree_distance(lat, units))
  ]
end

def compass_point(bearing, points = COMPASS_POINTS)
  seg_size = 360 / points.size
  points[((bearing + (seg_size / 2)) % 360) / seg_size]
end

def coordinates_present?(*args)
  args.each do |a|
    # note that Float::NAN != Float::NAN
    # still, this could probably be improved:
    return false if (!a.is_a?(Numeric) or a.to_s == "NaN")
  end
  true
end

def distance_between(point1, point2, options = {})
  # set default options
  options[:units] ||= Geocoder.config.units
  # convert to coordinate arrays
  point1 = extract_coordinates(point1)
  point2 = extract_coordinates(point2)
  # convert degrees to radians
  point1 = to_radians(point1)
  point2 = to_radians(point2)
  # compute deltas
  dlat = point2[0] - point1[0]
  dlon = point2[1] - point1[1]
  a = (Math.sin(dlat / 2))**2 + Math.cos(point1[0]) *
      (Math.sin(dlon / 2))**2 * Math.cos(point2[0])
  c = 2 * Math.atan2( Math.sqrt(a), Math.sqrt(1-a))
  c * earth_radius(options[:units])
end

def distance_to_radians(distance, units = nil)
  units ||= Geocoder.config.units
  distance.to_f / earth_radius(units)
end

def earth_radius(units = nil)
  units ||= Geocoder.config.units
  case units
    when :km; EARTH_RADIUS
    when :mi; to_miles(EARTH_RADIUS)
    when :nm; to_nautical_miles(EARTH_RADIUS)
  end
end

def endpoint(start, heading, distance, options = {})
  options[:units] ||= Geocoder.config.units
  radius = earth_radius(options[:units])
  start = extract_coordinates(start)
  # convert degrees to radians
  start = to_radians(start)
  lat = start[0]
  lon = start[1]
  heading = to_radians(heading)
  distance = distance.to_f
  end_lat = Math.asin(Math.sin(lat)*Math.cos(distance/radius) +
                Math.cos(lat)*Math.sin(distance/radius)*Math.cos(heading))
  end_lon = lon+Math.atan2(Math.sin(heading)*Math.sin(distance/radius)*Math.cos(lat),
                Math.cos(distance/radius)-Math.sin(lat)*Math.sin(end_lat))
  to_degrees [end_lat, end_lon]
end

def extract_coordinates(point)
  case point
  when Array
    if point.size == 2
      lat, lon = point
      if !lat.nil? && lat.respond_to?(:to_f) and
        !lon.nil? && lon.respond_to?(:to_f)
      then
        return [ lat.to_f, lon.to_f ]
      end
    end
  when String
    point = Geocoder.coordinates(point) and return point
  else
    if point.respond_to?(:to_coordinates)
      if Array === array = point.to_coordinates
        return extract_coordinates(array)
      end
    end
  end
  [ NAN, NAN ]
end

def geographic_center(points)
  # convert objects to [lat,lon] arrays and convert degrees to radians
  coords = points.map{ |p| to_radians(extract_coordinates(p)) }
  # convert to Cartesian coordinates
  x = []; y = []; z = []
  coords.each do |p|
    x << Math.cos(p[0]) * Math.cos(p[1])
    y << Math.cos(p[0]) * Math.sin(p[1])
    z << Math.sin(p[0])
  end
  # compute average coordinate values
  xa, ya, za = [x,y,z].map do |c|
    c.inject(0){ |tot,i| tot += i } / c.size.to_f
  end
  # convert back to latitude/longitude
  lon = Math.atan2(ya, xa)
  hyp = Math.sqrt(xa**2 + ya**2)
  lat = Math.atan2(za, hyp)
  # return answer in degrees
  to_degrees [lat, lon]
end

def km_in_mi
  KM_IN_MI
end

def km_in_nm
  KM_IN_NM
end

def latitude_degree_distance(units = nil)
  units ||= Geocoder.config.units
  2 * Math::PI * earth_radius(units) / 360
end

def longitude_degree_distance(latitude, units = nil)
  units ||= Geocoder.config.units
  latitude_degree_distance(units) * Math.cos(to_radians(latitude))
end

def mi_in_km
  1.0 / KM_IN_MI
end

def nm_in_km
  1.0 / KM_IN_NM
end

def radians_to_distance(radians, units = nil)
  units ||= Geocoder.config.units
  radians * earth_radius(units)
end

def random_point_near(center, radius, options = {})
  # set default options
  options[:units] ||= Geocoder.config.units
  # convert to coordinate arrays
  center = extract_coordinates(center)
  earth_circumference = 2 * Math::PI * earth_radius(options[:units])
  max_degree_delta =  360.0 * (radius / earth_circumference)
  # random bearing in radians
  theta = 2 * Math::PI * rand
  # random radius, use the square root to ensure a uniform
  # distribution of points over the circle
  r = Math.sqrt(rand) * max_degree_delta
  delta_lat, delta_long = [r * Math.cos(theta), r * Math.sin(theta)]
  [center[0] + delta_lat, center[1] + delta_long]
end

def to_degrees(*args)
  args = args.first if args.first.is_a?(Array)
  if args.size == 1
    (args.first * 180.0) / Math::PI
  else
    args.map{ |i| to_degrees(i) }
  end
end

def to_kilometers(mi)
  mi * mi_in_km
end

def to_miles(km)
  km * km_in_mi
end

def to_nautical_miles(km)
  km * km_in_nm
end

def to_radians(*args)
  args = args.first if args.first.is_a?(Array)
  if args.size == 1
    args.first * (Math::PI / 180)
  else
    args.map{ |i| to_radians(i) }
  end
end