C*********************************************************************** C* * C* PROGRAM TO COMPUTE THE DIRECT AND INVERSE PROBLEMS OF GEODESY * C* USING PUISSANTS FORMULAS (SEE SPECIFICATIONS IN SUBROUTINES DIRECT * C* AND INVRSE FOR SOURCE AND ACCURACIES) * C* * C* C. CHAMBERLAIN NOV. 1974 * C* * C* PROGRAM INPUT AND FORMAT SPECIFICATIONS * C* * C* CARD#1: JOB DESCRIPTION (20A4) * C* CARD#2: REFERENCE ELLIPSOID NAME, SEMI-MAJOR AND SEMI-MINOR AXES * C* OR SEMI-MAJOR AXIS AND RECIPROCAL OF THE FLATTENING * C* (5A4,2F15.0) * C* CARD#3: STATION NUMBERS OF THE TWO POINTS (2I9) * C* CARD#4: LATITUDE POINT 1 (2I2,F7.0 * C* LONGTITUDE POINT 1 I3,I2,F7.0 * C* LATITUDE POINT 2 2I2,F7.0 * C* LONGTITUDE POINT 2 I3,I2,F7.0 * C* AZIMUTH 1 - 2 I3,I2,F7.0 * C* DISTANCE 1 - 2 T60,F10.0) * C* NOTE 1: CARDS 3 AND 4 ARE REPEATED FOR EACH LINE TO BE COMPUTED * C* NOTE 2: ALL ANGLES ARE ENTERED IN DEGREES MINUTES AND SECONDS * C* FOR NEGATIVE ANGLES ONLY PUT MINUS SIGN ON THE DEGREES * C* NOTE 3: FOR DIRECT PROBLEM LEAVE LATITUDE 2 AND LONGTITUDE 2 * C* BLANK * C* NOTE 4: FOR INVERSE PROBLEM LEAVE AZ 1 - 2 AND DISTANCE BLANK * C* NOTE 5: ALL INPUT IS TERMINATED BY A BLANK CARD FOR A NORMAL END * C* * C* * C*********************************************************************** IMPLICIT REAL*8(A-H,O-$) REAL*4 S1,S2,S3,S4,S5 DIMENSION NAME(20) C READ IN AND WRITE OUT JOB DESCRIPTION C READ(5,1)(NAME(I),I=1,20) 1 FORMAT(20A4) WRITE(6,2)(NAME(I),I=1,20) 2 FORMAT('1'//20X,20A4//) C C C READ IN AND WRITE OUT ELLIPSOID INFORMATION C READ(5,3)(NAME(I),I=1,5),A,B 3 FORMAT(5A4,2F15.0) IF(B.LT.1000.)WRITE(6,5)(NAME(I),I=1,5),A,B IF(B.GT.1000.)WRITE(6,4)(NAME(I),I=1,5),A,B 4 FORMAT(5X,5A4,5X,'A=',F12.2,5X,'B=',F12.2) 5 FORMAT(5X,5A4,5X,'A=',F12.2,5X,'1/F=',F12.6) WRITE(6,6) 6 FORMAT('-',5X,'DIRECT AND INVERSE PROBLEMS USING PUISSANTS FORMULA .S'///) C C READ IN A NEW SET OF POINTS FOR COMPUTATION C 8 READ(5,12)NUM1,NUM2 12 FORMAT(2I9) C C CHECK IF CARD IS BLANK C IF(NUM1.EQ.0.AND.NUM2.EQ.0)GO TO 50 READ(5,7)I1,I2,S1,J1,J2,S2,K1,K2,S3,L1,L2,S4,M1,M2,S5,SIJ 7 FORMAT(2(2I2,F7.0,I3,I2,F7.0),I3,I2,F7.0,T60,F10.0) C C CHANGE LATITUDE AND LONGTITUDE OF POINT 1 TO RADIANS C CALL ARCRAD(I1,I2,S1,PHII) CALL ARCRAD(J1,J2,S2,XLONI) C C CHECK IF DISTANCE = 0 IF TRUE DO INVERSE PROBLEM C IF(SIJ.EQ.0.0)GO TO 20 C CHANGE AZIMUTH 1 - 2 TO RADIANS AND DO DIRECT PROBLEM C CALL ARCRAD(M1,M2,S5,AIJ) CALL DIRECT(A,B,PHII,XLONI,AIJ,SIJ,PHIJ,XLONJ,AJI) C C CHANGE LATITUDE AND LONGTITUDE OF POINT 2 PLUS AZIMUTH 1 - 2 TO C DEGREES MINUTES AND SECONDS AND WRITE OUT THE SOLUTION C CALL RADARC(PHIJ,I1,I2,S1) CALL RADARC(XLONJ,J1,J2,S2) CALL RADARC(AJI,K1,K2,S3) WRITE(6,13)NUM1,NUM2 13 FORMAT('0',35X,'STATION ',I9,' TO STATION ',I9 ./'+',35X,'_______ __ _______') WRITE(6,9)I1,I2,S1,J1,J2,S2,K1,K2,S3 9 FORMAT('0','DIRECT',5X,'LATITUDE 2=',I3,2X,I2,2X,F7.3,5X,'LONGTITU .DE 2=',I4,2X,I2,2X,F7.3,5X,'BACK AZ=',I4,2X,I2,2X,F7.3) C C CHECK SOLUTION BY DOING INVERSE SOLUTION AND WRITE OUT THE RESULTS C CALL INVRSE(A,B,PHII,XLONI,PHIJ,XLONJ,AIJ,AJI,SIJ) CALL RADARC(AIJ,I1,I2,S1) WRITE(6,10)I1,I2,S1,SIJ 10 FORMAT('0',5X,'CHECK BY INVERSE SOLUTION',5X,'AZIMUTH 1 - 2=',I4,2 .X,I2,2X,F7.3,5X,'DISTANCE =',F14.3,'M'///) C C GO BACK AND READ IN A NEW SET OF POINTS C GO TO 8 C C CHANGE LATITUDE AND LONGTITUDE OF 2 TO RADIANS AND DO THE INVERSE C SOLUTION C 20 CALL ARCRAD(K1,K2,S3,PHIJ) CALL ARCRAD(L1,L2,S4,XLONJ) CALL INVRSE(A,B,PHII,XLONI,PHIJ,XLONJ,AIJ,AJI,SIJ) C C CHANGE RESULTS TO DEGREES MINUTES AND SECONDS AND WRITE OUT C CALL RADARC(AIJ,I1,I2,S1) CALL RADARC(AJI,J1,J2,S2) WRITE(6,13)NUM1,NUM2 WRITE(6,21)I1,I2,S1,J1,J2,S2,SIJ 21 FORMAT('0','INVERSE',5X,'AZIMUTH 1-2 =',I4,2X,I2,2X,F7.3, .5X,'BACK AZ=',I4,2X,I2,2X,F7.3,5X,'DISTANCE =',F14.3,'M') C C CHECK SOLUTION BY DIRECT PROBLEM C CALL DIRECT(A,B,PHII,XLONI,AIJ,SIJ,PHIJ,XLONJ,AJI) C C CHANGE RESULTS TO DEGREES MINUTES AND SECONDS AND WRITE OUT C CALL RADARC(PHIJ,I1,I2,S1) CALL RADARC(XLONJ,J1,J2,S2) WRITE(6,22)I1,I2,S1,J1,J2,S2 22 FORMAT('0',5X,'CHECK BY DIRECT SOLUTION',5X,'LATITUDE 2=',I4,2X,I2 .,2X,F7.3,5X,'LONGTITUDE 2=',I4,2X,I2,2X,F7.3///) C C GO READ IN ANOTHER SET OF POINTS C GO TO 8 50 STOP END SUBROUTINE DIRECT(A,B,PHII,XLONI,AIJ,SIJ,PHIJ,XLONJ,AJI) C*********************************************************************** C* * C* SUBROUTINE DIRECT SOLVES THE DIRECT PROBLEM IN GEODESY, THAT * C* IS, GIVEN THE LATITUDE AND LONGTITUDE OF A POINT I, PLUS THE * C* AZIMUTH AND DISTANCE TO A POINT J, COMPUTE THE LATITUDE AND * C* LONGTITUDE OF J, PLUS THE BACKWARD AZIMUTH FROM J TO I. * C* * C* THE EQUATIONS PROGRAMMED ARE 13, 15, 90, 99, 105, 101, 115 * C* AND 109(A) IN 'GEODETIC POSITION COMPUTATIONS' EXCEPT THAT IN ALL * C* CASES WHERE IT APPLIES, THE COMPUTATIONS ARE IN RADIANS RATHER * C* THAN IN SECONDS OF ARC. THESE EQUATIONS ARE KNOWN AS 'PUISSANT'S * C* FORMULA FOR SHORT LINES' WHICH GIVE A REPUTED ACCURACY OF 1 PPM AT * C* 100KM. (IE. THEY ARE BASED ON A SPHERICAL APPROXIMATION TO THE * C* ELLIPSOID.) * C* * C* THE SOLUTION IS ITERATED UNTIL THE CHANGE IN DPHI IS LESS THAN * C* 1*10-10 RADIANS. * C* * C* NOTE: THE EQUIVALENT INVERSE PROGRAM IS SUBROUTINE 'INVRSE' * C* * C* INPUT DATA: * C* A=SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID IN METRES * C* B=SEMI-MINOR AXIS OF REFERENCE ELLIPSOID OR THE * C* RECIPROCAL OF THE FLATTENING IN METRES OR UNITLESS * C* PHII=GEODETIC LATITUDE OF POINT I IN RADIANS * C* XLONI=GEODETIC LONGTITUDE OF POINT I IN RADIANS * C* AIJ=GEODETIC AZIMUTH FROM I TO J IN RADIANS * C* NOTE: AZIMUTH ON ELLIPSOID * C* SIJ=DISTANCE FROM I TO J IN METRES * C* NOTE: DISTANCE ON ELLIPSOID * C* * C* OUTPUT DATA: * C* PHIJ=GEODETIC LATITUDE OF POINT J IN RADIANS * C* XLONJ=GEODETIC LONGTITUDE OF POINT J IN RADIANS * C* AJI=GEODETIC AZIMUTH FROM J TO I IN RADIANS * C* NOTE: AZIMUTH ON ELLIPSOID * C* * C* ( C. CHAMBERLAIN NOV. 1974 ) * C* * C*********************************************************************** C C WORK IN DOUBLE PRECISION C IMPLICIT REAL*8(A-H,O-$) C C USE FUNCTION STATEMENTS FOR COMPUTATION OF RADII OF CURVATURE C RNORM(X)=A/DSQRT(1.0D0-E2*DSIN(X)**2) RMER(X)=A*(1.0D0-E2)/DSQRT(1.0D0-E2*DSIN(X)**2)**3 C C CHECK IF ELLIPSOID DEFINED BY A & B OR A & 1/F AND COMPUTE C ECCENTRICITY C B1=B IF(B1.LT.1000.)B1=A-A/B1 E2=(A*A-B1*B1)/A/A C C COMPUTE RADII OF CURVATURE A POINT I C RNI=RNORM(PHII) RMI=RMER(PHII) C C COMPUTE DPHI AND ITERATE UNTIL DESIRED ACCURACY IS MET C C COMPUTE SOME CONSTANTS C S2=SIJ**2 S3=SIJ**3 CA=DCOS(AIJ) SA=DSIN(AIJ) SA2=SA**2 SP=DSIN(PHII) SP2=SP**2 CP=DCOS(PHII) TP=DTAN(PHII) TP2=TP**2 ACC=1.0D-10 C C COMPUTE FIRST APPROXIMATION TO DPHI C DP1=SIJ/RNI*CA-S2/2.0D0/RNI**2*TP*SA2-S3/6.0D0/RNI**3*CA*SA2* .(1.0D0+3.0D0*TP2) C C COMPUTE NEW VALUE OF DPHI C 1 DP2=(SIJ*CA/RMI-S2*TP*SA2/2.0D0/RMI/RNI-S3*CA*SA2*(1.0D0+3.0D0*TP2 .)/6.0D0/RMI/RNI**2)*(1.0D0-(3.0D0*E2*SP*CP*DP1/(2.0D0-2.0D0*E2*SP2 .))) C C CHECK IF DESIRED ACCURACY MET C DIFF=DABS(DP2-DP1) IF(DIFF.LT.ACC)GO TO 2 DP1=DP2 C C RETURN TO COMPUTE A NEW VALUE OF DPHI C GO TO 1 C C COMPUTE LATITUDE OF POINT J C 2 PHIJ=PHII+DP2 C C COMPUTE DLAMDA C RNJ=RNORM(PHIJ) DL=SIJ/RNJ*SA/DCOS(PHIJ)*(1.0D0-S2/6.0D0/RNJ**2*(1.0D0-SA2/ .DCOS(PHIJ)**2)) C C COMPUTE LONGTITUDE OF J C XLONJ=XLONI+DL C C COMPUTE THE CONVERGENCE DAZIMUTH C PHIM=(PHII+PHIJ)/2.0D0 SPM=DSIN(PHIM) SDP=1.0D0/DCOS(DP2/2.0D0) DA=DL*SPM*SDP+DL**3/12.0D0*(SPM*SDP-SPM**3*SDP**3) C C COMPUTE THE BACKWARD AZIMUTH AJI AND RETURN C PI=3.141592653589793 AJI=AIJ+DA+PI IF(AJI.GT.(2.0D0*PI))AJI=AJI-2.0D0*PI RETURN END SUBROUTINE INVRSE(A,B,PHII,XLONI,PHIJ,XLONJ,AIJ,AJI,SIJ) C*********************************************************************** C* * C* SUBROUTINE INVRSE SOLVES THE INVERSE PROBLEM IN GEODESY, THAT * C* IS, GIVEN THE LATITUDE AND LONGTITUDE OF TWO POINTS I AND J, * C* COMPUTE THE FORWARD AND BACKWARD AZIMUTHS AIJ AND AJI AND THE * C* DISTANCE BETWEEN I AND J. * C* * C* THE EQUATIONS PROGRAMMED ARE 13, 15, 118, 115, 109(A), 116 * C* AND 117 'GEODETIC POSITION COMPUTATIONS' EXCEPT THAT IN ALL CASES * C* WHERE IT APPLIES, THE COMPUTATIONS ARE IN RADIANS RATHER THAN * C* SECONDS OF ARC. THESE EQUATIONS ARE KNOWN AS 'PUISSANT'S FORMULA * C* FOR SHORT LINES' WHICH GIVE A REPUTED ACCURACY OF 1 PPM AT 100KM * C* (IE. THEY ARE BASED ON A SPHERICAL APPROXIMATION TO THE ELLIPSOID) * C* * C* THE SOLUTION IS ITERATED UNTIL THE CHANGE IN DALPHA IS LESS * C* THAN 1*10-10 RADIANS AND THE CHANGE IN SIJ IS LESS THAN 0.1 MM. * C* * C* NOTE: THE EQUIVALENT DIRECT PROGRAM IS SUBROUTINE 'DIRECT' * C* * C* INPUT DATA: * C* A=SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID IN METRES * C* B=SEMI-MINOR AXIS OF REFERENCE ELLIPSOID OR THE * C* RECIPROCAL OF THE FLATTENING IN METRES OR UNITLESS * C* PHII=GEODETIC LATITUDE OF POINT I IN RADIANS * C* XLONI=GEODETIC LONGTITUDE OF POINT I IN RADIANS * C* PHIJ=GEODETIC LATITUDE OF POINT J IN RADIANS * C* XLONJ=GEODETIC LONGTITUDE OF POINT J IN RADIANS * C* * C* OUTPUT DATA: * C* AIJ=GEODETIC AZIMUTH FROM I TO J IN RADIANS * C* NOTE: AZIMUTH ON ELLIPSOID * C* AJI=GEODETIC AZIMUTH FROM J TO I IN RADIANS * C* NOTE: AZIMUTH ON ELLIPSOID * C* SIJ=DISTANCE FROM I TO J (J TO I) IN METRES * C* NOTE: DISTANCE ON ELLIPSOID * C* * C* ( C. CHAMBERLAIN NOV. 1974 ) * C* * C*********************************************************************** C C WORK IN DOUBLE PRECISION C IMPLICIT REAL*8(A-H,O-$) C C COMPUTE RADII OF CURVATURE IN FUNCTION STATEMENTS C RNORM(X)=A/DSQRT(1.0D0-E2*DSIN(X)**2) RMER(X)=A*(1.0D0-E2)/DSQRT(1.0D0-E2*DSIN(X)**2)**3 C C CHECK IF ELLIPSOID DEFINED BY A & B OR A & 1/F AND COMPUTE C ECCENTRICITY C B1=B IF(B1.LT.1000.0)B1=A-A/B1 E2=(A*A-B1*B1)/A/A C C COMPUTE THE RADII OF CURVATURE AT I AND J C RMI=RMER(PHII) RNI=RNORM(PHII) RNJ=RNORM(PHIJ) C C COMPUTE THE FIRST APPROXIMATION TO AIJ AND SIJ C ACC1=1.0D-10 ACC2=1.0D-4 RNJ2=RNJ**2 DP=PHIJ-PHII DL=XLONJ-XLONI TPI=DTAN(PHII) TPI2=TPI**2 SPJ=1.0D0/DCOS(PHIJ) SPJ2=SPJ**2 XT=1.0D0-3.0D0*DSIN(PHII)*DCOS(PHII)*DP/2.0D0/(1.0D0-E2*DSIN(PHII) .**2)*E2 AIJ1=DATAN2(DL*RNJ/SPJ,DP*RMI/XT) PI=3.141592653589793D0 PIE=PI/2.0D0 TSO=PI*2.0D0 IF(DABS(AIJ1).LT.ACC1)S1=DP*RMI/XT/DCOS(AIJ1) IF(DABS(AIJ1).LT.ACC1)GO TO 1 IF(DABS(AIJ1-PIE).LT.ACC1)S1= DL*RNJ/SPJ/DSIN(AIJ1) IF(DABS(AIJ1-PIE).LT.ACC1)GO TO 1 IF(DABS(AIJ1-PI).LT.ACC1)S1=DP*RMI/XT/DCOS(AIJ1) IF(DABS(AIJ1-PI).LT.ACC1)GO TO 1 IF(DABS(AIJ1+PIE).LT.ACC1)S1=DL*RNJ/SPJ/DSIN(AIJ1) IF(DABS(AIJ1+PIE).LT.ACC1)GO TO 1 S1=DP*RMI/DCOS(AIJ1)/XT C C COMPUTE THE NEXT APPROXIMATION TO AIJ AND SIJ C 1 EQ117=DL*RNJ/SPJ+S1**3*DSIN(AIJ1)/6.0D0/RNJ2-S1**3*DSIN(AIJ1)**3*S .PJ2/6.0D0/RNJ2 EQ116=RMI*DP/XT+S1**2*TPI*DSIN(AIJ1)**2/2.0D0/RNI+(S1**3*DCOS(AIJ1 .) *DSIN(AIJ1)**2*(1.0D0+3.0D0*TPI2)/6.0D0/RNI/RNI) AIJ2=DATAN2(EQ117,EQ116) IF(DABS(AIJ2).LT.ACC1)S2=EQ116/DCOS(AIJ2) IF(DABS(AIJ2).LT.ACC1)GO TO 3 IF(DABS(AIJ2-PIE).LT.ACC1)S2=EQ117/DSIN(AIJ2) IF(DABS(AIJ2-PIE).LT.ACC1)GO TO 3 IF(DABS(AIJ2-PI).LT.ACC1)S2=EQ116/DCOS(AIJ2) IF(DABS(AIJ2-PI).LT.ACC1)GO TO 3 IF(DABS(AIJ2+PIE ).LT.ACC1)S2=EQ117/DSIN(AIJ2) IF(DABS(AIJ2+PIE ).LT.ACC1)GO TO 3 S2=EQ116/DCOS(AIJ2) C C CHECK IF ACCURACY HAS BEEN MET C 3 DIFF1=DABS(AIJ2-AIJ1) DIFF2=DABS(S2-S1) IF(DIFF1.LT.ACC1.AND.DIFF2.LT.ACC2)GO TO 2 S1=S2 AIJ1=AIJ2 C C RETURN AND COMPUTE A NEW APPROXIMATION TO AIJ AND SIJ C GO TO 1 C C COMPUTE THE FINAL VALUES OF AIJ AND SIJ C 2 AIJ=AIJ2 IF(AIJ.GT.TSO)AIJ=AIJ-TSO IF(AIJ.LT.0.0)AIJ=AIJ+TSO IF(DABS(AIJ).LT.ACC1)SIJ=EQ116/DCOS(AIJ) IF(DABS(AIJ).LT.ACC1)GO TO 4 IF(DABS(AIJ-PIE).LT.ACC1)SIJ=EQ117/DSIN(AIJ) IF(DABS(AIJ-PIE).LT.ACC1)GO TO 4 IF(DABS(AIJ-PI).LT.ACC1)SIJ=EQ116/DCOS(AIJ) IF(DABS(AIJ-PI).LT.ACC1)GO TO 4 IF(DABS(AIJ+PIE ).LT.ACC1)SIJ=EQ117/DSIN(AIJ) IF(DABS(AIJ+PIE ).LT.ACC1)GO TO 4 SIJ=(EQ116/DCOS(AIJ)+EQ117/DSIN(AIJ))/2.0D0 C C COMPUTE THE CONVERGENCE DA C 4 PHIM=(PHII+PHIJ)/2.0D0 SPM=DSIN(PHIM) SDP=1.0D0/DCOS(DP/2.0D0) DA=DL*SPM*SDP+DL**3/12.0D0*(SPM*SDP-SPM**3*SDP**3) C C COMPUTE THE FINAL VALUE OF AJI AND RETURN C AJI=AIJ+DA+PI IF(AJI.GT.TSO)AJI=AJI-TSO RETURN END SUBROUTINE RADARC (A,I,J,S) C * SUBROUTINE 'RADARC' CONVERTS DOUBLE PRECISION RADIANS TO SINGLE* 6260 C * PRECISION DEGREES,MINUTES AND SECONDS * 6270 C ****************************************************************** 6280 REAL*8 A,SEC,AD,AJ,RHO/206264.8062470964D0/ 6290 IF(A)2,1,1 6300 1 SIGN=1. 6310 GO TO 3 6320 2 SIGN=-1. 6330 A=-A 6340 3 SEC=A*RHO+0.0005 6350 I=SEC/3600. 6360 AD=I 6370 J=SEC/60.-AD*60. 6380 AJ=J 6390 S=SEC-AD*3600.-AJ*60.-0.0005 6400 IF(I)5,4,5 6410 4 IF(J)6,7,6 6420 6 J=J*SIGN 6430 GO TO 10 6440 7 S=S*SIGN 6450 GO TO 10 6460 5 I=I*SIGN 6470 10 CONTINUE 6480 IF(SIGN.LT.0.)A=-A RETURN 6490 END 6500 SUBROUTINE ARCRAD (I,M,S,RADS) C *$**************************************************************** 1390 C * SUBROUTINE 'ARCRAD' CONVERTS SINGLE PRECISION DEGREES,MINUTES * 1400 C * AND SECONDS TO DOUBLE PRECISION RADIANS * 1410 C *$**************************************************************'' 1420 REAL*8 RADS,D,Q,SEC,DABS,RHO/57.29577951308233D0/ 1430 D=I 1440 Q=M 1450 SEC=S 1460 IF(I)1,2,3 1470 1 SIGN=-1. 1480 GO TO 7 1490 2 IF(M)4,5,3 1500 4 SIGN=-1. 1510 GO TO 7 1520 5 IF(S)6,99,3 1530 6 SIGN=-1. 1540 GO TO 7 1550 3 SIGN=1. 1560 7 RADS=((DABS(D)+DABS(Q/60.)+DABS(SEC/3600.))/RHO)*SIGN 1570 GO TO 999 1580 99 RADS=0.D0 1590 999 CONTINUE 1600 RETURN 1610 END 1620