SUBROUTINE VINDI(AE,F,XLAT1,XLONG1,FAZ,LINE, $ LAT2,LON2,BA) DOUBLE PRECISION A,ABS,AE,ATAN,B,BA,BE,C,CA,C2SM,CS,COS,CU1,CFAZ, $ DSIN,DCOS,DATAN,DSQRT,DABS,DSIG,DL,DLON,FAZ, $ LAT2,LINE,LON2,NEWSIG,TAN,TU1,TWOSM,XLAT1, $ XLONG1,U1,USQ,PI,F,FL,SIN,SQRT,SU1,SIG1, $ SFAZ,SA,SIG,SS,ATAN2,DATAN2 DATA PI/3.1415926535D0/ C STATEMENT FUNCTIONS: SIN(A) = DSIN(A) COS(A) = DCOS(A) TAN(A) = DSIN(A) / DCOS(A) SQRT(A) = DSQRT(A) ABS(A) = DABS(A) ATAN(A) = DATAN(A) ATAN2(A,B) = DATAN2(A,B) C$ C NAME VINDI C C TYPE SUBROUTINE C C PURPOSE DIRECT GEODETIC PROBLEM C C AUTHOR/DATE P.DELORME 780830 USING ALGORITHM OF T. VINCENTY C C EXTERNALS DSIN,DCOS,DSQRT,DABS,DATAN,DATN2 C C CALLING CALL VINDI(AE,F,XLAT1,XLONG1,FAZ,LINE,LAT2,LON2,BA) C C PARAMETERS AE = ELLIPSOID SEMI-MAJOR AXIS (M) C F = INVERSE OF ELLIPSOID FLATTENING C XLAT1 = INPUT LAT OF WESTERN POINT (DEG) C XLONG1 = INPUT LONG OF WESTERN POINT (DEG) C FAZ = FORWARD AZIMUTH(DEG) C LINE = GEODESIC DISTANCE (M) C LAT2 = OUTPUT LAT OF EASTERN POINT (DEG) C LON2 = OUTPUT LONG OF EASTERN POINT (DEG) C BA = BACKWARD AZIMUTH (DEG) C C LANGUAGE FORTRAN C$ C C MISCELLANEOUS VALUES: C FL = 1D0 / F BE = AE * (1D0 - FL) C C CALCULATE THE REDUCED LATITUDE AND ITS TRIG FUNCTIONS: C TU1 = (1D0 - FL) * TAN(XLAT1 * PI / 180D0) U1 = ATAN(TU1) SU1 = SIN(U1) CU1 = COS(U1) C C CALCULATE THE ANGULAR DISTANCE ON THE SPHERE FROM THE EQUATOR TO P1: C CFAZ = COS(FAZ * PI / 180D0) SIG1 = ATAN(TU1 / CFAZ) C C CALCULATE THE AZIMUTH OF THE GEODESIC AT THE EQUATOR: C SFAZ = SIN(FAZ * PI / 180D0) SA = CU1 * SFAZ CA = SQRT(1D0 - (SA**2)) C C CALCULATE U SQUARE: C USQ = (CA**2) * ((AE**2 - BE**2) / (BE**2)) C C CALCULATE A: C A = 1D0 + ((USQ / 256D0) * (64D0 + USQ * (-12D0 + 5D0 * USQ))) C C CALCULATE B: C B = (USQ / 512D0) * (128D0 + USQ * (-64D0 + 37D0 * USQ)) C C SIGMA'S FIRST APPROXIMATION: C SIG = LINE / (BE * A ) C C C **************** BEGINNING OF ITERATION ******************* C C 100 CONTINUE TWOSM = (2D0 * SIG1) + SIG C2SM = COS(TWOSM) CS = COS(SIG) SS = SIN(SIG) C C CALCULATE DELTA SIGMA: C DSIG = B * SS * (C2SM + 0.25D0 * B * CS * (-1D0 + 2D0 $ * (C2SM**2))) NEWSIG = (LINE / (BE * A)) + DSIG IF(ABS(NEWSIG - SIG) .LE. 10D-10) GOTO 70 SIG = NEWSIG GOTO 100 70 CONTINUE C C C ******************* END OF ITERATION *************************** C C C CALCULATE THE LATITUDE OF EASTERN POINT: C CS = COS(NEWSIG) SS = SIN(NEWSIG) LAT2 = ATAN2(((SU1 * CS) + (CU1 * SS * CFAZ)),((1D0 - FL) $ * SQRT((SA**2) + ((SU1 * SS) - (CU1 * CS * CFAZ)) $ **2))) * 180D0 / PI C C CALCULATE THE DIFF IN LONGITUDE ON THE AUXILIARY SPHERE: C DL = ATAN2((SS * SFAZ),((CU1*CS) - (SU1 * SS * CFAZ))) C C CALCULATE C: C C = (FL / 16D0) * (CA**2) * (4D0 + FL * (4D0 - 3D0 * (CA**2))) C C CALCULATE THE DIFF IN LONGITUDE ON THE ELLIPSOID: C DLON = (DL - (1D0 - C) * FL * SA * (NEWSIG + C * SS * (C2SM $ + C * CS * (-1D0 + 2D0 * (C2SM**2))))) * 180D0 / PI C C CALCULATE THE LONGITUDE OF EASTERN POINT: C LON2 = XLONG1 + DLON C C C CALCULATE THE BACKWARD AZIMUTH: C BA = ATAN2((-1D0 * SA),((SU1 * SS) - (CU1 * CS * CFAZ))) $ * 180D0 / PI IF(BA .GE. 0D0) GOTO 80 BA = BA + 360D0 80 CONTINUE C C C C RETURN END