How Far Do the Birds Fly?

Date: 08/12/2002 at 17:16:24
From: John Vance
Subject: Radius of the Earth

Dear Dr Math,

I am a racing pigeon enthusiast, and one problem I wanted to solve for 
myself and others in the sport is how to calculate the survey distance 
our birds fly from point A to point B, where A is the release point 
and B is the home loft of a returning bird.

We have the latitude and longitude of each home loft and the release 
points that we race from, and up until now have used a surveying 
service to calculate these distances for us. With the advent of the 
Internet and Web sites, more and more information is available to the 
curious and I have gone on a quest to create an area on my site where 
other flyers can figure out the distance between two points.

From your site, I have been able to program a Web page to do the 
calculations, and that is when I came to the problem I now face.

The formula I use when written in expressions my server software can 
understand is as follows:

mylength = ACOS((COS(a1)*COS(b1)*COS(a2)*COS(b2)) + (COS(a1)*SIN(b1)
*COS(a2)*SIN(b2)) + (SIN(a1)*SIN(a2)))*r 

the value "r" however is the great unknown.

I live at 48.8236666667 -122.5888333333

Since the radius of the earth is ever-changing from equator to pole, 
and since I am somewhat near the halfway mark in Latitude, I chose to 
use the volumetric radius of the earth as given by NASA (3958.76 
miles) since I thought that this radius would most closely 
approximate my local conditions.

When I read this entry in the Dr. Math archives:

   Using Longitude and Latitude to Determine Distance
   http://mathforum.com/library/drmath/view/51711.html 

Where Nat Keller uses a value for "r" of 3963.1 miles it only further 
points out that few measurements would agree, since the likelihood 
that two calculations from two different solutions would use the same 
radius values is small.

We all agree on the formula, just not the appropriate radius values.  
I think there is a 22-mile difference between the radius of the earth 
at the equator and at the poles. Is there a formula that derives a 
close approximation of the average radius between two latitude 
readings? For example, if the two latitude readings I was using were 
48.8236666667 and 44.9386544, is there a formula I could use to 
calculate an average value for "r"?

Thank you,
John Vance

Date: 08/13/2002 at 09:42:57
From: Doctor Rick
Subject: Re: Radius of the Earth

Hi, John.

Recently I have been learning a little about working with ellipsoids -
a closer approximation to the shape of the earth, with its slight 
polar flattening. Maybe I've only learned enough to be dangerous, but 
here goes!

I think the correct radius to use is the radius of curvature of the 
earth in the direction of interest, at the point of interest. Over a 
long course the radius of curvature would change, so this can be only 
an approximation. However, we're only talking about a maximum error 
of +-11/3960 = 0.28% or so even when we pretend the earth is a 
sphere, so we aren't going to get too large an error by pretending 
it's an ellipsoid but the radius of curvature doesn't change, unless 
the course is, say, 1000 miles or more.

The radius of curvature at a point on the ellipsoid depends on the 
direction of travel. Two radii are defined: M, the radius of curvature 
in the meridian (north-south), and N, the radius of curvature in the 
prime vertical (east-west). They are

  M = a(1 - e^2)/(1 - e^2*sin^2(phi))^3/2

  N = a/sqrt(1 - e^2*sin^2(phi))

where a is the radius of the earth at the equator, e is the 
eccentricity of the earth (e = ), and phi is the latitude. 
(Specifically it's the geodetic latitude; there are several ways to 
define latitude on an ellipsoid, but this one is what we measure 
astronomically.) The radius of curvature along a course at bearing 
beta is expressed in terms of the two principal radii of curvature:

  R = M*N/(N*cos^2(beta)+M*sin^2(beta))

The mean radius at a given latitude - that is, averaging over all 
directions at a given point - is

  sqrt(MN) = a*sqrt(1-e^2)/(1 - e^2*sin^2(phi))

I got these formulas from the following site:

   Geodesy Course (The Hong Kong Polytechnic University)
   http://www.lsgi.polyu.edu.hk/Cyber-Class/geodesy/index.htm  

From the same Web site, the values of a and e for the formulas I gave 
are:

  a = 6378137 m = 3963.191 miles
  e = 0.08181919084

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

Date: 08/13/2002 at 12:15:43
From: John Vance
Subject: Radius of the Earth

Thank you for the help.

sqrt(MN) = a*sqrt(1-e^2)/(1 - e^2*sin^2(phi))

I have tested this formula in an Excel spreadsheet and only get a 
correct answer at phi=0 degrees, and then it is the polar radius 
6356752.314 (shouldn't the polar radius be at phi=90 degrees?).

I believe that I am not processing the sin(phi) properly.  Below are 
my values.

a=6378137 meters
e^2=0.00669438
1-e^2=0.99330562
sqrt(1-e^2)=0.996647189
a*sqrt(1-e^2)=6356752.314

When phi=45 degrees the below values are what I get:

sin(phi)=0.707106781
sin^2(phi)=0.5 or sin(phi)*sin(phi)=0.5
1-e^2*sin^2(phi)=0.99665281

If I set phi to 90 degrees, I get Sqrt(MN)=6399593.626 or 21456.62576 
meters in excess of the equatorial radius.

What am I doing wrong?

Thanks.

Date: 08/13/2002 at 13:53:50
From: Doctor Rick
Subject: Re: Radius of the Earth

Hi, John.

I think your only error is in your expectations. Remember that we are 
calculating the radius of curvature, and this has nothing necessarily 
to do with the distance from the center of the earth.

The earth is flattened at the poles. The flatter a curve (or a surface 
containing a curve), the greater its radius of curvature - the radius 
of the circle that most closely approximates the curve. So the radius 
of curvature at the poles will be greater than the equatorial radius, 
and indeed we find that it is a^2/b. It's the same in all directions, 
due to the symmetry at the poles, so the mean radius of curvature is 
also a^2/b.

At the equator, it's easy to see that the radius of curvature in an 
east-west direction is a, because the equator *is* a circle of this 
radius - no approximation is required. In the north-south direction 
we'll get a smaller radius of curvature; somehow it works out that 
the mean is exactly b.

The geodetic latitude is the complement of the angle between a line 
perpendicular to the surface and the polar axis. This perpendicular 
line does not necessarily go through the center of the earth. If a 
north-south path on the ellipsoid is short enough, it can be 
approximated by a circle of radius equal to the radius of curvature. 
Thus the arc length is approximately the radius of curvature times 
the difference in geodetic latitudes of the endpoints. This is as 
much as I understand so far regarding the calculation of arc lengths 
on an ellipsoid, but I think it is all you need.

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

Date: 08/13/2002 at 18:10:40
From: John Vance
Subject: Radius of the Earth

Thanks so much for the help.

My quick fix solutition to the problem is as follows:

a=equatorial radius
b=polar radius
f=(a-b)/a (flattening)

R(phi)=a-(a*f*sin(phi))

at phi=90 degrees
sin(90)=1
R(phi)=a-(a*f)
R(phi)=a-(a*(a-b)/a)
R(phi)=a-(a-b)
R(phi)=b

at phi=0 degrees
sin(0)=0
R(phi)=a-(a*f*sin(phi))
R(phi)=a

at phi=45
sin(45)=0.707106781
R(phi)=a-(a*f*sin(phi))
R(phi)=6363015.744 meters or 3953.79 statute miles

Whereas, (a+b)/2 = 6367444.657 meters or 3956.55 staute miles
which is the result if phi were equal to 30 degrees.

Not the perfect model but closer than using "a".

Now if I was calculating the distance of two points, one at 43 
degrees latitude and the other at 48 degrees latitude, I could 
approximate the average radius over the distance as:
phi=Lat(ave)=(43+48)/2
phi=45.5
and plug this value into:
R(phi)=a-(a*f*sin(phi))

One of several questions remaining would be if radians would be a 
better appoximation for phi.  I am not really all that good at trig 
to know when radians apply.

Your comments are welcome as to this solution as generalized as it is.

Thank you,

John Vance

Date: 08/14/2002 at 08:30:34
From: Doctor Rick
Subject: Re: Radius of the Earth

Hi, John.

I take it that you don't want to make use of the radius of curvature 
at all, and instead you are going with the approach of interpolating 
between a and b. Perhaps the following web page will convince you 
that the radius of curvature is really the way to go:

   GIS FAQ Question 5.1
   http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1 

This page has a table of radius of curvature versus latitude (just 
over halfway down the page; r = mean radius of curvature, R' = my M). 
You will see there that the radius of curvature is GREATEST at the 
poles and LEAST at the equator, contrary to your calculation based on 
distance from the center of the earth.

In reference to radians, it's not a matter of which gives a better 
approximation for phi; that would be like saying that meters give a 
better approximation for a distance than feet. Radians are just a 
different unit for measuring angles: 1 radian = 180/pi degrees, or 
57.295779 degrees. The angle MUST be in radians, however, before you 
multiply it by the radius (of curvature) to get an arc length. If you 
use the angle in degrees, you won't just have a worse approximation, 
you'll be off by a factor of 57.295779 - not even close!

I hope this helps.

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

Date: 08/14/2002 at 14:54:25
From: John Vance
Subject: Radius of the Earth

Hello Dr Math,

My assumption was that delta R (change in the length of the radius 
from the surface to the center of the earth) would be greater moving 
away from the equator then moving away from the pole and that is why 
I used a sin function.  For example, delta R from 90 degrees to 89 
degrees would be less than from 0 degrees to 1 degree.

If I have this backward, then possibly I should be using the cosine.

You may wonder why I even care about this.

We race our birds from distances of 75 miles out to 500 and 600 miles.  
Each bird carries a band with an RFID chip that holds a unique ID 
number. Each loft has an antenna that reads the chip upon the bird's 
return and logs the ID number and the time of arrival to 1/10th of a 
second. Through GPS we know the Lat and Long of the release point and 
each loft. We calculate the distance, divide that by the time of 
flight and derive a velocity for each bird in the race. The velocity 
is printed on our race results to 1/1000th of a ypm. Some lofts are 
70 - 80 miles farther than other lofts in the race so the better we 
are able to calculate the straight line distance the better we are 
able to declare the race winner based on the ypm criteria.

It is true that the birds do not travel in a straight line, but we 
must calculate distances as though they do.

As it is now, we pay for a surveyor to calculate our distances from 
each loft to each release point. I want to develop a Web site that 
would do the calculations for free. Also, we don't know what value the 
surveyor is using for r so I must be able to justify why my value 
for r is as accurate as is possible under the circumstances.

Sorry it has taken so much of your time to get this far. I have 
learned a lot from all your input. Thank you!

One more thing, and this is possibly why I was confused about the 
Radii of Curvature: 

latitude...........r...................R'..................N..........
 
00 degrees . 6357 km (3950 mi) . 6336 km (3937 mi) . 6378 km (3963 mi) 
15 degrees . 6360 km (3952 mi) . 6340 km (3940 mi) . 6379 km (3964 mi) 
30 degrees . 6367 km (3957 mi) . 6352 km (3947 mi) . 6383 km (3966 mi) 
45 degrees . 6378 km (3963 mi) . 6367 km (3957 mi) . 6389 km (3970 mi) 
60 degrees . 6388 km (3970 mi) . 6383 km (3966 mi) . 6394 km (3973 mi) 
75 degrees . 6396 km (3974 mi) . 6395 km (3974 mi) . 6398 km (3975 mi) 
90 degrees . 6399 km (3976 mi) . 6399 km (3976 mi) . 6399 km (3976 mi) 

Would this mean that if I am calculating distance near Lat 45, I 
should assume the radius of the earth is 6378 km, and if near the 
equator that the radius is 6399 km?

Thanks again for all the time you have spent,

John V

Date: 08/15/2002 at 13:42:19
From: Doctor Rick
Subject: Re: Radius of the Earth

Hi, John.

If you don't want to include the bearing in the calculation, then the 
average radius of curvature is the best you can do; that means the 
figures you quote. I gave you a formula to calculate the average 
radius of curvature of the ellipsoid exactly, rather than choosing an 
arbitrary means of interpolation as you have been doing. 

You see from the table, though, that there is still a variation in 
the actual radius of curvature at lat 45 degrees of plus or minus 
11 km from the average radius of curvature. That is, if the course is 
north-south, the radius of curvature is 11 km less than the average, 
and if the course is east-west, the radius of curvature is 11 km 
greater than the average. 

You'll be better off if you can use the formula I gave for the radius 
of curvature as a function of both latitude and bearing. If you can 
calculate the bearing between the two points (and I can show you 
formulas for this if you don't have one), you can then use the 
formula to get a more accurate radius for use in the distance formula.

In response to your previous comments about the need for this 
precision: I don't know the speed of a pigeon in a 600-mile race, but 
I'd estimate that you're measuring the time to a relative precision 
of about one millionth (10^-6) and printing the velocity to a relative 
precision of a millionth or even better. For the velocity to have an 
accuracy of this magnitude, you'd need to calculate the distance to an 
accuracy of plus or minus a yard or so. I doubt that you can get 
anywhere near that accuracy. Over 600 miles the radius of curvature 
could change by about 4 miles if it's a north-south course, and I 
would expect the sphere-based distance formula to show errors from 
other sources as well over such a distance. I suspect you'd need a 
much more complicated numerical algorithm to get the accuracy you 
seek. This goes beyond my basic knowledge of geodesy and ellipsoids.

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

Date: 08/21/2002 at 14:11:21
From: John Vance
Subject: Radius of the Earth

I suppose at this point I am incredibly buried in the terminology of 
this solution.

If I have the Lat and Lon of the release point and the home loft, 
what exactly should I do to arrive at the most accurate distance 
calculation?

For example, if the release point is at (I have converted to decimal 
notation):

Lat: 43.801216 Lon: -123.048701

and the home loft at:

Lat: 48.823666 Lon: -122.588833

What would I do according to the instructions you have been giving me?  
I must admit that I think I do not understand how to properly identify 
the vars and plug them into a formula to give the proper distance 
based on the Lat - Lon - direction of flight.

This would probably clear up my misunderstanding if I saw how you 
calculated the results.

Thank you,
John Vance

Date: 08/21/2002 at 15:15:15
From: Doctor Rick
Subject: Re: Radius of the Earth

Hi, John.

Here is one way to get an approximation. First find the heading from 
point 1, using this formula (from the Aviation Formulary on the Web):

  heading = mod(atan2(-sin(lon1-lon2)*cos(lat2), 
                      cos(lat1)*sin(lat2)
                      -sin(lat1)*cos(lat2)*cos(lon1-lon2)),
                2*pi)

I get 0.060276182 radians, or 3.453570825 degrees east of north. (If 
your atan2 function takes arguments (x,y) like Excel rather that 
(y,x) like C, then switch the arguments.)

Now we can use this heading in the calculation of the radius of 
curvature:

  M = a(1 - e^2)/(1 - e^2*sin^2(phi))^3/2
  N = a/sqrt(1 - e^2*sin^2(phi))
  R = M*N/(N*cos^2(beta)+M*sin^2(beta))

I get:

  M = 3937.942347 miles
  N = 3963.621319 miles
  R = 3938.03493 miles

Because the route is nearly north-south, the radius along the route is 
close to the radius of curvature in the meridian (north-south).

You have the formula for the distance; use the radius R.

If you want even better accuracy, you can calculate the latitude and 
longitude of the midpoint of the route, using the formula here:

   Latitude and Longitude of a Point Halfway between Two Points
   http://mathforum.org/library/drmath/view/51822.html 

Then calculate the heading at the midpoint, using the formula I used 
above, and use the latitude and heading at the midpoint in the 
calculation of R. It's still an approximation, but it's about as good 
as I could do.

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