SUBROUTINE VININ (AE,F,XLAT1,XLONG1,XLAT2,XLONG2,DIST,AZ1) IMPLICIT REAL * 8 (A - H, O - Z) C STATEMENT FUNCTIONS: SIN(A) = DSIN(A) COS(A) = DCOS(A) ACOS(A)=DARCOS(A) TAN(A) = DSIN(A)/DCOS(A) ATAN(A) = DATAN(A) ATAN2(A,B)=DATAN2(A,B) SQRT(A) = DSQRT(A) ABS(A) = DABS(A) PI=ACOS(-1.D0) C$ C NAME VININ C C TYPE SUBROUTINE C C PURPOSE INVERSE GEODETIC PROBLEM C C AUTHOR/DATE P.DELORME 780830 USING ALGORITHM OF T.VINCENTY C C EXTERNALS DSIN,DCOS,DARCOS,DATAN,DATN2,DSQRT,DABS C C CALLING CALL VININ(AE,F,XLAT1,XLONG1,XLAT2,XLONG2,DIST,AZ1) C C PARAMETERS AE = ELLIPSOID SEMI-MAJOR AXIS (M) C F = INVERSE OF ELLIPSOID FLATTENING C XLAT1 = INPUT LAT OF POINT 1 (DEG) C XLONG1 = INPUT LONG OF POINT 1 (DEG) C XLAT2 = INPUT LAT OF POINT 2 (DEG) C XLONG2 = INPUT LONG OF POINT 2 (DEG) C DIST = OUTPUT GEODESIC DISTANCE (M) C AZ1 = OUTPUT FORWARD AZIMUTH (DEG) C AZ2 = OUTPUT BACK AZIMUTH (DEG) (IF REQUIRED ADD AZ2 C TO END OF ARGUMENT LIST AND ACTIVATE APPENDIX) C C LANGUAGE FORTRAN C$ C FL = 1D0 / F FUZ = 20.0D0 * DEPS(DARG) BE = AE *(1.0 - FL) C C CONVERT LATITUDE AND LONGITUDE OF POINTS 1 AND 2 TO RADIANS. C ALAT1 = (XLAT1 * PI) / 180D0 ALON1 = (XLONG1 * PI) / 180D0 ALAT2 = (XLAT2 * PI) / 180D0 ALON2 = (XLONG2 * PI) / 180D0 C C CALCULATE THE REDUCED LATITUDE AND ITS TRIG FUNCTIONS: C TU1 = (1D0 - FL) * TAN(ALAT1) TU2 = (1D0 - FL) * TAN(ALAT2) U1 = ATAN(TU1) U2 = ATAN(TU2) SU1 = SIN(U1) SU2 = SIN(U2) CU1 = COS(U1) CU2 = COS(U2) C C 1ST APPROX: DIFF IN LONG ON ELLIPSOID = DIFF IN LONG ON SPHERE: C DL = ALON2 - ALON1 XDL = DL C C C ********************* BEGINNING OF ITERATION ***************** C C 100 CONTINUE SDL = SIN(DL) CDL = COS(DL) CS = (SU1 * SU2) + (CU1 * CU2 * CDL) SS = SQRT(1D0 - (CS**2)) SIG = ATAN2(SS , CS) IF(ABS(SS) .LT. FUZ) SS = FUZ SA = (CU1 * CU2 * SDL)/SS CA = SQRT(1D0 - (SA**2)) C2SM = CS - (2D0 * SU1 * SU2) / (CA**2) C C CALCULATE THE DIFFERENCE IN LONGITUDE ON THE AUXILIARY SPHERE: C C = (FL / 16D0) * (CA**2) * (4D0 + (FL * (4D0 - 3D0 * (CA**2)))) DL1 = XDL + (1D0 - C) * FL * SA * (SIG + C * SS * (C2SM + C $ * CS * (-1D0 + 2D0 * (C2SM**2)))) IF(ABS(DL1 - DL) .LE. 10D-10) GOTO 45 DL = DL1 GOTO 100 45 CONTINUE C C C ********************* END OF ITERATION ********************* C C C C CALCULATE U SQUARE: C U = (CA**2) * ((AE**2) - (BE**2)) / (BE**2) C C CALCULATE A: C A = 1D0 + (U / 256D0) * (64D0 + U * (-12D0 + 5D0 *U)) C C CALCULATE B: C B = (U / 512D0) * (128D0 + U * (-64D0 + (37D0 * U))) C C CALCULATE DELTA SIGMA: C DSIG = B * SS * (C2SM + 0.25D0 * B * CS * (-1D0 +2D0*(C2SM**2))) C C CALCULATE THE GEODESIC: C DIST = BE * A * (SIG - DSIG) C C CALCULATE THE FORWARD AZIMUTH: C SDL1 = SIN(DL1) CDL1 = COS(DL1) AZ1 = ATAN2((CU2 * SDL1),(CU1 * SU2 - SU1 * CU2 * CDL1)) $ * 180D0 / PI IF (AZ1 .GT. 0D0) GOTO 70 AZ1 = AZ1 + 360D0 70 CONTINUE C C C APPENDIX ************** CALCULATE THE BACKWARD AZIMUTH: ***** C C AZ2 = ATAN2((-1D0 * CU1 * SDL1),(SU1 * CU2 - CU1 * SU2 * CDL1 C $ * 180D0 / PI C IF(AZ2 .GT. 0D0) GOTO 75 C AZ2 = AZ2 + 360D0 C75 CONTINUE C C C ******************************************************************* C C C RETURN END