Finds the intersection point on a sphere between the line (x1,y1),(x2,y2) and the perpendicular line passing through (x3,y3), and compares the result with an approximated approach considering the sphere as an euclidean space

vector.rb

def rad(a)
  a * Math::PI / 180.0
end

def deg(a)
  a * 180.0 / Math::PI
end

class Vector

  attr_reader :x, :y, :z, :r

  def initialize(x, y, z)
    @x, @y, @z = x.to_f, y.to_f, z.to_f
    @r = Math.sqrt(x ** 2 + y ** 2 + z ** 2)
  end
  
  def self.from_spherical(r, azimuth, inclination)
    new(r * Math.sin(inclination) * Math.cos(azimuth),
        r * Math.sin(inclination) * Math.sin(azimuth),
        r * Math.cos(inclination))
  end
  
  def to_s
    cartesian_str
  end
  
  def cartesian_str
    "[#{@x}, #{@y}, #{@z}]"
  end
  
  def spherical_str(radians = true)
    if radians
      "[#{@r} <#{azimuth} <#{inclination}]"
    else
      "[#{@r} <#{deg(azimuth)} <#{deg(inclination)}]"
    end
  end
  
  def azimuth
    Math.atan2(@y, @x)
  end
  
  def inclination
    Math.acos(@z / @r)
  end
  
  def rotate_x(angle)
    Vector.new(@x,
               @y * Math.cos(angle) - @z * Math.sin(angle),
               @y * Math.sin(angle) + @z * Math.cos(angle))
  end
  
  def rotate_y(angle)
    Vector.new(@x * Math.cos(angle) + @z * Math.sin(angle),
               @y,
               - @x * Math.sin(angle) + @z * Math.cos(angle))
  end
  
  def rotate_z(angle)
    Vector.new(@x * Math.cos(angle) - @y * Math.sin(angle),
               @x * Math.sin(angle) + @y * Math.cos(angle),
               @z)
  end
  
  def distance_to(other)
    phi_2 = Math.sin((other.inclination - inclination) / 2) ** 2
    lambda_2 = Math.sin((other.azimuth - azimuth) / 2) ** 2
    x = phi_2 + Math.sin(other.inclination) * Math.sin(inclination) * lambda_2
    2 * Math.asin(Math.sqrt(x))
  end
  
end

if __FILE__ == $0
  #v1 = Vector.from_spherical(1, rad(30), rad(70))
  #v2 = Vector.from_spherical(1, rad(60), rad(80))
  #v3 = Vector.from_spherical(1, rad(30), rad(90))

  v1 = Vector.from_spherical(1, rad(-49.5), rad(90 - 25.2))
  v2 = Vector.from_spherical(1, rad(-49.0), rad(90 - 25.3))
  v3 = Vector.from_spherical(1, rad(-49.1), rad(90 - 25.1))
  
  #v1 = Vector.from_spherical(1, rad(-42), rad(90))
  #v2 = Vector.from_spherical(1, rad(-60), rad(90))
  #v3 = Vector.from_spherical(1, rad(-50), rad(90 - 1))
  
  gama, beta = -v1.azimuth, Math::PI / 2 - v1.inclination
  v1_a = v1.rotate_z(gama).rotate_y(beta)
  v2_a = v2.rotate_z(gama).rotate_y(beta)
  v3_a = v3.rotate_z(gama).rotate_y(beta)
  alpha = - Math.atan2(v2_a.z, v2_a.y)
  v1_b = v1_a.rotate_x(alpha)
  v2_b = v2_a.rotate_x(alpha)
  v3_b = v3_a.rotate_x(alpha)
  puts v1_b.spherical_str(false)
  puts v2_b.spherical_str(false)
  puts v3_b.spherical_str(false)
  vi = Vector.from_spherical(v3_b.r, v3_b.azimuth, Math::PI / 2)
  puts vi.spherical_str(false)

  puts v1_b.rotate_x(-alpha).rotate_y(-beta).rotate_z(-gama).spherical_str(false)
  puts v2_b.rotate_x(-alpha).rotate_y(-beta).rotate_z(-gama).spherical_str(false)
  puts v3_b.rotate_x(-alpha).rotate_y(-beta).rotate_z(-gama).spherical_str(false)
  print "\n"
  v4 = vi.rotate_x(-alpha).rotate_y(-beta).rotate_z(-gama)
  puts v4.spherical_str(false)
  
  x1, y1 = deg(v1.azimuth), deg(v1.inclination)
  x2, y2 = deg(v2.azimuth), deg(v2.inclination)
  x3, y3 = deg(v3.azimuth), deg(v3.inclination)
  dx12 = x1 - x2
  dy12 = y1 - y2
  k = (dx12 * (y3 - y1) - dy12 * (x3 - x1)) / (dx12 ** 2 + dy12 ** 2)
  x4 = x3 + k * dy12
  y4 = y3 - k * dx12
  print x4, " ", y4, "\n"
  dist_results = Vector.from_spherical(1, rad(x4), rad(y4)).distance_to(v4)
  dist_ref_points = v1.distance_to(v2)
  print dist_results, " ", dist_ref_points, "\n"
  print "error: ", 100 * (dist_results / dist_ref_points), "%"
end