WSJT-X/GeoDist.f
Joe Taylor 9590ff5f43 More code cleanup: silenced many g95 warning messages about variables "set
but not used".


git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/trunk@275 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
2006-09-06 19:51:52 +00:00

103 lines
2.8 KiB
Fortran

subroutine geodist(Eplat, Eplon, Stlat, Stlon,
+ Az, Baz, Dist)
implicit none
real eplat, eplon, stlat, stlon, az, baz, dist
C JHT: In actual fact, I use the first two arguments for "My Location",
C the second two for "His location"; West longitude is positive.
c
c
c Taken directly from:
c Thomas, P.D., 1970, Spheroidal geodesics, reference systems,
c & local geometry, U.S. Naval Oceanographic Office SP-138,
c 165 pp.
c
c assumes North Latitude and East Longitude are positive
c
c EpLat, EpLon = End point Lat/Long
c Stlat, Stlon = Start point lat/long
c Az, BAz = direct & reverse azimuith
c Dist = Dist (km); Deg = central angle, discarded
c
real BOA, F, P1R, P2R, L1R, L2R, DLR, T1R, T2R, TM,
+ DTM, STM, CTM, SDTM,CDTM, KL, KK, SDLMR, L,
+ CD, DL, SD, T, U, V, D, X, E, Y, A, FF64, TDLPM,
+ HAPBR, HAMBR, A1M2, A2M1
real AL,BL,D2R,Pi2
data AL/6378206.4/ ! Clarke 1866 ellipsoid
data BL/6356583.8/
c real pi /3.14159265359/
data D2R/0.01745329251994/ ! degrees to radians conversion factor
data Pi2/6.28318530718/
BOA = BL/AL
F = 1.0 - BOA
c convert st/end pts to radians
P1R = Eplat * D2R
P2R = Stlat * D2R
L1R = Eplon * D2R
L2R = StLon * D2R
DLR = L2R - L1R ! DLR = Delta Long in Rads
T1R = ATan(BOA * Tan(P1R))
T2R = ATan(BOA * Tan(P2R))
TM = (T1R + T2R) / 2.0
DTM = (T2R - T1R) / 2.0
STM = Sin(TM)
CTM = Cos(TM)
SDTM = Sin(DTM)
CDTM = Cos(DTM)
KL = STM * CDTM
KK = SDTM * CTM
SDLMR = Sin(DLR/2.0)
L = SDTM * SDTM + SDLMR * SDLMR * (CDTM * CDTM - STM * STM)
CD = 1.0 - 2.0 * L
DL = ACos(CD)
SD = Sin(DL)
T = DL/SD
U = 2.0 * KL * KL / (1.0 - L)
V = 2.0 * KK * KK / L
D = 4.0 * T * T
X = U + V
E = -2.0 * CD
Y = U - V
A = -D * E
FF64 = F * F / 64.0
Dist = AL*SD*(T -(F/4.0)*(T*X-Y)+FF64*(X*(A+(T-(A+E)
+ /2.0)*X)+Y*(-2.0*D+E*Y)+D*X*Y))/1000.0
TDLPM = Tan((DLR+(-((E*(4.0-X)+2.0*Y)*((F/2.0)*T+FF64*
+ (32.0*T+(A-20.0*T)*X-2.0*(D+2.0)*Y))/4.0)*Tan(DLR)))/2.0)
HAPBR = ATan2(SDTM,(CTM*TDLPM))
HAMBR = Atan2(CDTM,(STM*TDLPM))
A1M2 = Pi2 + HAMBR - HAPBR
A2M1 = Pi2 - HAMBR - HAPBR
1 If ((A1M2 .ge. 0.0) .AND. (A1M2 .lt. Pi2)) GOTO 5
If (A1M2 .lt. Pi2) GOTO 4
A1M2 = A1M2 - Pi2
GOTO 1
4 A1M2 = A1M2 + Pi2
GOTO 1
c
c all of this gens the proper az, baz (forward and back azimuth)
c
5 If ((A2M1 .ge. 0.0) .AND. (A2M1 .lt. Pi2)) GOTO 9
If (A2M1 .lt. Pi2) GOTO 8
A2M1 = A2M1 - Pi2
GOTO 5
8 A2M1 = A2M1 + Pi2
GOTO 5
9 Az = A1M2 / D2R
BAZ = A2M1 / D2R
c
c Fix the mirrored coords here.
c
az = 360.0 - az
baz = 360.0 - baz
end