Find a point on the sphere minimizing the average distance to a given set of points, within a given degree of precision.

earth_midpoint.rb

class Point
  attr_reader :lat, :long

  def initialize(lat, long)
    @lat = lat.to_f
    @long = long.to_f
    
    if (-Math::PI/2 > @lat || @lat > Math::PI/2)
      raise "lat must between -PI/2 and PI/2, inclusive"
    end
    if (0 > @long || @long >= Math::PI*2)
      raise "long must be between 0 and 2*PI (including 0 but not 2*PI)"
    end
  end

  def distance(point)
    Math.acos(Math.sin(lat)*Math.sin(point.lat) + Math.cos(lat)*Math.cos(point.lat)*Math.cos((long-point.long).abs))
  end
end

###########################

class Algorithm
  def initialize(points, precision)
    @points = points
    @precision = precision
  end

  def compute
    distance = Float::INFINITY
    lat = -Math::PI/2
    long = 0
    point = nil
    while (lat <= Math::PI/2) do
      while (long < 2*Math::PI) do
        new_point = Point.new(lat, long)
        new_distance = (@points.map{|pt| pt.distance(new_point)}).inject(:+)
        if new_distance < distance
          point = new_point
          distance = new_distance
        end
        long = long + @precision
      end
      long = 0
      lat = lat + @precision
    end
    point
  end
end

###########################

point_1 = Point.new(-Math::PI/2, 0)
point_2 = Point.new(0, 0)
point_3 = Point.new(0, Math::PI/2)
point_4 = Point.new(-1.5, 1)
point_5 = Point.new(0.1, 4)
point_6 = Point.new(1.6, 6)
point_7 = Point.new(-0.4, 2)

points = [point_1, point_2, point_3, point_4, point_5, point_6, point_7]
precision = 1.0/360

# This will give you the latitude and longitude of "the" point minimizing
# the average distance to the 7 points above, good to within one 360th
# of a radian.  This means it tests about 2.5 million points.
p Algorithm.new(points, precision).compute