C*********************************************************************** C*********************************************************************** C* * C* GEODETIC DIRECT AND INVERSE D.WELLS JAN 1972 * C* BASED ON CLARKE'S BEST FORMULA (NORMAL SECTION APPROACH) * C* SEE ROBBINS 'LONG LINES ON THE SPHEROID' E.S.R. #125,P.301,1962 * C* CLAIMED TO BE GOOD TO 1 CM AT 1600 KM * C* * C*********************************************************************** IMPLICIT REAL*8(A-H,O-Z) REAL*8 NAME(3) PI=3.141592653589793D0 C C ...READ AND PRINT ELLIPSOID NAME, A, F C READ 100,NAME,A,F 100 FORMAT(10X,3A8,5X,F15.5,5X,F10.5) PRINT 110,NAME,A,F 110 FORMAT('1',10X,3A8,5X,F15.5,5X,F15.5,//) IF(F.GT.1000.) F=A/(A-F) F=1D0/F E2=F*(2D0-F) EP=E2/(1D0-E2) C C ...READ AND PRINT LAT1,LONG1,LAT2,LONG2,AZ12,AZ21,DIST C PRINT 115 115 FORMAT(7X,' LATITUDE-1',4X,'LONGITUDE-1',4X,' LATITUDE-2',4X, 1 'LONGITUDE-2',4X,'AZIMUTH 1-2',4X,'AZIMUTH 2-1 DISTANCE(M)') 120 CONTINUE READ 150,IP1,MP1,SP1,IL1,ML1,SL1,IP2,MP2,SP2,IL2,ML2,SL2, 1 IA12,MA12,SA12,IA21,MA21,SA21,DIST 150 FORMAT(2I2,F2.0,I4,I2,F2.0,2I3,F7.3,I4,I3,F7.3,2I3,F6.2,I4,I3, 1 F6.2,F11.2) IF(IP1.EQ.99) GOTO 300 PRINT 151,IP1,MP1,SP1,IL1,ML1,SL1,IP2,MP2,SP2,IL2,ML2,SL2 151 FORMAT(' IN ' ,2(2I3,F7.3,2X),2(2I3,F5.1,4X)) C C ...CONVERT TO RADIANS C PHI1=(IP1 +MP1 /60D0+SP1 /3600D0)*PI/180D0 XLNG1=(IL1 +ML1 /60D0+SL1 /3600D0)*PI/180D0 APHI2=(IP2 +MP2 /60D0+SP2 /3600D0)*PI/180D0 XLNG2=(IL2 +ML2 /60D0+SL2 /3600D0)*PI/180D0 ALF12=(IA12+MA12/60D0+SA12/3600D0)*PI/180D0 ALF21=(IA21+MA21/60D0+SA21/3600D0)*PI/180D0 C C ...THIS PROGRAM ASSUMES EAST LONGITUDES POSITIVE C AND AZIMUTHS CLOCKWISE FROM NORTH C XLNG1=-XLNG1 XLNG2=-XLNG2 C C ...CALL INVERSE SUBROUTINE, CONVERT RESULTS TO D/M/S AND PRINT C CALL INVERS(A,E2,EP,PHI1,XLNG1,APHI2,XLNG2 ,AZ12,AZ21,ADIST) C IAZ1=IDEG(AZ12,MAZ1,SAZ1) IAZ2=IDEG(AZ21,MAZ2,SAZ2) C PRINT 210,IAZ1,MAZ1,SAZ1,IAZ2,MAZ2,SAZ2,ADIST 210 FORMAT(' INV ' ,60X,2(2I3,F7.3,2X),F11.2,//) C ...CALL DIRECT SUBROUTINE, CONVERT RESULTS TO D/M/S AND PRINT C CALL DIRECT(A,E2,EP,PHI1,XLNG1,AZ12,ADIST,APHI2,XLONG2,AZ21) C IPH2=IDEG(APHI2,MPH2,SPH2) ILG2=IDEG(-XLONG2,MLG2,SLG2) IAZ2=IDEG(AZ21,MAZ2,SAZ2) C PRINT 200,IPH2,MPH2,SPH2,ILG2,MLG2,SLG2,IAZ2,MAZ2,SAZ2 200 FORMAT(' DIR ' ,30X,2(2I3,F7.3,2X),15X,2I3,F7.3,///) C C C ...RETURN FOR NEXT POINT C GOTO 120 300 PRINT 310 310 FORMAT(10X,'END OF DATA') RETURN END FUNCTION IDEG(RAD,MIN,SEC) C C ...CONVERT RADIANS TO DEGREES,MINUTES,SECONDS C IMPLICIT REAL*8(A-H,O-Z) PI=3.141592653589793D0 DEG=RAD*180D0/PI+1D-7 IDEG=DEG MIN=(DEG-IDEG)*60D0 SEC=((DEG-IDEG)*60D0-MIN)*60D0 RETURN END SUBROUTINE DIRECT(A,E2,EP,PHI1,XLNG1,ALF12,DIST,PHI2,XLNG2,ALF21) C*********************************************************************** C* * C* GEODETIC DIRECT PROBLEM D.WELLS JAN 1972 * C* GIVEN PHI1,XLNG1,ALF12,DIST COMPUTE PHI2,XLNG2,ALF21 * C* * C* SPHEROID PARAMETERS - EQUATORIAL RADIUS 'A' * C* - SQUARE OF FIRST ECCENTRICITY 'E2' * C* - SQUARE OF SECOND ECCENTRICITY 'EP' * C* ALL ANGLES IN RADIANS, ALL DISTANCES SAME UNITS AS 'A' * C* * C* USES ROBBINS PRESENTATION OF CLARKE'S BEST FORMULA, WHICH IS * C* BASED ON COMPUTATIONS ON THE NORMAL SECTION SPHERE. * C* SEE ROBBINS 'LONG LINES ON THE SPHEROID', EMPIRE SURVEY REVIEW * C* NO.125 PAGES 301-309, APRIL 1962 * C* * C*********************************************************************** IMPLICIT REAL*8(A-H,O-Z) PI=3.141592653589793D0 XN1=A/DSQRT(1D0-E2*DSIN(PHI1)**2) G2=EP*DSIN(PHI1)**2 H2=EP*(DCOS(PHI1)*DCOS(ALF12))**2 GH=EP*DSIN(PHI1)*DCOS(PHI1)*DCOS(ALF12) THETA=DIST/XN1 SIGMA=THETA*(1D0+THETA**2*H2*(1D0-H2)/6D0-THETA**3*GH* 1 (1D0-2D0*H2)/8D0-THETA**4*(H2*(4D0-17D0*H2)-3D0*G2* 2 (1D0-7D0*H2))/120D0+THETA**5*GH/48D0) SPSI2=DSIN(PHI1)*DCOS(SIGMA)+DCOS(PHI1)*DCOS(ALF12)*DSIN(SIGMA) PSI2=DARSIN(SPSI2) XLNG2=XLNG1+DARSIN(DSIN(SIGMA)*DSIN(ALF12)/DCOS(PSI2)) XN1OR2=1D0+EP*(SPSI2-DSIN(PHI1))**2/2D0 PHI2=DATAN(DTAN(PSI2)*(1D0+EP)*(1D0-E2*XN1OR2*DSIN(PHI1)/SPSI2)) AZ21P=DATAN2(DSIN(ALF12),DCOS(SIGMA)*DCOS(ALF12)-DSIN(SIGMA)* 1 DTAN(PHI1)) ALF21=AZ21P-(PHI2-PSI2)*DSIN(AZ21P)*DTAN(SIGMA/2D0) +PI RETURN END SUBROUTINE INVERS(A,E2,EP,PHI1,XLNG1,PHI2,XLNG2,ALF12,ALF21,DIST) C*********************************************************************** C* * C* GEODETIC INVERSE PROBLEM D.WELLS JAN 1972 * C* GIVEN PHI1,XLNG1,PHI2,XLNG2 COMPUTE ALF12,ALF21,DIST * C* * C* SPHEROID PARAMETERS - EQUATORIAL RADIUS 'A' * C* - SQUARE OF FIRST ECCENTRICITY 'E2' * C* - SQUARE OF SECOND ECCENTRICITY 'EP' * C* ALL ANGLES IN RADIANS, ALL DISTANCES SAME UNITS AS 'A' * C* * C* USES ROBBINS PRESENTATION OF CLARKE'S BEST FORMULA, WHICH IS * C* BASED ON COMPUTATIONS ON THE NORMAL SECTION SPHERE. * C* SEE ROBBINS 'LONG LINES ON THE SPHEROID', EMPIRE SURVEY REVIEW * C* NO.125 PAGES 301-309, APRIL 1962 * C* * C*********************************************************************** IMPLICIT REAL*8(A-H,O-Z) PI=3.141592653589793D0 XN1=A/DSQRT(1D0-E2*DSIN(PHI1)**2) XN2=A/DSQRT(1D0-E2*DSIN(PHI2)**2) TPSI1=DTAN(PHI1)*(1D0-E2+E2*XN2*DSIN(PHI2)/(XN1*DSIN(PHI1))) TPSI2=DTAN(PHI2)*(1D0-E2+E2*XN1*DSIN(PHI1)/(XN2*DSIN(PHI2))) DLONG=XLNG2-XLNG1 ALF12=DATAN2(DSIN(DLONG),DCOS(PHI1)*TPSI2-DSIN(PHI1)*DCOS(DLONG)) ALF21=DATAN2(DSIN(-DLONG),DCOS(PHI2)*TPSI1-DSIN(PHI2)*DCOS(-DLONG 1 )) IF(ALF12.LT.0.D0)ALF12=ALF12+2.*PI IF(ALF21.LT.0.) ALF21=2D0*PI+ALF21 CHAP=DCOS(PHI1)*TPSI2-DSIN(PHI1)*DCOS(DLONG) SIGMA=DARSIN(DSQRT((DSIN(DLONG)**2+CHAP**2)/(1D0+TPSI2**2))) G2=EP*DSIN(PHI1)**2 H2=EP*(DCOS(PHI1)*DCOS(ALF12))**2 GH=EP*DSIN(PHI1)*DCOS(PHI1)*DCOS(ALF12) DIST=XN1*SIGMA*(1D0-SIGMA**2*H2*(1D0-H2)/6D0+SIGMA**3*GH*(1D0-2D0 1 *H2)/8D0+SIGMA**4*(H2*(4D0-7D0*H2)-3D0*G2*(1D0-7D0*H2))/120D0 2 -SIGMA**5*GH/48D0) RETURN END