title: Python calculate point of intersection of two great circles
body: <p>I am trying to calculate the point of intersection (lat and lon in degrees) of two great circles that are each defined by two points on the circle. I have been trying to follow method outlined <a href="https://stackoverflow.com/questions/2954337/great-circle-rhumb-line-intersection">here</a>.
  But the answer I get is incorrect, my code is below does anyone see where I went wrong?</p>
  
  <pre><code>import numpy as np
  from numpy import cross
  from math import cos, sin, atan2, asin, asinh
  
  ################################################
  #### Intersection of two great circles.
  # Points on great circle 1.
  glat1 = 54.8639587
  glon1 = -8.177818
  
  glat2 = 52.65297082
  glon2 = -10.78064876
  
  # Points on great circle 2.
  cglat1 = 51.5641564
  cglon1 = -9.2754284
  
  cglat2 = 53.35422063
  cglon2 = -12.5767799
  
  # 1. Put in polar coords.
  
  x1 = cos(glat1) * sin(glon1)
  y1 = cos(glat1) * cos(glon1)
  z1 = sin(glat1)
  
  x2 = cos(glat2) * sin(glon2)
  y2 = cos(glat2) * cos(glon2)
  z2 = sin(glat2)
  
  
  cx1 = cos(cglat1) * sin(cglon1)
  cy1 = cos(cglat1) * cos(cglon1)
  cz1 = sin(cglat1)
  
  cx2 = cos(cglat2) * sin(cglon2)
  cy2 = cos(cglat2) * cos(cglon2)
  cz2 = sin(cglat2)
  
  
  # 2. Get normal to planes containing great circles.
  #    It's the cross product of vector to each point from the origin.
  
  N1 = cross([x1, y1, z1], [x2, y2, z2])
  N2 = cross([cx1, cy1, cz1], [cx2, cy2, cz2])
  
  
  # 3. Find line of intersection between two planes.
  #    It is normal to the poles of each plane.
  
  L = cross(N1, N2)
  
  
  # 4. Find intersection points.
  
  X1 = L / abs(L)
  X2 = -X1
  
  
  ilat = asin(X1[2]) * 180./np.pi
  ilon = atan2(X1[1], X1[0]) * 180./np.pi
  </code></pre>
  
  <p>I should also mention this is on the Earth's surface (assuming a sphere).</p>