Can anyone help me what's the algorythm to calculate points of intersection whet is given one point p1 [lat,long] , bearing x (let say 137 degrees), arc center [lat,long] and radius? I woulde like to get p2[lat,long] and p3[lat,long].

For better imagination what I mean something like on this picture but on sphere with only one point and bearing:

enter image description here


1 Answer 1


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:

Ionosphere pierce point

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.

import bpy

EARTH_EQUATOR = 6378137.0 # meters
IONO_HEIGHT = 350000.0 # meters

# scale down the radius by 1E6 so that it fits in blender's world


from mathutils import Vector

# weight    
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:

Ionosphere intersections

  • I'm glad to hear it!
    – katahdin
    Commented Mar 7, 2013 at 15:48
  • @user1275513: In addition to accepting an answer, you can also upvote helpful answers (such as this one) by clicking the upward pointing arrow tip next to the answer.
    – Jake
    Commented Mar 7, 2013 at 16:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.