[OSRS-PROJ] Donald Elliptic Projection for AT&T V&H Coordinates?
Duncan Agnew
agnew at bilby.ucsd.edu
Sun Dec 1 15:16:56 EST 2002
Being interested in puzzles, I worked through the code (more or less)
and rewrote it in a language I prefer. I append the result. I've added
some explanatory comments, although some mysteries remain--the original
Donald notes would be needed to clear these up. Apparently John Snyder
had a copy--anyone know where it might now be?
I agree with a comment that, for something like this, the procdure
is what it is, correct or not--it is purely conventional. So here
are some results, with what comparison I could make.
lat long V H Comparison
37 42 14.69 82 39 15.27 6363.235 2250.700 (one end)
41 02 55.53 112 03 39.35 7519.455 7107.410 (other end)
39 04 41.68 77 00 04.58 5587.002 1600.999 5587 1601 (1)
39 06 09.26 76 50 59.97 5568.001 1582.999 5568 1583 (1)
38 53 50.30 77 01 36.47 5621.999 1582.002 5622 1582 (1)
30 16 28.82 97 44 25.19 9003.475 3996.792 9004 3995 (2)
(1) - htlt.com V&H converter (examples on their Web page)
(2) - Peter H. Dana's coordinate systems Web page example
note that it appears to be conventional to give V and H as integers
It is possible that Dana did not have an "official" version--or, perhaps
his is correct, and the htlt one is not, but matches the code here
(they may have a common source). If there is an official answer, it
would be with Telcordia, which is part of what used to be Bell Labs.
Not a bad recreation for a holiday weekend.
Duncan Agnew
------------------
subroutine tovh(rlat,rlong,v,h)
implicit double precision(a-h,o-z)
save ifrst,dr,rotc,rots,rk9,rk10
c
c Transforms from geographic lat and long to the V[ertical] and H[orizontal]
c coordinates used for setting tariffs for wired telephony, originally by AT&T.
c
c Projection for USA and southern Canada developed by J. K. Donald of Bell
c Labs in 1957, known as the "Donald Elliptical Projection".
c
c This program written by D. C Agnew, 1 December 2002, based on C code
c embedded in the BBN Technologies openmap software. Original C code
c source was:
c
c Article: 54928 of comp.dcom.telecom
c Date: Tue, 29 Aug 1995 00:16:38 -0800
c From: stu at shell.portal.com (Stu Jeffery)
c Subject: Re: V&H Questions
c Message-ID: <telecom15.362.11 at eecs.nwu.edu>
c X-Telecom-Digest: Volume 15, Issue 362, Message 11 of 11
c
c This is an implementation of the Donald Elliptical Projection,
c a Two-Point Equidistant projection developed by Jay K. Donald
c of AT&T in 1956 to establish long-distance telephone rates.
c (ref: V-H Coordinate Rediscovered, Eric K. Grimmelmann, Bell
c Labs Tech. Memo, 9/80. (References Jay Donald notes of Jan 17, 1957.))
c Ashok Ingle of Bellcore also wrote an internal memo on the subject.
c
c Derived from a program obtained from an anonymous author
c within Bellcore by way of the National Exchange Carrier
c Association. Cleaned up and improved a bit by
c Tom Libert (tom at comsol.com, libert at citi.umich.edu).
c
c A partial description of the math involved is given in comments provided
c for the inverse code (same openmap source) by Col. G. L. Sicherman, with
c some corrections and modifications:
c
c Subject: V & H to Latitude and Longitude
c From: sicherman at lucent.com (Col. G.L. Sicherman)
c Date: 1997/03/05
c Message-ID: <telecom17.57.10 at massis.lcs.mit.edu>
c Newsgroups: comp.dcom.telecom
c
c Since I work for Bell Labs, I had no trouble getting a copy of
c Erik Grimmelmann's legendary memorandum. (Don't get your hopes
c up - Bell Labs has no intention of releasing it to the public!)
c
c V&H is a system of coordinates (V and H) for describing
c locations of rate centers in the United States. The
c projection, devised by J. K. Donald, is an elliptical,
c or "doubly equidistant" projection, scaled down by a factor
c of 0.003 to balance errors.
c
c The foci of the projection, from which distances are
c measured accurately (except for the scale correction),
c are at 37d 42m 14.69s N, 82d 39m 15.27s W (in Floyd Co.,
c Ky.) and 41d 02m 55.53s N, 112d 03m 39.35 W (in Webster
c Co., Utah). They are just 0.4 radians apart.
c
c Here is the transformation from latitude and longitude to V&H:
c First project the earth from its ellipsoidal surface
c to a sphere. This alters the latitude.
c Also subtract 52 degrees from the longitude.
c
c Compute the arc distances of the given point to the reference
c points, and transform them to the coordinate
c system in which the line through the reference points is the
c X-axis and the origin is the eastern reference point.
c Reduce by three-tenths of a percent, rotate by 76.597497
c degrees, and add 6363.235 to V and 2250.7 to H.
c
data rk1/0.99435487d0/,rk2/0.00336523d0/,rk3/-0.00065596d0/,
1 rk4/0.00005606d0/,rk5/-0.00000188d0/
c
c The above are the constants for the latitude transformation.
c John Snyder's description of this (p 235 of his _Flattening the
c Earth_, drawing on the original Donald manuscript) states
c "The Clarke 1866 ellipsoid was first transformed
c onto an authalic sphere, but the latitude adjustment was altered
c empirically to maintian true distance from the equator along the
c meridians".
c
c
c direction cosines of the E and W points (after transforming to the
c sphere and subtracting 52 degrees from the longitude)
c
data ex/0.40426992d0/,ey/0.68210848d0/,ez/0.60933887d0/
data wx/0.65517646d0/,wy/0.37733790d0/,wz/0.65449210d0/
c
c direction cosines of the perpendicular to the line between the
c points
c
data px/-0.555977821730048699d0/,py/-0.345728488161089920d0/,
1 pz/0.755883902605524030d0/
c
c rotation and false v and h
c
data rot/76.597497064d0/,transv/6363.235d0/,transh/2250.700d0/
c
c scale is set by Earth radius, which is supposed to be in statute miles,
c diminished by 0.3 percent. For the Clarke 1866 ellipsoid, using the
c equatorial radius and converting to statute miles using 39.37 inches/m,
c sqrt(10) times this number is 12495.3, so it is not clear how the value
c given here was actually derived. (The sqrt(10) comes from the distances
c being computed as sqrt(((v1-v2)**2+(h1-h2)**2)/10). )
c
data radius/12481.103d0/
c
c set a few more precomputed constants
c
data pi/3.14159265358979d0/
data ifrst/0/
if(ifrst.eq.0) then
ifrst=1
dr=180.d0/pi
rotc = dcos(rot/dr)
rots = dsin(rot/dr)
rk9 = radius*rotc
rk10 = radius*rots
endif
c
c convert to radians
c
rlat0=rlat/dr
rlon1=(rlong-52.d0)/dr
c
c transform the latitude
c
sqla = rlat0*rlat0
rlat1 = rlat0*(rk1 +(rk2 +(rk3 +(rk4 + rk5*sqla)*sqla)*sqla)*sqla)
c
c find direction cosines of point
c
clat = dcos(rlat1)
x = clat*dsin(rlon1)
y = clat*dcos(rlon1)
z = dsin(rlat1)
c
c e and w are the cosine of the angular distance (radians) between our
c point and the east and west centers.
c
e = ex*x + ey*y + ez*z
w = wx*x + wy*y + wz*z
if(e.lt.0.d0) e=0.d0
if(w.lt.0.d0) w=0.d0
if(e.gt.1.d0) e=1.d0
if(w.gt.1.d0) w=1.d0
c
c convert e and w to arc distance in radians
c
e=dacos(e)
w=dacos(w)
c original method using atan function
c if(e<1) e = pi/2 - atan(e/sqrt(1 - e*e))
c if(w<1) w = pi/2 - atan(w/sqrt(1 - w*w))
c if(e>=1) e=0
c if(w>=1) w=0
c
c find v and h in radians; the .16 and .8 are 0.4**2 and 2*0.4
c
ht = (e*e - w*w + 0.16d0)/0.8d0
vt = dsqrt(dabs(e*e - ht*ht))
c
c adjust the sign of v if need be
c
vs = px*x + py*y + pz*z
if(vs.lt.0) vt=-vt
c
c rotate and translate to get final v and h.
c
v = transv + rk9*ht - rk10*vt
h = transh + rk10*ht + rk9*vt
return
end
----------------------------------------
PROJ.4 Discussion List
See http://www.remotesensing.org/proj for subscription, unsubscription
and other information.
More information about the Proj
mailing list