SUBROUTINE SPOINTS
A (FLAT1, ELON1, FLAT2, ELON2, N, FLAT, ELON)
C S/R to determine the latitude and longitude of N equally spaced
C points between the points (FLAT1, ELON1) and (FLAT2, ELON2).
C All values specified in degrees. End points are excluded from
C the set of returned points.
C Written by G. IRLAM 8-Apr-87
C Last modified 23-mar-88
C INPUT DATA FLAT1 North latitude of first end point
C ELON1 East longitude of first end point
C FLAT2 North latitude of second end point
C ELON2 East longitude of second end point
C N No. of intermediate points
C OUTPUT DATA FLAT N values of latitude of intermediate points
C FLON N values of longitude of intermediate points
DIMENSION FLAT (N), ELON (N)
CALL DISTBEAR (FLAT1, ELON1, FLAT2, ELON2, DIST, BEAR)
D = DIST / (N + 1)
DO 100 I = 1, N
CALL LATLON (FLAT1, ELON1, I * D, BEAR, FLAT (I), ELON (I))
100 continue
RETURN
END
SUBROUTINE DISTBEAR
A (X1, Y1, X2, Y2, D, B)
C S/R to determine the distance and bearing from point
C (X1, Y1) to point (X2, Y2). All values specified in degrees.
C All distances in kilometres.
C Written by G. IRLAM 8-Apr-87
c Last modified 15-june-87 david barker
c formulae courtesy of: "Fundamental Formulas of Physics"
c ed: Donald H. Menzel p. 8, section 2.11
c and considering the spherical triangle formed by
c the two given vertices and the North Pole.
C INPUT DATA X1 North latitude of first end point
C Y1 East longitude of first end point
C X2 North latitude of second end point
C Y2 East longitude of second end point
C OUTPUT DATA D Shortest distance between points
C B Bearing of second point from first point
C along along shortest path. Degrees East
C of North , 0 < B <= 360.
PI = 4.0 * ATAN (1.0)
RAD = PI / 180
DEG = 180 / PI
R0 = 6372.1
c modified 16-oct-87 to agree with QPAR
r0 = 6370.
YY = Y2 - Y1
IF (YY .LE. -180.0) YY = YY + 360.0
IF (YY .GT. 180.0) YY = YY - 360.0
SX1 = SIN (X1 * RAD)
CX1 = COS (X1 * RAD)
SX2 = SIN (X2 * RAD)
CX2 = COS (X2 * RAD)
CYY = COS (YY * RAD)
CD = SX1 * SX2 + CX1 * CX2 * CYY
IF (ABS (CD) .GT. 1.0) CD = SIGN (1.0, CD)
D = R0 * ACOS (CD)
SD = SQRT (1.0 - CD * CD)
c modified 23-mar-88 to account for co-located points, or point 1
c being N or S pole - in this case, bearing is set to zero
if(cx1*sd.eq.0.) then
b = 0.
else
CBN = (SX2 - SX1 * CD) / (CX1 * SD)
if (abs(cbn) .gt. 1.0) cbn = sign(1.0, cbn)
bn = acos(cbn) *deg
if (yy .lt. 0.0) bn = 360. -bn
b = bn
endif
RETURN
END
SUBROUTINE LATLON
A (X1, Y1, D, B, X2, Y2)
C S/R to determine the point (X2, Y2) at a distance D and bearing
C B from the point (X1, Y1).
C All values specified in degrees. All distances in kilometres.
C Written by G. IRLAM 8-Apr-87
c Last modified 15-june-87 david barker
c Formulae courtesy of: "Fundamental Formulas of Physics."
c ed: Donald H. Menzel p. 8 section 2.11
c and considering the spherical triangle formed by
c the bearing point, the North Pole, and the side
c described by bearing B and distance D.
C INPUT DATA X1 North latitude of bearing point
C Y1 East longitude of bearing point
C D Shortest distance to new points
C B Bearing of new point from bearing point
C along along shortest path. Degrees East
C of North, 0 < B <= 360.
C OUTPUT DATA X2 North latitude of new point
C Y2 East longitude of new point
PI = 4.0 * ATAN (1.0)
RAD = PI / 180
DEG = 180 / PI
R0 = 6372.1
c modified 16-oct-87 to agree with QPAR
c modified 7-jun-88 to have longitudes between 0 and 359.999999
r0 = 6370.
BN = B
SX1 = SIN (X1 * RAD)
CX1 = COS (X1 * RAD)
SD = SIN (D / R0)
CD = COS (D / R0)
CBN = COS (BN * RAD)
sx2 = sx1*cd + cx1*sd*cbn
IF (ABS (sX2) .GT. 1.0) sX2 = SIGN (1.0, sX2)
x2 = asin(sx2)*deg
CX2 = SQRT (1.0 - sX2 * sX2)
CYY = (CD - SX1 * SX2) / (CX1 * CX2)
IF (ABS (CYY) .GT. 1.0) CYY = SIGN (1.0, CYY)
YY = ACOS (CYY) * DEG
if (bn .gt. 180.0) yy = -yy
y2 = y1 + yy
if (y2 .lt. 0.0) y2 = y2 + 360.0
RETURN
END