title: Intersections between Geodesics (shortest distance paths) on the surface of a sphere
body: <p>I've searched far and wide but have yet to find a suitable answer to this problem.  Given two lines on a sphere, each defined by their start and end points, determine whether or not and where they intersect.  I've found this site (<a href="http://mathforum.org/library/drmath/view/62205.html" rel="nofollow">http://mathforum.org/library/drmath/view/62205.html</a>) which runs through a good algorithm for the intersections of two great circles, although I'm stuck on determining whether the given point lies along the finite section of the great circles.</p>
  
  <p>I've found several sites which claim they've implemented this, Including some questions here and on stackexchange, but they always seem to reduce back to the intersections of two great circles.</p>
  
  <p>The python class I'm writing is as follows and seems to almost work:</p>
  
  <pre><code>class Geodesic(Boundary):
    def _SecondaryInitialization(self):
      self.theta_1 = self.point1.theta
      self.theta_2 = self.point2.theta
      self.phi_1 = self.point1.phi
      self.phi_2 = self.point2.phi
  
      sines = math.sin(self.phi_1) * math.sin(self.phi_2)
      cosines = math.cos(self.phi_1) * math.cos(self.phi_2)
      self.d = math.acos(sines - cosines * math.cos(self.theta_2 - self.theta_1))
  
      self.x_1 = math.cos(self.theta_1) * math.cos(self.phi_1)
      self.x_2 = math.cos(self.theta_2) * math.cos(self.phi_2)
      self.y_1 = math.sin(self.theta_1) * math.cos(self.phi_1)
      self.y_2 = math.sin(self.theta_2) * math.cos(self.phi_2)
      self.z_1 = math.sin(self.phi_1)
      self.z_2 = math.sin(self.phi_2)
  
      self.theta_wraps = (self.theta_2 - self.theta_1 &gt; PI)
      self.phi_wraps = ((self.phi_1 &lt; self.GetParametrizedCoords(0.01).phi and
          self.phi_2 &lt; self.GetParametrizedCoords(0.99).phi) or (
          self.phi_1 &gt; self.GetParametrizedCoords(0.01).phi) and
          self.phi_2 &gt; self.GetParametrizedCoords(0.99))
  
    def Intersects(self, boundary):
      A = self.y_1 * self.z_2 - self.z_1 * self.y_2
      B = self.z_1 * self.x_2 - self.x_1 * self.z_2
      C = self.x_1 * self.y_2 - self.y_1 * self.x_2
      D = boundary.y_1 * boundary.z_2 - boundary.z_1 * boundary.y_2
      E = boundary.z_1 * boundary.x_2 - boundary.x_1 * boundary.z_2
      F = boundary.x_1 * boundary.y_2 - boundary.y_1 * boundary.x_2
  
      try:
        z = 1 / math.sqrt(((B * F - C * E) ** 2 / (A * E - B * D) ** 2)
            + ((A * F - C * D) ** 2 / (B * D - A * E) ** 2) + 1)
      except ZeroDivisionError:
        return self._DealWithZeroZ(A, B, C, D, E, F, boundary)
  
      x = ((B * F - C * E) / (A * E - B * D)) * z
      y = ((A * F - C * D) / (B * D - A * E)) * z
  
      theta = math.atan2(y, x)
      phi = math.atan2(z, math.sqrt(x ** 2 + y ** 2))
  
      if self._Contains(theta, phi):
        return point.SPoint(theta, phi)
  
      theta = (theta + 2* PI) % (2 * PI) - PI
      phi = -phi
  
      if self._Contains(theta, phi):
        return spoint.SPoint(theta, phi)
  
      return None
  
    def _Contains(self, theta, phi):
      contains_theta = False
      contains_phi = False
  
      if self.theta_wraps:
        contains_theta = theta &gt; self.theta_2 or theta &lt; self.theta_1
      else:
        contains_theta = theta &gt; self.theta_1 and theta &lt; self.theta_2
  
      phi_wrap_param = self._PhiWrapParam()
      if phi_wrap_param &lt;= 1.0 and phi_wrap_param &gt;= 0.0:
        extreme_phi = self.GetParametrizedCoords(phi_wrap_param).phi
        if extreme_phi &lt; self.phi_1:
          contains_phi = (phi &lt; max(self.phi_1, self.phi_2) and
              phi &gt; extreme_phi)
        else:
          contains_phi = (phi &gt; min(self.phi_1, self.phi_2) and
              phi &lt; extreme_phi)
      else:
        contains_phi = (phi &gt; min(self.phi_1, self.phi_2) and
            phi &lt; max(self.phi_1, self.phi_2))
  
      return contains_phi and contains_theta
  
    def _PhiWrapParam(self):
      a = math.sin(self.d)
      b = math.cos(self.d)
      c = math.sin(self.phi_2) / math.sin(self.phi_1)
      param = math.atan2(c - b, a) / self.d
  
      return param
  
    def _DealWithZeroZ(self, A, B, C, D, E, F, boundary):
      if (A - D) is 0:
        y = 0
        x = 1
      elif (E - B) is 0:
        y = 1
        x = 0
      else:
        y = 1 / math.sqrt(((E - B) / (A - D)) ** 2 + 1)
        x = ((E - B) / (A - D)) * y
  
      theta = (math.atan2(y, x) + PI) % (2 * PI) - PI
      return point.SPoint(theta, 0)
  
  def GetParametrizedCoords(self, param_value):
      A = math.sin((1 - param_value) * self.d) / math.sin(self.d)
      B = math.sin(param_value * self.d) / math.sin(self.d)
  
      x = A * math.cos(self.phi_1) * math.cos(self.theta_1) + (
      B * math.cos(self.phi_2) * math.cos(self.theta_2))
      y = A * math.cos(self.phi_1) * math.sin(self.theta_1) + (
          B * math.cos(self.phi_2) * math.sin(self.theta_2))
      z = A * math.sin(self.phi_1) + B * math.sin(self.phi_2)
  
      new_phi = math.atan2(z, math.sqrt(x**2 + y**2))
      new_theta = math.atan2(y, x)
  
      return point.SPoint(new_theta, new_phi)  
  </code></pre>
  
  <p>EDIT: I forgot to specify that if two curves are determined to intersect, I then need to have the point of intersection.</p>