Computations involving spheres

This page describes routines from the JHU/APL/S1R IDL Library (called the S1R library below) that deal with computations involving spheres.

Several areas are covered in this page

Distance between two points.
Range and bearing from one point to another
Finding a point with a given range and bearing from a reference point.
Find the intersection point of two circles.
Some applications.

Return to the Spheres main page

Distance between two points

Two routines are of interest for this computation: sphdist and ll2rb. The first is a function and gives only the distance between the points. The second is a procedure and gives both distance (range) and azimuth (bearing). The built-in help for both is given below, followed by some examples.

IDL> d = sphdist(/help)
 Angular distance between points on a sphere.
 d = sphdist(long1, lat1, long2, lat2)
   long1 = longitude of point 1.         in
   lat1 = latitude of point 1.           in
   long2 = longitude of point 2.         in
   lat2 = latitude of point 2.           in
   d = angular distance between points.  out
 Keywords:
   /DEGREES means angles are in degrees, else radians.
 Notes: points 1 and 2 may be arrays.

IDL> ll2rb, /help
 From latitude, longitude compute range, bearing.
 ll2rb, lng0, lat0, lng1, lat1, dist, azi
   lng0, lat0 = long, lat of reference point (deg).    in
   lng1, lat1 = long, lat of point of interest (deg).  in
   dist = range to point point of interest (radians).  out
   azi = azimuth to point of interest (degrees).       out
 Notes: A unit sphere is assumed, thus dist is in radians
   so to get actual distance multiply dist by radius.
   Useful constants:
   Radius of Earth (mean) = 6371.23 km = 3958.899 miles.

Examples

Let lng1 = 0
    lat1 = 0
    lng2 = 0
    lat2 = 40

IDL> print,sphdist(lng1,lat1,lng2,lat2,/deg)
      40.0000

IDL> ll2rb,lng1,lat1,lng2,lat2,r,b & r=r*!radeg & help,r,b
R               FLOAT     =       40.0000

Note: in sphdist if /degrees is specified then both the input positions and output distance are in degrees. In ll2rb the input positions are always in degrees and the output distance is always in radians (sometimes more convenient).

Let lng1 = -90
    lat1 = 40
    lng2 = 35
    lat2 = 8.5

IDL> print,sphdist(lng1,lat1,lng2,lat2,/deg)
      109.849

IDL> ll2rb,lng1,lat1,lng2,lat2,r,b & r=r*!radeg & help,r,b
R               FLOAT     =       109.849

Range and bearing from one point to another

This computation is done using the routine ll2rb, same as used in the above section to compute the distance between two points (see the built-in help listed above).

Examples

Let lng1 = 0
    lat1 = 0
    lng2 = 0
    lat2 = 40

IDL> ll2rb,lng1,lat1,lng2,lat2,r,b & r=r*!radeg & help,r,b
R               FLOAT     =       40.0000
B               FLOAT     =       0.00000

Let lng1 = -90
    lat1 = 40
    lng2 = 35
    lat2 = 8.5

IDL> ll2rb,lng1,lat1,lng2,lat2,r,b & r=r*!radeg & help,r,b
R               FLOAT     =       109.849
B               FLOAT     =       59.4659

Finding a point with a given range and bearing from a reference point

This computation is done using the routine rb2ll. The built-in help is:

IDL> rb2ll, /help
 From range, bearing compute latitude, longitude .
 rb2ll, lng0, lat0, dist, azi, lng1, lat1
   lng0, lat0 = long, lat of starting point (deg).     in
   dist = range to point of interest in RADIANS.       in
   azi = azimuth to point of interest (degrees).       in
   lng1, lat1 = long, lat of point of interest (deg).  out
 Notes: A unit sphere is assumed, thus dist is in radians.
   Useful constants:
   Radius of Earth (mean) = 6371.23 km = 3958.899 miles.
   Distance to horizon from height H above surface:
     For small H: dist = sqrt(2*H/R) in Radians
     For large H: dist = acos(R/(R+H)) in Radians
   To plot horizon from lat0, lng0, H:
     rb2ll,lng0,lat0,dist,makex(0,360,1),plng,plat
     plots,plng-360.,plat,psym=3

Find the intersection point of two circles

sphic
sphgc

Some applications

Plotting great circle routines in any IDL map projection

The above routines make it easy to plot a great circle between two points. In the following example the range and bearing of the second point from the first are found using ll2rb. Then the latitude and longitude of points of increasing range with that bearing are computed using rb2ll. The results are then plotted.

lat1 = 40.
lng1 = -90.
lat2 = 8.5
lng2 = 35.

map_set,/cont                             ; Plot map. 
ll2rb,lng1,lat1,lng2,lat2,r,b             ; Find r & b of pt2.
rb2ll,lng1,lat1,maken(0,r,100),b,lng,lat  ; Step along track.
plots,lng,lat, col=3                      ; Plot track.
plots,lng1,lat1,psym=2, col=1             ; Plot point 1.
plots,lng2,lat2,psym=2, col=2             ; Plot point 2.

Circles of visibility

Let R be the radius of the earth. For an object of height H above the surface, the circle on the earth's surface where the object appears ALT degrees above the horizon may be found with the following algorithm.

   Know:
   R = Radius of earth = 6371.23 km.
   H = object height above surface.
   ALT = desired altitude of object above horizon (degrees).

   Want to find:
   D = distance along surface from point directly beneath the object.
   Use to find points on a circle of this radius.

   Algorithm:
   b = (ALT + 90)/!radeg
   D = !pi - b - asin(sin(b)*R/(R+H))

The chemical releases in space discussed in the CRRES example will be used in this example.

 DATE    TIME      LAT.     LONG.   ALT.
                            (WEST)  (KM)
JAN 18   0503      7.4      61.0    33565
JAN 25   0843     16.6 S    38.9     5799

First set up an array of bearings:

  a = makex(0,360,5)

Now find angle b in the above algorithm for ALT=0, 30, and 60 degrees:

  b = ([0,30,60]+90)/!radeg

The computation will be shown only for chemical release point 1, but both release points will be included in the accompanying figure. (the color table will have to be customized or the following will not show the altitude circles).

  lat1 = 7.4
  lng1 = -61.0
  r = 6371.23
  h = 33565.
  d = !pi - b - asin(sin(b)*r/(r+h))

  map_set, /cont
  plots, makex(-180,180,90),0
  plots, lng1, lat1, psym=2, col=1
  for i=0, 2 do begin  $
    rb2ll, lng1, lat1, d(i), a, lng, lat  &$
    plots, lng, lat, col=2  &$
  endfor

Where the symbols are plotted the chemical release would appear directly overhead. It would appear on the horizon (ignoring refraction) for the largest circles, 30 degrees above for the next circle inward, and 60 for the smallest circle. Some of the predictions stated that both the size and brightnesses would be like that of the full moon. The lithium released were expected to be bright red. The Persian Gulf War was in progress about the time of these releases and it was interesting to speculate what thoughts such resleases would trigger if seen in that area. However all releases were below the horizon for that part of the world.