SUBROUTINE W3FB10(DLAT1, DLON1, DLAT2, DLON2, BEARD, GCDKM)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C . . . .
C SUBPROGRAM: W3FB10 LAT/LONG PAIR TO COMPASS BEARING, GCD
C PRGMMR: CHASE ORG: NMC421 DATE:88-10-26
C
C ABSTRACT: GIVEN A PAIR OF POINTS (1) AND (2) GIVEN BY LATITUDE AND
C LONGITUDE, W3FB10 COMPUTES THE BEARING AND GREAT CIRCLE DISTANCE
C FROM POINT (1) TO POINT (2) ASSUMING A SPHERICAL EARTH. THE
C NORTH AND SOUTH POLES ARE SPECIAL CASES. IF LATITUDE OF POINT
C (1) IS WITHIN 1E-10 DEGREES OF THE NORTH POLE, BEARING IS THE
C NEGATIVE LONGITUDE OF POINT (2) BY CONVENTION. IF LATITUDE OF
C POINT (1) IS WITHIN 1E-10 DEGREES OF THE SOUTH POLE, BEARING IS
C THE LONGITUDE OF POINT (2) BY CONVENTION. IF POINT (2) IS WITHIN
C 1E-6 RADIANS OF THE ANTIPODE OF POINT (1), THE BEARING WILL BE
C SET TO ZERO. IF POINT (1) AND POINT (2) ARE WITHIN 1E-10 RADIANS
C OF EACH OTHER, BOTH BEARING AND DISTANCE WILL BE SET TO ZERO.
C
C PROGRAM HISTORY LOG:
C 88-08-29 CHASE, P.
C 88-09-23 CHASE, P. FIX DUMB SOUTH POLE ERROR
C 88-10-05 CHASE, P. FIX BEARING AMBIGUITY
C 90-04-12 R.E.JONES CONVERT TO CFT77 FORTRAN
C
C USAGE: CALL W3FB10(DLAT1, DLON1, DLAT2, DLON2, BEARD, GCDKM)
C
C INPUT ARGUMENT LIST:
C DLAT1 - REAL LATITUDE OF POINT (1) IN DEGREES NORTH.
C DLON1 - REAL LONGITUDE OF POINT (1) IN DEGREES EAST.
C DLAT2 - REAL LATITUDE OF POINT (2) IN DEGREES NORTH.
C DLON2 - REAL LONGITUDE OF POINT (2) IN DEGREES EAST.
C
C OUTPUT ARGUMENT LIST:
C BEARD - REAL BEARING OF POINT (2) FROM POINT (1) IN
C COMPASS DEGREES WITH NORTH = 0.0, VALUES FROM
C -180.0 TO +180.0 DEGREES.
C GCDKM - REAL GREAT CIRCLE DISTANCE FROM POINT (1) TO
C POINT (2) IN KILOMETERS.
C
C SUBPROGRAMS CALLED: NONE
C
C REMARKS: ACCORDING TO THE NMC HANDBOOK, THE EARTH'S RADIUS IS
C 6371.2 KILOMETERS. THIS IS WHAT WE USE, EVEN THOUGH THE VALUE
C RECOMMENDED BY THE SMITHSONIAN METEOROLOGICAL HANDBOOK IS
C 6371.221 KM. (I WOULDN'T WANT YOU TO THINK THAT I DIDN'T KNOW
C WHAT THE CORRECT VALUE WAS.)
C METHOD: THE POLES ARE SPECIAL CASES, AND HANDLED SEPARATELY.
C OTHERWISE, FROM SPHERICAL TRIGONOMETRY, THE LAW OF COSINES IS USED
C TO CALCULATE THE THIRD SIDE OF THE SPHERICAL TRIANGLE HAVING
C SIDES FROM THE POLE TO POINTS (1) AND (2) (THE COLATITUDES).
C THEN THE LAW OF SINES IS USED TO CALCULATE THE ANGLE AT POINT
C (1). A TEST IS APPLIED TO SEE WHETHER THE ARCSINE RESULT MAY BE
C BE USED AS SUCH, GIVING AN ACUTE ANGLE AS THE BEARING, OR WHETHER
C THE ARCSINE RESULT SHOULD BE SUBTRACTED FROM PI, GIVING AN OBTUSE
C ANGLE AS THE BEARING. THIS TEST IS DERIVED BY CONSTRUCTING A
C RIGHT SPHERICAL TRIANGLE USING THE POLE, POINT (2), AND THE
C MERIDIAN THROUGH POINT(1). THE LATITUDE OF THE RIGHT-ANGLED
C VERTEX THEN PROVIDES A TEST--IF LATITUDE (1) IS GREATER THAN THIS
C LATITUDE, THE BEARING ANGLE MUST BE OBTUSE, OTHERWISE ACUTE.
C IF THE TWO POINTS ARE WITHIN 1E-6 RADIANS OF EACH OTHER
C A FLAT EARTH IS ASSUMED, AND THE FOUR-QUADRANT ARCTANGENT
C FUNCTION IS USED TO FIND THE BEARING. THE Y-DISPLACEMENT IS
C THE DIFFERENCE IN LATITUDE AND THE X-DISPLACEMENT IS THE
C DIFFERENCE IN LONGITUDE TIMES COSINE LATITUDE, BOTH IN RADIANS.
C DISTANCE IS THEN THE DIAGONAL.
C FUNDAMENTAL TRIGONOMETRIC IDENTITIES ARE USED FREELY, SUCH
C AS THAT COS(X) = SIN(PI/2 - X), ETC. SEE ALMOST ANY MATHEMATICAL
C HANDBOOK, SUCH AS THE C.R.C. STANDARD MATH TABLES UNDER 'RELATIONS
C IN ANY SPHERICAL TRIANGLE', OR THE NATIONAL BUREAU OF STANDARDS
C 'HANDBOOK OF MATHEMATICAL FUNCTIONS' UNDER SECTION 4.3.149,
C FORMULAS FOR SOLUTION OF SPHERICAL TRIANGLES.
C DOUBLE PRECISION IS USED INTERNALLY BECAUSE OF THE WIDE
C RANGE OF GEOGRAPHIC VALUES THAT MAY BE USED.
C
C ATTRIBUTES:
C LANGUAGE: CRAY CFT77 FORTRAN
C MACHINE: CRAY Y-MP8/832
C
C$$$
C
C *** IMPLICIT TYPE DEFAULTS.....
C
IMPLICIT REAL (A-H,O-Z)
C
C *** CONSTANTS......
C
REAL PI
REAL HALFPI
REAL DR
REAL RD
REAL TDEG, TRAD, TPOD, TFLT
REAL EARTHR
REAL WHOLCD, HALFCD, QUARCD
C
C *** VARIABLES......
C
REAL RLAT1, RLAT2, COSLA1, COSLA2, SINLA1, SINLA2
REAL DLOND, RLOND, COSLO, SINLO, SANGG, ABEAR
REAL YDISP, XDISP, DDLAT1, DDLAT2, DBANG
REAL DLAT1, DLAT2, DLON1, DLON2, BEARD, GCDKM
C
C *** CONVERT LATITUDES AND LONGITUDE DIFFERENCE TO RADIANS.
C
DATA PI /3.141592653589793238462643/
DATA HALFPI/1.570796326794896619231322/
DATA DR /0.017453292519943295769237/
DATA RD /57.295779513082320876798155/
DATA TDEG /1E-10/, TRAD/1E-10/, TPOD/1E-6/, TFLT/1E-6/
DATA EARTHR/6371.2/
DATA WHOLCD/360.0/, HALFCD/180.0/, QUARCD/90.0/
DDLAT1 = DLAT1
DDLAT2 = DLAT2
RLAT1 = DR * DDLAT1
RLAT2 = DR * DDLAT2
DLOND = DLON2 - DLON1
IF (DLOND .GT. HALFCD) DLOND = DLOND - WHOLCD
IF (DLOND .LT. -HALFCD) DLOND = DLOND + WHOLCD
RLOND = DR * DLOND
C
C *** FIRST WE ATTACK THE CASES WHERE POINT 1 IS VERY CLOSE TO THE
C *** NORTH OR SOUTH POLES.
C *** HERE WE USE CONVENTIONAL VALUE FOR BEARING.. - LONG (2) AT THE
C *** NORTH POLE, AND + LONG (2) AT THE SOUTH POLE.
C
IF (ABS(DDLAT1-QUARCD) .LT. TDEG) THEN
IF (ABS(DDLAT2-QUARCD) .LT. TDEG) THEN
DBANG = 0.0
SANGG = 0.0
ELSE IF (ABS(DDLAT2+QUARCD) .LT. TDEG) THEN
DBANG = 0.0
SANGG = PI
ELSE
DBANG = -DLON2
SANGG = HALFPI - RLAT2
ENDIF
ELSE IF (ABS(DDLAT1+QUARCD) .LT. TDEG) THEN
IF (ABS(DDLAT2-QUARCD) .LT. TDEG) THEN
DBANG = 0.0
SANGG = PI
ELSE IF (ABS(DDLAT2+QUARCD) .LT. TDEG) THEN
DBANG = 0.0
SANGG = 0.0
ELSE
DBANG = +DLON2
SANGG = HALFPI + RLAT2
ENDIF
C
C *** NEXT WE ATTACK THE CASES WHERE POINT 2 IS VERY CLOSE TO THE
C *** NORTH OR SOUTH POLES.
C *** HERE BEARING IS SIMPLY 0 OR 180 DEGREES.
C
ELSE IF (ABS(DDLAT2-QUARCD) .LT. TDEG) THEN
DBANG = 0.0
SANGG = HALFPI - RLAT1
ELSE IF (ABS(DDLAT2+QUARCD) .LT. TDEG) THEN
DBANG = HALFCD
SANGG = HALFPI + RLAT1
C
C *** THE CASE REMAINS THAT NEITHER POINT IS AT EITHER POLE.
C *** FIND COSINE AND SINE OF LATITUDES AND LONGITUDE DIFFERENCE
C *** SINCE THEY ARE USED IN MORE THAN ONE FORMULA.
C
ELSE
COSLA1 = COS(RLAT1)
SINLA1 = SIN(RLAT1)
COSLA2 = COS(RLAT2)
SINLA2 = SIN(RLAT2)
COSLO = COS(RLOND)
SINLO = SIN(RLOND)
C
C *** FOLLOWING IS FORMULA FOR GREAT CIRCLE SUBTENDED ANGLE BETWEEN
C *** POINTS IN RADIAN MEASURE.
C
SANGG = ACOS(SINLA1*SINLA2 + COSLA1*COSLA2*COSLO)
C
C *** IF THE GREAT CIRCLE SUBTENDED ANGLE IS VERY SMALL, FORCE BOTH
C *** BEARING AND DISTANCE TO BE ZERO.
C
IF (ABS(SANGG) .LT. TRAD) THEN
DBANG = 0.0
SANGG = 0.0
C
C *** IF THE GREAT CIRCLE SUBTENDED ANGLE IS JUST SMALL, ASSUME A
C *** FLAT EARTH AND CALCULATE Y- AND X-DISPLACEMENTS. THEN FIND
C *** BEARING USING THE ARCTANGENT FUNCTION AND DISTANCE USING THE
C *** SQUARE ROOT.
C
ELSE IF (ABS(SANGG) .LT. TFLT) THEN
YDISP = RLAT2-RLAT1
XDISP = RLOND*COSLA2
ABEAR = ATAN2(XDISP, YDISP)
DBANG = RD*ABEAR
SANGG = SQRT(YDISP**2 + XDISP**2)
C
C *** IF THE ANGLE IS RATHER CLOSE TO PI RADIANS, FORCE BEARING TO
C *** BE ZERO AND DISTANCE TO BE PI.
C *** THE TEST FOR 'CLOSE TO PI' IS MORE RELAXED THAN THE TEST FOR
C *** 'CLOSE TO ZERO' TO ALLOW FOR GREATER RELATIVE ERROR.
C
ELSE IF (ABS(SANGG-PI) .LT. TPOD) THEN
DBANG = 0.0
SANGG = PI
C
C *** OTHERWISE COMPUTE THE PRINCIPAL VALUE OF THE BEARING ANGLE
C *** USING THE LAW OF SINES. THE DIVISION BY THE SINE FORCES US TO
C *** LIMIT THE DOMAIN OF THE ARCSINE TO (-1,1).
C
ELSE
ABEAR = ASIN(AMAX1(-1.0,AMIN1(+1.0,COSLA2*SINLO/
& SIN(SANGG))))
C
C *** IF THE LONGITUDE DIFFERENCE IS LESS THAN PI/2 IT IS NECESSARY
C *** TO CHECK WHETHER THE BEARING ANGLE IS ACUTE OR OBTUSE BY
C *** COMPARING LATITUDE (1) WITH THE LATITUDE OF THE GREAT CIRCLE
C *** THROUGH POINT (2) NORMAL TO MERIDIAN OF LONGITUDE (1). IF
C *** LATITUDE (1) IS GREATER, BEARING IS OBTUSE AND THE ACTUAL
C *** BEARING ANGLE IS THE SUPPLEMENT OF THE ANGLE CALCULATED ABOVE.
C
IF (0.0 .LE. COSLA1*SINLA2 .AND. COSLA1*SINLA2 .LE.
& COSLA2*SINLA1*COSLO .OR. COSLA1*SINLA2 .LE. 0.0 .AND.
& COSLA2*SINLA1*COSLO .GE. COSLA1*SINLA2) ABEAR =
& SIGN(PI,ABEAR) - ABEAR
DBANG = RD * ABEAR
ENDIF
ENDIF
C
C *** THIS FINISHES THE CASE WHERE POINTS ARE NOT AT THE POLES.
C *** NOW CONVERT BEARING TO DEGREES IN RANGE -180 TO +180 AND FIND
C *** GREAT CIRCLE DISTANCE IN KILOMETERS.
C
IF (DBANG .LE. -HALFCD) DBANG = DBANG + WHOLCD
IF (DBANG .GT. HALFCD) DBANG = DBANG - WHOLCD
GCDKM = EARTHR * SANGG
BEARD = DBANG
RETURN
END