Internal: Routines to calculate great-circle distances on an ellipsoidal Earth
Anonym
Update Dec 14 2012: "map_2points.pro" is probably better for this.
Topic:
Date: Fri, 14 Feb 92 18:45:19 -0700
From: Dan O'Connell <geomagic@seismo.do.usbr.gov>
</geomagic@seismo.do.usbr.gov>
<geomagic@seismo.do.usbr.gov>Routines to calculate great-circle distances on an ellipsoidal Earth.</geomagic@seismo.do.usbr.gov>
<geomagic@seismo.do.usbr.gov>Discussion:</geomagic@seismo.do.usbr.gov>
Here are some idl routines to calculate great-circle distances on an ellipsoidal Earth. One gets the distance between two geographic lat-lon points and the other does the inverse.Solution:
<geomagic@seismo.do.usbr.gov>
PRO Latlon_disaz,eth,eph,th,phi,dist,az
;+
; NAME:
; LATLON_DISAZ
; PURPOSE:
; Return the great-circle distances (km) and azimuths (degrees) between
; the geographic lat-lon position eth-eph and the vector
; of geographic lat-lon positions th-phi, in dist-az.
;
; CALLING SEQUENCE:
; Latlon_disaz(eth,eph,th,phi,dist,az)
; INPUTS:
; eth = source latitude (signed decimal degrees)
; eph = source longitude (signed decimal degrees)
; th = vector of receiver latitudes (signed decimal degrees)
; phi = vector of receiver longitudes (signed decimal degrees)
; OUTPUTS:
; dist = double precision vector of source-receivers distances in km
; az = double precision vector of source-receivers azimuths in degrees
;
; COMMON BLOCKS:
; None.
; SIDE EFFECTS:
; None.
;
; PROCEDURE:
; The distances are calculated on an ellipsoidal earth using
; equidistant latitude following Brown (1984) Geophy. J. R. astrn. Soc.,
; v. 70, 445-459, "On the determination of source-receiver distances
; using a new equidistant latitude".
; MODIFICATION HISTORY:
; DRHO, BOR, Feb. 1992.
;-
;
on_error,2 ;Return to caller if an error occurs
n = n_elements(th) ;# of points.
if(n lt 1) then return
radco=1.745329251994329d-02 ;degree->radian
degco=57.29577951308232d0 ;radian->degree
flattn = 1.d0/298.257d0 ;earth flattening
onflat = 1.d0 - flattn
ellip0 = 1.001119d0 ;Brown (1984) equation (36) coefficients
ellip1 = 0.001687d0
degkm = 111.19504d0 ;km/degree for equidistant latitude
fltfac = onflat*sqrt(onflat) ;geographic -> equidistant flattening factor
pio2 = radco*90.d0
thk = eth*radco ;convert to radians
phk = eph*radco
;convert geographic latitude to equidistant latitude
if(abs(eth) ne 90.d0) then thg = atan(fltfac*tan(thk)) else thg = thk
thg = pio2 - thg ;convert to colatitude
d = sin(phk) ;now calculate spherical trig. quantities needed to
e = cos(phk) ;make spherical distance and azimuth calculations
f = sin(thg)
a = e*f
b = d*f
c = cos(thg)
g = c*e
h = c*d
e = -e
f = -f
az=dblarr(n) ;allocate return values
dist=dblarr(n)
for i=0,n-1,1 do begin ;calculate distances and azimuths for vector of
thc = th(i)*radco ;receiver points
phc = phi(i)*radco
if(abs(th(i)) ne 90.d0) then thg = atan(fltfac*tan(thc)) else thg = thc
thg = pio2 -thg ;convert to colatitude
d1 = sin(phc)
e1 = cos(phc)
f1 = sin(thg)
a1 = f1*e1
b1 = d1*f1
c1 = cos(thg)
xdeg = acos(a*a1+b*b1+c*c1) ;sherical triangle rule
ad = a1 - d
be = b1 - e
ag = a1 - g
bh = b1 - h
ck = c1 - f
az(i) = atan(ad*ad+be*be+c1*c1-2.d0,ag*ag+bh*bh+ck*ck-2.d0) ;Napier's rule
bh = sin(acos(-e*sin(az(i)))) ;Brown's great ellipse
xdeg = xdeg * degco * (ellip0 - ellip1*bh*bh) ;correction
az(i) = degco * az(i) ;convert to degrees
dist(i) = degkm * xdeg ;convert to km
endfor
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
PRO Disaz_latlon,eth,eph,dist,az,th,phi
;+
; NAME:
; DISAZ_LATLON
; PURPOSE:
; Calculates the lat-lon position of a point specified
; by the distance and great-circl azimuth from the point, ETH, EPH.
; The equidistant latitude method of Brown (1984) GJRAS is
; used to account for an ellipsoidal earth.
;
; CALLING SEQUENCE:
; Disaz_latlon(eth,eph,dist,az,th,phi)
; INPUTS:
; eth = source latitude (signed decimal degrees)
; eph = source longitude (signed decimal degrees)
; dist = vector of source-receivers distances in km
; az = vector of source-receivers azimuths in degrees
; OUTPUTS:
; th = double precision vector of receiver latitudes
; (signed decimal degrees)
; phi = double precision vector of receiver longitudes
; (signed decimal degrees)
;
; COMMON BLOCKS:
; None.
;
; SIDE EFFECTS:
; None.
;
; PROCEDURE:
; The distances are calculated on an ellipsoidal earth using
; equidistant latitude following Brown (1984) Geophy. J. R. astrn. Soc.,
; v. 70, 445-459, "On the determination of source-receiver distances
; using a new equidistant latitude".
; MODIFICATION HISTORY:
; DRHO, BOR, Feb. 1992.
;-
;
on_error,2 ;Return to caller if an error occurs
n = n_elements(dist) ;# of points.
if(n lt 1) then return
radco=1.745329251994329d-02 ;degree->radian
degco=57.29577951308232d0 ;radian->degree
flattn = 1.d0/298.257d0 ;earth flattening
onflat = 1.d0 - flattn
ellip0 = 1.001119d0 ;Brown (1984) equation (36) coefficients
ellip1 = 0.001687d0
degkm = 111.19504d0 ;km/degree for equidistant latitude
degkmi = 1.d0/degkm
fltfac = onflat*sqrt(onflat) ;geographic -> equidistant flattening factor
pio2 = radco*90.D0
thk = eth*radco ;convert to radians
phk = eph*radco
;convert source geographic latitude to equidistant latitude
if(abs(thk) ne 90.d0) then thg = atan(fltfac*tan(thk)) else thg = thk
thg = pio2 - thg ;convert to colatitude
sa = sin(thg) ;calculate source trig quantities
ca = cos(thg)
cp = cos(phk)
th=dblarr(n) ;allocate return values
phi=dblarr(n)
for i=0,n-1,1 do begin
azim = radco * az(i) ;convert azimuth to radians
;calculate quantities for great-ellipse correction (brown (36))
bh = sin(acos(cp*sin(azim)))
;convert distance to radians including great-ellipse correction
delta = (radco * dist(i)*degkmi)/(ellip0 - ellip1*bh*bh)
sd = sin(delta) ;calculate projection trig quantities
sz = sin(azim)
cz = cos(azim)
cd = cos(delta)
th(i) = acos(sa*sd*cz+cd*ca) ;calculate colatitude of new point
phi(i) = atan(sd*sz,sa*cd-sd*cz*ca) ;calculate latitude difference
;convert to longitude of new point in decimal degrees
phi(i) = (phi(i) * degco) + eph
if(phi(i) ge 360.d0) then phi(i) = phi(i) - 360.d0 ;for longitude between
if(phi(i) le -360.d0) then phi(i) = phi(i) + 360.d0 ; -360 & 360
;convert from equidistant latitude to geographic latitude
th(i) = pio2 - th(i) ;convert from colatitude to latitude
sd = th(i) * degco ; convert to latitude in decimal degrees
if(abs(sd) ne 90.d0) then th(i) = atan(tan(th(i))/fltfac)
th(i) = th(i) * degco
endfor
return
end
</geomagic@seismo.do.usbr.gov>