Circle on a Sphere Equation

Date: 05/05/2000 at 00:51:18
From: Stephen Walker
Subject: Geocentric Trig

I've read several of your answers on finding the distance between two 
points using latitude and longitude. What I need to do is write a 
program that finds all points that fall within a given radius from a 
single latitude and longitude.

Is there a way to find a maximum latitude/longitude from a specific 
point so that all latitudes/longitudes fall within a specific radius?

I've tried everything I can think of and haven't come up with anything 
useful. I think this is an interesting problem.

Thanks for your help.
Steve

Date: 05/05/2000 at 08:19:29
From: Doctor Jerry
Subject: Re: Geocentric Trig

Hi Steve,

Because there are infinitely many points falling within a given great 
circle distance D of a given point, I assume that you mean a test as 
to whether an arbitrary point is or is not within D of the given 
point.

Because the distance D is expressed as D = alpha*a, where alpha, which 
is an angle with vertex at the center of the earth, subtends the great 
circle arc along which the distance D is measured, you can test 
whether alpha*a is less or equal to the given D. The expression for 
alpha will involve the latitudes and longitudes of the given point and 
the test point.

As to the maximum latitude/longitude, wouldn't the maximum latitude 
and minimum latitude points be on the meridian through the given 
point? And wouldn't the maximum longitude and minimum longitude points 
be on the parallel through the given point? I don't have a proof of 
this but it sure seems likely.

- Doctor Jerry, The Math Forum
  http://mathforum.org/dr.math/   

Date: 05/05/2000 at 09:12:27
From: Stephen Walker
Subject: Re: Geocentric Trig

Wow, what you said in your first paragraph is what I'm looking for. So 
if I understand this right, then I can do the following:

     alpha*a < D

where alpha is the angle between the two lines that form the 
coordinates, D is the distance I'm testing for, and alpha is some 
combination of the two points' latitude/longitude. Any idea what this 
combination is? Or what the complete formula would be, given the two 
points?

Alpha wouldn't be:

     acos[sin(alpha1)*sin(alpha2)
      + cos(alpha1)*cos(alpha2)*cos(beta1-beta2)] 

where alpha is the latitude of the 2 points in radians and beta is the 
longitude of the two points in radians, would it? Which then, I 
believe, gives me the arc distance in radians?

Wow my trig is really rusty. :) I think I'm completely off here.

I really appreciate the quick reply!

Thanks,
Steve

Date: 05/06/2000 at 08:31:34
From: Doctor Jerry
Subject: Re: Geocentric Trig

Hi Stephen,

Working in spherical coordinates a, p (for phi, the co-latitude), and 
t (for theta, related to longitude) for a moment, let (a_0,p_0,t_0) be 
the spherical coordinates of the center point (the point from which D 
would be measured), and (a,p,t) any other point. Then the rectangular 
coordinates of these points are:

     (a*sin(p_0)*cos(t_0),a*sin(p_0)*sin(t_0),a*cos(p_0))
and
     (a*sin(p_1)*cos(t_1),a*sin(p_1)*sin(t_1),a*cos(p_1))

Now think of these as position vectors from the center of the earth of 
the two points:

     v_0 = {a*sin(p_0)*cos(t_0),a*sin(p_0)*sin(t_0),a*cos(p_0)}
and
     v_1 = {a*sin(p_1)*cos(t_1),a*sin(p_1)*sin(t_1),a*cos(p_1)}

Using the properties of dot product, we see that

     v_0.v_1 = a*a*cos(alpha)

Writing out the dot product, removing the common factor of a^2 from 
both sides, and simplifying a bit with a trig identity,

     cos(alpha) = cos(p_0)*cos(p_1) + sin(p_0)*sin(p_1)*cos(t_1-t_0)

The distance D between (a,p_0,t_0) and (a,p_1,t_1) is:

     alpha*a = a*arccos[cos(p_0)*cos(p_1)+
                         sin(p_0)*sin(p_1)*cos(t_1-t_0)]

So, if (a,p_1,t_1) is to be less than or equal to D from (a,p_0,t_0), 
as measured on a great circle through these two points, we must have:

   a*arccos[cos(p_0)*cos(p_1)+sin(p_0)*sin(p_1)*cos(t_1-t_0)] <= D

This looks a lot like the formula you gave. Taking the simplest case, 
where 0 < p_0 < 90 and 0 < t_0 < 90, and not allowing (a,p_1,t_1) to 
over the pole or any kind of special case like that, we can convert 
from spherical to geographical angles. Letting L be latitude, we know 
that L = 90-p. Also, because in the above formula the longitudes only 
enter as a difference, I think we can take t_1-t_0 = long_1-long_0. 
So, the above formula becomes:

   a*arccos[sin(p_0)*sin(p_1)+cos(p_0)*cos(p_1)*cos(t_1-t_0)] <= D

This is, I think, your formula, in different notation. So, given a 
fixed center (a,p_0,t_0) and a test point (a,p_1,t_1) you can find out 
if it is within D or not be testing the above inequality.  

As I think you said, if in the above we replace <= by =, we have a 
formula defining a relation between p_1 and t_1, which would be a 
formula for the curve of the points exactly D from (a,p_0,t_0). 
Ignoring special cases again, I think you can solve the equality for 
t_1 in terms of p_1. Suppose this done and we find t_1 = f(p_1).

Then, parametric equations for the curve of equal distance would be

     x = a*sin(p_1)*cos(f(p_1))

     y = a*sin(p_1)*sin(f(p_1))

     z = a*cos(p_1)

I'm thinking about p_1 as a parameter. Because f would involve arccos, 
one would have to take some care to graph the entire curve.

I think your comments are correct. Inside a trig function like sine or 
cosine, the angles can be in radians or degrees. When you do an 
arccosine, however, you want the angle to be in radians for the 
formula s = a*b, where b is the central angle, a is the radius of the 
circle, and s is the arc length.

- Doctor Jerry, The Math Forum
  http://mathforum.org/dr.math/   

Date: 05/07/2000 at 10:38:38
From: Stephen Walker
Subject: Re: Geocentric Trig

Yeah, I got it to work using the formula that I sent back to you. 
Thanks so much for your help! If there's anything I can do for you, 
let me know.

Thanks again.
Steve