In the Specification for the Wide Area Augmentation System (WAAS) on page 99 there are two formulae for calculating the "pierce point" of a GPS signal through the ionosphere. This may not be your application area, but may be able to adapt these to your application and they happen to work with longitude, latitude, azimuth, and elevation angle. Here's a picture from the document:
On the site http://paulbourke.net/geometry/circlesphere/ there is a section called "Intersection of a Line and a Sphere (or circle)" with more general equations for intersections. At the end of that section you will find some implementations in C, VB, Python, and LISP.
As a bonus (this is a limited time offer, you are one of a select few to receive this offer ;-) ) here is a short Python script for visualizing line-sphere intersections in Blender (adapted from this blog post. Note that is takes Earth-centered, Earth-fixed coordinates, not longitude/latitude/azimuth/elevation angles.
EARTH_EQUATOR = 6378137.0 # meters
IONO_HEIGHT = 350000.0 # meters
# scale down the radius by 1E6 so that it fits in blender's world
radius = (EARTH_EQUATOR + IONO_HEIGHT) / 1E6
from mathutils import Vector
w = 10
def MakePolyLine(objname, curvename, cList):
curvedata = bpy.data.curves.new(name=curvename, type='CURVE')
curvedata.dimensions = '3D'
objectdata = bpy.data.objects.new(objname, curvedata)
objectdata.location = (0,0,0) #object origin
polyline = curvedata.splines.new('POLY')
for num in range(len(cList)):
polyline.points[num].co = (cList[num])+(w,)
polyline.order_u = len(polyline.points)-1
polyline.use_endpoint_u = True
# here's a sample vector
listOfVectors = [(0.81,3.443,-5.723),(3.723,0.111,-5.603)]
MakePolyLine("NameOfMyCurveObject", "NameOfMyCurve", listOfVectors)
With more vectors, this would output something like the image below: