IMPLICIT REAL*8(A-H,O-Z) DIMENSION C(4,4) ,SLL(4,4),C1(6,6),C2(7,7) ,B4(5,4),C3(5,5) DIMENSIONB3(6,4),B1(6,4),B2(7,6),PL(4,4) DIMENSION B5(1,3),VC(3,3) DO 99 III=1,3 PI=DARCOS(-1.0D0) RHO=6.48D5/PI A=6378206.4D0 B=6356583.8D0 C C C READ IN HEIGHTS OF TWO POINTS AND V-C INFORMATION C C READ(5,*)H1,H2 READ(5,*)SH1,SH2,SH3 C C C READ IN DISTANCE AND STANDARD DEVIATION C C READ(5,*)DIST READ(5,*)DS,SD C C C READ IN ASTRO AZIMUTH ,ZENITH DISTANCE AND DEFLECTIONS AND CONVERT C THEM TO RADIANS C C READ(5,*)I1,I2,SEC READ(5,*)SA CALL DEGRAD(I1,I2,SEC,AAZ12) READ(5,*)I1,I2,S CALL DEGRAD(I1,I2,S,ZENANG) READ(5,*)I1,I2,S ATE=S CALL DEGRAD(I1,I2,S,ETA) READ(5,*)I1,I2,S XIX=S CALL DEGRAD(I1,I2,S,XI) WRITE(6,65) C C C READ IN COORDINATES OF POINT 1 AND CONVERT TO RADIANS C READ IN V-C INFORMATION OF POINT 1 C C READ(5,*)I1,I2,S CALL DEGRAD(I1,I2,S,PHI1) READ(5,*)I1,I2,S CALL DEGRAD(I1,I2,S,AMB1) READ(5,*)SPI,SLI,SPL SD=DS+SD*DIST WRITE(6,159)H1,H2 WRITE(6,160)ATE,XIX C C C COMPUTE APPROXIMATE COORDINATES OF POINT 2 IN RADIANS C C C CALL APPROX(PHI1,AMB1,AAZ12,DIST,A,B,PHI2,AMB2) C C COMPUTATION OF CORRECTIONS.IN ALL CASES CORRECTIONS ARE ADDED C C ALP21=AAZ12+PI ALP12=AAZ12 30 CONTINUE PHI2P=PHI2 AMB2P=AMB2 C C COMPUTE CORRECTION TO ZENITH DISTANCE AND ADD TO OBSERVED ZENITH C DISTANCE C C CALL RADEG(ZENANG,IDZ1,IMZ1,SEZ1) CALL ZENCOR(ALP12,XI,ETA,CORREC) ZEN=CORREC*RHO ZENANG=ZENANG+CORREC CALL RADEG(ZENANG,IDZ,IMZ,SEZ) C C C COMPUTE GRAVITY CORRECTION C C C CALL GRACOR(PHI1,ETA,XI,ALP12,ZENANG,CORREC) GRA=CORREC*RHO ALP12=AAZ12+CORREC CALL RADEG(ALP12,ID1,IM1,SE1) C C COMPUTE SKEW NORMAL CORRECTION C C CALL SKECOR(H2,ALP12,PHI1,PHI2,CORREC,A,B) SKE=CORREC*RHO ALP12=ALP12+CORREC CALL RADEG(ALP12,ID2,IM2,SE2) C C C COMPUTE ELLIPSOIDAL DISTANCE FROM SPATIAL DISTANCE AND ASSOCIATED C V-C MATRIX C CALL ELLCOR(H1,H2,DIST,PHI1,PHI2,ALP12,EDIST,A,B,SD,SH1,SH2,SH3, *VD,B5,VC,ALP21) C C C C COMPUTE NORMAL SECTION CORRECTION C C C CALL NSGCOR(PHI1,PHI2,ALP12,EDIST,A,B,CORREC) CNS=CORREC*RHO ALP12=ALP12+CORREC CALL RADEG(ALP12,ID3,IM3,SE3) C C PRINT LATITUDE AND LONGITUDE OF POINT 1 ,THE APPROXIMATE COORDINATES C OF POINT 2 AND THE CORRECTED AZIMUTH C WRITE(6,108) CALL RADEG(PHI1,IDEG,IMIN,SEC) WRITE(6,105)IDEG,IMIN,SEC CALL RADEG(AMB1,IDEG,IMIN,SEC) WRITE(6,106)IDEG,IMIN,SEC CALL RADEG(PHI2,IDEG,IMIN,SEC) WRITE(6,103)IDEG,IMIN,SEC CALL RADEG(AMB2,IDEG,IMIN,SEC) WRITE(6,104)IDEG,IMIN,SEC CALL RADEG(ALP12,IDEG,IMIN,SEC) WRITE(6,100)IDEG,IMIN,SEC C C PRINT OUT CORRECTIONS C C WRITE(6,150) WRITE(6,161)IDZ1,IMZ1,SEZ1 WRITE(6,151)ZEN WRITE(6,162)IDZ,IMZ,SEZ WRITE(6,152)GRA WRITE(6,163)ID1,IM1,SE1 RO=6.48D5/PI WRITE(6,153)SKE WRITE(6,163)ID2,IM2,SE2 WRITE(6,154)CNS WRITE(6,163)ID3,IM3,SE3 WRITE(6,158) WRITE(6,155) CALL MOUTD(B5,1,1,3) WRITE(6,156) CALL MOUTD(VC,3,3,3) WRITE(6,102)EDIST WRITE(6,157)VD C C C PUISSANT'S DIRECT PROBLEM C C CALL ELLDIR(PHI1,AMB1,ALP12,EDIST,ALP21,PHI2P,AMB2P,A,B) C C C CHECK TO SEE IF INTIAL APPROXIMATIONS WERE WITHIN 1SECOND C OF FINAL ANSWER C C DIFF1=DABS(PHI2-PHI2P) DIFF2=DABS(AMB2-AMB2P) PHI2=PHI2P AMB2=AMB2P IF((DIFF1.GT.4.848D-6).OR .(DIFF2.GT.4.848D-6))GO TO 30 C C C PRINT OUT SOLUTION OF PUISSANT'S DIRECT PROBLEM C C C WRITE(6,165) CALL RADEG(PHI1,IDEG,IMIN,SEC) WRITE(6,105)IDEG,IMIN,SEC CALL RADEG(AMB1,IDEG,IMIN,SEC) WRITE(6,106)IDEG,IMIN,SEC CALL RADEG(PHI2,IDEG,IMIN,SEC) WRITE(6,103)IDEG,IMIN,SEC CALL RADEG(AMB2,IDEG,IMIN,SEC) WRITE(6,104)IDEG,IMIN,SEC CALL RADEG(ALP12,IDEG,IMIN,SEC) WRITE(6,100)IDEG,IMIN,SEC CALL RADEG(ALP21,IDEG,IMIN,SEC) WRITE(6,101)IDEG,IMIN,SEC WRITE(6,102)EDIST C C PUISSANT'S INVERSE PROBLEM C C C CALL ELLINV(PHI1,AMB1,PHI2,AMB2,ALP12,EDIST,ALP21,A,B) C C PRINT OUT SOLUTION OF PUISSANT'S INVERSE PROBLEM C C WRITE(6,111) CALL RADEG(PHI1,IDEG,IMIN,SEC) WRITE(6,105)IDEG,IMIN,SEC CALL RADEG(AMB1,IDEG,IMIN,SEC) WRITE(6,106)IDEG,IMIN,SEC CALL RADEG(PHI2,IDEG,IMIN,SEC) WRITE(6,103)IDEG,IMIN,SEC CALL RADEG(AMB2,IDEG,IMIN,SEC) WRITE(6,104)IDEG,IMIN,SEC CALL RADEG(ALP12,IDEG,IMIN,SEC) WRITE(6,100)IDEG,IMIN,SEC CALL RADEG(ALP21,IDEG,IMIN,SEC) WRITE(6,101)IDEG,IMIN,SEC WRITE(6,102)EDIST C C C ERROR PROPAGATION OF DIRECT PROBLEM C C CALL ERDML(PHI1,PHI2,ALP12,ALP21,EDIST,A,B,SPI,SLI,SPL,SA,VD,C, *AMB1,AMB2,B3,C3,B4,PL) DO 15 K=1,4 DO 15 J=1,4 SLL(K,J)=C(K,J) 15 CONTINUE C C C ERROR PROPAGATION OF INVERSE PROBLEM C C CALL ERIML(PHI1,AMB1,PHI2,AMB2,ALP12,ALP21,A,B,B1,C2,SLL,B2,C1) C C C PRINT OUT ERROR PROPAGATION - DIRECT CASE C C WRITE(6,65) 65 FORMAT('1') WRITE(6,115) WRITE(6,112) CALL MOUTD(B3,4,4,4) WRITE(6,114) CALL MOUTD(B4,5,5,4) WRITE(6,65) WRITE(6,124) CALL MOUTD(PL,4,4,4) WRITE(6,113) WRITE(6,120) CALL MOUTD(C,4,4,4) WRITE(6,121) CALL MOUTD(C3,5,5,5) C C C PRINT OUT ERROR PROPAGATION - INVERSE CASE C C WRITE(6,65) WRITE(6,116) WRITE(6,112) CALL MOUTD(B1,6,6,4) WRITE(6,114) CALL MOUTD(B2,7,7,6) WRITE(6,65) WRITE(6,113) WRITE(6,122) CALL MOUTD(C1,6,6,6) WRITE(6,123) CALL MOUTD(C2,7,7,7) 100 FORMAT('0','ALPHA-12 ='/'0','DEG',I4,5X,'MIN',I3,5X,'SEC',F10.4) 101 FORMAT('0','ALPHA-21 ='/'0','DEG',I4,5X,'MIN',I3,5X,'SEC',F10.4) 102 FORMAT('0','ELLIPSOID DISTANCE',3X,F15.4) 103 FORMAT('0','PHI-2 '/'0','DEG',I4,5X,'MIN',I3,5X,'SEC',F10.4) 104 FORMAT('0','LAMBDA-2 '/'0','DEG',I4,5X,'MIN',I3,5X,'SEC',F10.4) 105 FORMAT('0','PHI-1 '/'0','DEG',I4,5X,'MIN',I3,5X,'SEC',F10.4) 106 FORMAT('0','LAMBDA-1 '/'0','DEG',I4,5X,'MIN',I3,5X,'SEC',F10.4) 108 FORMAT('-',25X,'APPROXIMATIONS') 111 FORMAT('-','**** INVERSE PROBLEM****') 113 FORMAT('-',30X,' VARIANCE-COVARIANCE MATRIX') 112 FORMAT('-',30X,'FIRST DESIGN MATRIX') 114 FORMAT('-',30X,'SECOND DESIGN MATRIX') 115 FORMAT('-',30X,'ERROR PROPAGATION-DIRECT CASE') 116 FORMAT('-',30X,'ERROR PROPAGATION-INVERSE CASE') 120 FORMAT('-',5X,'V-C MATRIX FOR PHI1,RLAM1,PHI2,RLAM2, IN SECONDS *SQUARED') 121 FORMAT('-',5X,'V-C MATRIX FOR PHI1,RLAM1,PHI2,RLAM2,A21 IN SECONDS *SQUARED') 122 FORMAT(' ',5X,'V-C MATRIX FOR PHI1,RLAM1,PHI2,RLAM2,A12,A21 IN *SECONDS SQUARED') 123 FORMAT('-',5X,'V-C MATRIX FOR PHI1,RLAM1,PHI2,RLAM2,A12,A21,S IN *SECONDS SQUARED,METERS SQUARED AND SECOND METERS') 124 FORMAT('-',5X,'WEIGHT MATRIX FOR PHI1,RLAM1,A12,S IN SECONDS *SQUARED AND METERS SQUARED') 150 FORMAT('0',20X,'CORRECTIONS TO AZIMUTH') 151 FORMAT('0',25X,'ZENITH DISTANCE CORRECTION=',F7.3) 152 FORMAT('0',25X,'GRAVITY CORRECTION=',F7.3) 153 FORMAT('0',25X,'SKEW NORMAL CORRECTION=',F7.3) 154 FORMAT('0',25X,'NORMAL SECTION CORRECTION=',F7.3) 155 FORMAT('0',20X,'DESIGN MATRIX') 156 FORMAT('0',20X,'INPUT WEIGHT MATRIX') 157 FORMAT('0',20X,'STANDARD DEVIATION OF ELLIPSOIDAL DISTANCE=',F8.3) 158 FORMAT('0',10X,'MATRICES FOR ERROR PROPAGATION IN THE REDUCTION *OF THE SPATIAL DISTANCE TO ELLIPSOIDAL DISTANCE') 159 FORMAT('0',10X,'HEIGHT OF POINT 1 =',F12.4,10X,'HEIGHT OF POINT 2 *=',F12.4) 160 FORMAT('0',15X,'DEFLECTION COMPONENTS ARE',10X,'ETA =',F8.3,10X, *'XI =',F8.3) 161 FORMAT('0',20X,'MEASURED ZENITH DISTANCE',2X,I3,2X,I3,2X,F7.3) 162 FORMAT('0',20X,'CORRECTED ZENITH DISTANCE',2X,I3,2X,I3,2X,F7.3) 163 FORMAT('0',20X,'CORRECTED AZIMUTH',2X,I4,2X,I3,2X,F7.3) 165 FORMAT('0',20X,'DIRECT PROBLEM') 20 CONTINUE 99 CONTINUE STOP END SUBROUTINE APPROX(PHI1,AMB1,AAZ12,DIST,A,B,PHI2,AMB2) C C SUBROUTINE TO COMPUTE APPROXIMATE COORDINATES OF A SECOND PT. C C *** INPUT *** C C AMB1----GEODETIC LONGITUDE OF PT. 1 (RADIANS) C PHI1----GEODETIC LATITUDE OF PT.1 (RADIANS) C AAZ12----OBSERVED ASTONOMIC AZIMUTH (RADIANS) C DIST----SPATIAL DISTANCE (METRES) C A----SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID(METRES) C B----SEMI-MINOR AXIS OF REFERENCE ELLIPSOID(METRES) C C *** OUTPUT *** C C PHI2----GEODETIC LATITUDE OF PT.2 (RADIANS) C AMB2----GEODETIC LONGITUDE OF PT. 2 (RADIANS) C IMPLICIT REAL*8(A-H,O-Z) E=(A**2-B**2)/(A**2) RN1=A/((1.0D0-E*((DSIN(PHI1))**2))**.5D0) RM1=A*(1.0D0-E)/((1.0D0-E*((DSIN(PHI1))**2))**1.5D0) PHI2=PHI1+DIST*DCOS(AAZ12)/RM1 PHIM=(PHI1+PHI2)/2.0D0 AMB2=AMB1+DIST*DSIN(AAZ12)/(RN1*DCOS(PHIM)) RETURN END SUBROUTINE ZENCOR(AAZ12,XI,ETA,CORREC) C C SUBROUTINE TO COMPUTE ZENITH ANGLE CORRECTION C C ***INPUT*** C C ETA----PRIME VERTICAL DEFLECTION AT PT.1 (RADIANS) C XI----MERIDIAN DEFLECTION AT PT.1 (RADIANS) C AAZ12----OBSERVED ASTRO AZIMUTH FROM PT.1 TO PT.2 (RADIANS) C C ***OUTPUT*** C C CORREC----CORRECTIVE TERM TO BE ADDED TO OBSERVATION(RADIANS) C IMPLICIT REAL*8(A-H,O-Z) CORREC=XI*DCOS(AAZ12)+ETA*DSIN(AAZ12) RETURN END SUBROUTINE GRACOR(PHI1,ETA,XI,AAZ12,ZENANG,CORREC) C C SUBROUTINE COMPUTES COMPLETE CORRECTION C TO OBSERVED ASTRO AZIMUTH RESULTING C FROM GRAVITY DEFLECTIONS AND ZENITH DISTANCE C C ***INPUT*** C C PHI1----GEODETIC LATITUDE OF PT.1 (RADIANS) C ETA----PRIME VERTICAL DEFLECTION AT PT.1 (RADIANS) C XI----MERIDIAN DEFLECTION AT PT.1 (RADIANS) C AAZ12----OBSERVED ASTRO AZIMUTH FROM PT.1 TO PT.2 (RADIANS) C C ***OUTPUT*** C C CORREC----CORRECTIVE TERM TO BE ADDED TO OBSERVATION(RADIANS) C IMPLICIT REAL*8(A-H,O-Z) CORREC=(-ETA*DTAN(PHI1))-(XI*DSIN(AAZ12)- CETA*DCOS(AAZ12))/DTAN(ZENANG) RETURN END SUBROUTINE SKECOR(H2,ALP12,PHI1,PHI2,CORREC,A,B) C C SUBROUTINE COMPUTES SKEW NORMAL CORRECTION. C C ***INPUT*** C C H2----HEIGHT OF PT. 2 ABOVE ELLIPSOID SURFACE(METRES) C ALP12----AZIMUTH FROM PT.1 TO PT.2(RADIANS) C PHI1----GEODETIC LATITUDE OF PT.1 (RADIANS) C PHI2----GEODETIC LATITUDE OF PT.2 (RADIANS) C A----SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID(METRES) C B----SEMI-MINOR AXIS OF REFERENCE ELLIPSOID(METRES) C C ***OUTPUT*** C C CORREC----CORRECTION TO BE ADDED TO OBSERVATION(RADIANS) C IMPLICIT REAL*8(A-H,M-Z) PI=DARCOS(-1.0D0) P=1.0D0 E=(A**2-B**2)/A**2 M1=A*(1.0D0-E)/((1.0D0-E*((DSIN(PHI1))**2))**1.5) M2=A*(1.0D0-E)/((1.0D0-E*((DSIN(PHI2))**2))**1.5) MM=(M1+M2)/2.0D0 CORREC=P*H2*E*DSIN(ALP12)*DCOS(ALP12)*((DCOS(PHI2))**2)/MM RETURN END SUBROUTINE NSGCOR(PHI1,PHI2,ALPH12,S12,A,B,CORREC) C C SUBROUTINE COMPUTES THE NORMAL SECTION TO GEODESIC CORRCTION. C C ***INPUT*** C C PHI1----GEODETIC LATITUDE OF PT.1 (RADIANS) C PHI2----GEODETIC LATITUDE OF PT.2 (RADIANS) C ALP12----NORMAL SECTION AZIMUTH FROM PT.1 TO PT.2(RADIANS) C S12----ELLIPSOID DISTANCE BETWEEN PTS. 1-2.(METRES) C A----SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID(METRES) C B----SEMI-MINOR AXIS OF REFERENCE ELLIPSOID(METRES) C C ***OUTPUT*** C C CORREC----CORRECTION TO BE ADDED TO OBSERVATION(RADIANS) C IMPLICIT REAL*8(A-H,N-Z) PI=DARCOS(-1.0D0) P=-1.0D0 E=(A**2-B**2)/(A**2) PHIM=(PHI1+PHI2)/2.0D0 C=2.0D0*ALPH12 NM=A/((1.0D0-E*((DSIN(PHIM))**2))**.5D0) CORREC=P*E*(S12**2)*((DCOS(PHIM))**2)*DSIN(C)/(12.0D0*(NM**2)) RETURN END SUBROUTINE RADEG(RAD,IDEG,IMIN,SEC) C C SUBROUTINE TO CONVERT FROM RADIANS TO DEGREES C C ***INPUT*** C C RAD----ANGLE (RADIANS) C C ***OUTPUT*** C C IDEG----DEGREES(INTEGER) C IMIN----MINUTES (INTEGER) C SEC----SECONDS (REAL*8) C IMPLICIT REAL*8(A-H,O-Z) PI=DARCOS(-1.0D0) A=RAD*180.D0/PI IDEG=A B=A-IDEG A=B*60.0D0 IMIN=A B=A-IMIN SEC=B*60.0D0 RETURN END SUBROUTINE DEGRAD(IDEG,IMIN,SEC,RAD) C C SUBROUTINE TO CONVERT FROM DEGREES TO RADIANS C C ***INPUT*** C C IDEG----DEGREES(INTEGER) C IMIN----MINUTES (INTEGER) C SEC----SECONDS (REAL*8) C C ***OUTPUT*** C C RAD----ANGLE (RADIANS) C IMPLICIT REAL*8(A-H,O-Z) PI=DARCOS(-1.0D0) SIGN=1.0D0 IF(IDEG.LT.0)SIGN=-1.0D0 RAD=(IABS(IDEG)+IMIN/DFLOAT(60)+SEC/3600.D0)*PI/180.0D0 RAD=RAD*SIGN RETURN END SUBROUTINE ELLCOR(H1,H2,D,PHI1,PHI2,ALP12,SO,A,BB,SD,SH1,SH2,SH3, *VD,B,VC,ALP21) C C SUBROUTINE COMPUTES THE CORRECTIONS TO REDUCE SPATIAL DISTANCE TO C THE DISTANCE ON THE ELLIPSOID SURFACE. C C ***INPUT*** C C H1----HEIGHT OF PT. 1 ABOVE ELLIPSOID SURFACE(METRES) C H2----HEIGHT OF PT. 2 ABOVE ELLIPSOID SURFACE(METRES) C D----SPATIAL DISTANCE(METRES) C ALP12----GEODETIC AZIMUTH FROM PT.1 TO PT.2(RADIANS) C PHI1----GEODETIC LATITUDE OF PT.1 (RADIANS) C PHI2----GEODETIC LATITUDE OF PT.2 (RADIANS) C A----SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID(METRES) C BB----SEMI-MINOR AXIS OF REFERENCE ELLIPSOID(METRES) C SD-- STANDARD DEVIATION OF SPATIAL DISTANCE D--METERS C SH1--STANDARD DEVIATION OF HEIGHT AT POINT 1--METERS C SH2--STANDARD DEVIATION OF HEIGHT AT POINT 2--METERS C SH3--COVARIANCE OF HEIGHTS BETWEEN POINTS 1 AND 2--METERS SQUARED C C ***OUTPUT*** C C SO----ELLIPSOID DISTANCE(METRES) C B--DESIGN MATRIX C VD--STANDARD DEVIATION OF ELLIPSOIDAL DISTANCE BETWEEN POINTS 1 AND 2 C IMPLICIT REAL*8(A-H,O-Z) DIMENSION B(1,3),VC(3,3),BT(3,1) DIMENSION F(1,1),C(1,3) PI=DARCOS(-1.0D0) ALP21=ALP12+PI E=(A**2-BB**2)/(A**2) RN=A/((1.0D0-E*((DSIN(PHI1))**2))**.5D0) RM=A*(1.0D0-E)/((1.0D0-E*((DSIN(PHI1))**2))**1.5D0) R1=(RN*RM)/(RM*((DSIN(ALP12))**2)+RN*((DCOS(ALP12))**2)) RN=A/((1.0D0-E*((DSIN(PHI2))**2))**.5D0) RM=A*(1.0D0-E)/((1.0D0-E*((DSIN(PHI2))**2))**1.5D0) R2=(RN*RM)/(RM*((DSIN(ALP21))**2)+RN*((DCOS(ALP21))**2)) R=(R1+R2)/2.0D0 DD=(((D**2)-((H1-H2)**2))/((1.0D0+H1/R)*(1.0D0+H2/R)))**.5D0 SO=2.0D0*R*DARSIN(DD/(2.0D0*R)) P1=1.D0+H1/R P2=1.D0+H2/R DO 1 I=1,3 DO 1 J=1,3 VC(I,J)=0.0D0 1 CONTINUE VC(1,1)=SD**2 VC(2,2)=SH1**2 VC(2,3)=SH3 VC(3,2)=VC(2,2) VC(3,3)=SH2**2 B(1,1)=-1.D0*D/(DD*P1*P2) DH=H2-H1 DH2=DH*DH B(1,2)=(.5D0/DD)*((-1.D0*((D**2.D0)-DH2)/(R*(P1**2.D0)*P2)) *+(2.*DH/(P1*P2))) B(1,3)=(.5D0/DD)*((-1.D0*((D**2.D0)-DH2)/(R*(P2**2.D0)*P1)) *-(2.*DH/(P1*P2))) DO 2 I=1,3 BT(I,1)=B(1,I) 2 CONTINUE CALL MMULD(C,1,B,1,VC,3,1,3,3) CALL MMULD(F,1,C,1,BT,3,1,3,1) VD=DSQRT(F(1,1)) RETURN END SUBROUTINE ELLDIR(PHI,AMB,ALP,DIS,ALP2,PHI2,AMB2,A,B) C C SUBROUTINE TO DO DIRECT PROBLEM ON THE ELLIPSOID USING C PUISSANTS FORMULAE C C ***INPUT*** C PHI----GEODETIC LATITUDE OF PT.1 (RADIANS) C AMB----GEODETIC LONGITUDE OF PT. 1 (RADIANS) C DIS----GIVEN ELLIPSOID DISTANCE(METRES) C ALP----GIVEN GEODETIC AZIMUTH FROM PT.1 TO PT.2 (RADIANS) C A----SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID(METRES) C B----SEMI-MINOR AXIS OF REFERENCE ELLIPSOID(METRES) C C ***OUTPUT*** C C PHI2----GEODETIC LATITUDE OF PT.2 (RADIANS) C AMB2----GEODETIC LONGITUDE OF PT. 2 (RADIANS) C ALP2----COMPUTED GEODETIC AZIMUTH FROM PT.2 TO PT.1 (RADIANS) C IMPLICIT REAL*8(A-H,O-Z) REAL*8 N1,N2,M1 WRITE(6,5) E=(A**2-B**2)/(A**2) PI=DARCOS(-1.0D0) RO=6.48D5/PI N1=A/((1.0D0-E*((DSIN(PHI))**2))**.5D0) M1=A*(1.0D0-E)/((1.0D0-E*((DSIN(PHI))**2))**1.5D0) DP1=DIS/N1*DCOS(ALP) DP2=(DIS**2)*DTAN(PHI)*((DSIN(ALP))**2)/(2.0D0*(N1**2)) DP3=(DIS**3)*DCOS(ALP)*((DSIN(ALP))**2)*(1.0D0+((DTAN(PHI))**2) C*3.0D0)/(6.0D0*(N1**3.0D0)) DP=(DP1-DP2-DP3)*RO DO 1 JJ=1,10 DP1=DIS*DCOS(ALP)/M1 DP2=(DIS**2)*DTAN(PHI)*((DSIN(ALP))**2)/(2.0D0*M1*N1) DP3=(DIS**3)*DCOS(ALP)*((DSIN(ALP))**2)*(1.0D0+3.0D0* F((DTAN(PHI))**2))/(6.0D0*M1*(N1**2)) DS=RO*(DP1-DP2-DP3) DP4=(1.0D0-3.0D0*E*DSIN(PHI)*DCOS(PHI)/(2.0D0*(1.0D0-E* ^((DSIN(PHI))**2)))*DP/RO) DS=DS*DP4 PHI2=PHI+DS/RO DIFF=DP-DS DIFF =DIFF/RO WRITE(6,3)DIFF CALL RADEG(DIFF,ID,IM,SE) WRITE(6,4)ID,IM,SE 4 FORMAT('0',20X,'DELTA PHI =',2X,I3,2X,I3,2X,F9.5) DIFF=DIFF*RO 3 FORMAT('0',20X,'DELTA PHI=',F12.9) 5 FORMAT('0',20X,'ITERATIONS IN PUISSANTS DIRECT PROBLEM') IF(DABS(DIFF).LT.1.0D-5)GO TO 2 DP=DS 1 CONTINUE 2 CONTINUE N2=A/((1.0D0-E*((DSIN(PHI2))**2))**.5D0) DL1=DIS*DSIN(ALP)/(N2*DCOS(PHI2)) DL2=1.0D0-(DIS**2)*(1.0D0-((DSIN(ALP))**2)/((DCOS(PHI2))**2)) B/(6.0D0*(N2**2)) DL=RO*DL1*DL2 AMB2=AMB+DL/RO PM=(PHI+PHI2)/2.0D0 DA1=DL/RO*DSIN(PM)/(DCOS(DS/(2.0D0*RO))) DA2=((DL/RO)**3)/12.0D0*(DSIN(PM)/(DCOS(DS/(2.0D0*RO))) N-((DSIN(PM))**3)/((DCOS(DS/(2.0D0*RO)))**3)) DA=(DA1+DA2)*RO ALP2=ALP+PI+DA/RO IF(ALP.GT.PI)ALP2=ALP-PI+DA/RO IF(JJ.GT.9)WRITE(6,100) 100 FORMAT('0','WARNING---MAXIMUM # OF ITERATIONS PERFORMED IN ELLDIR CSUBROUTINE---CHECK FOR ERRORS') RETURN END SUBROUTINE ELLINV(PHI1,AMB1,PHI2,AMB2,ALP,DIS,ALP2,A,B) IMPLICIT REAL*8(A-H,O-Z) REAL*8 M1,N1,N2 C C SUBROUTINE TO DO INVERSE PROBLEM ON THE ELLIPSOID USING C PUISSANTS FORMULAE C C ***INPUT*** C C PHI1----GEODETIC LATITUDE OF PT.1 (RADIANS) C AMB1----GEODETIC LONGITUDE OF PT. 1 (RADIANS) C PHI2----GEODETIC LATITUDE OF PT.2 (RADIANS) C AMB2----GEODETIC LONGITUDE OF PT. 2 (RADIANS) C A----SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID(METRES) C B----SEMI-MINOR AXIS OF REFERENCE ELLIPSOID(METRES) C C ***OUTPUT*** C C DIS----ELLIPSOID DISTANCE(METRES) C ALP----GEODETIC AZIMUTH FROM PT.1 TO PT.2 (RADIANS) C ALP2----COMPUTED GEODETIC AZIMUTH FROM PT.2 TO PT.1 (RADIANS) C WRITE(6,5) 5 FORMAT('0',20X,'ITERATIONS IN PUISSANTS INVERSE PROBLEM') PI=DARCOS(-1.0D0) RO=6.48D5/PI E=(A**2-B**2)/(A**2) N1=A/((1.0D0-E*((DSIN(PHI1))**2))**.5D0) M1=A*(1.0D0-E)/((1.0D0-E*((DSIN(PHI1))**2))**1.5D0) N2=A/((1.0D0-E*((DSIN(PHI2))**2))**.5D0) DP=PHI2-PHI1 DL=AMB2-AMB1 F1=DP*M1/(1.0D0-3.0D0*E*DSIN(PHI1)*DCOS(PHI1)*DP/(2.0D0*(1.0D0 D-E*((DSIN(PHI1))**2)))) F2=DL*N2*DCOS(PHI2) ALP=DATAN2(F2,F1) IF(ALP.LT.0.0D0)ALP=2.0D0*PI+ALP AAA=ALP DIS=F1/DCOS(ALP) DO 1 JJ=1,10 F11=(DIS**2)*DTAN(PHI1)*((DSIN(ALP))**2)/(2.0D0*N1) F111=(DIS**3)*DCOS(ALP)*((DSIN(ALP))**2)*(1.0D0+3.0D0*((DTAN( NPHI1))**2))/(6.0D0*(N1**2)) F22=(DIS**3)*DSIN(ALP)/(6.0D0*(N2**2)) F222=(DIS**3)*((DSIN(ALP))**3)/(6.0D0*(N2**2)*((DCOS(PHI2))**2)) F=F1+F11+F111 FP=F2+F22-F222 ALP=DATAN2(FP,F) IF(ALP.LT.0.0D0)ALP=2.0D0*PI+ALP DIS=F/DCOS(ALP) DIFF=ALP-AAA WRITE(6,3)DIFF 3 FORMAT('0',20X,'AZIMUTH DIFFERENCE=',F12.9) CALL RADEG(DIFF,ID,IM,SE) WRITE(6,4)ID,IM,SE 4 FORMAT('0',20X,'AZIMUTH DIFFERENCE=',2X,I4,2X,I3,2X,F9.5) CALL RADEG(ALP,IDEG,IMIN,SEC) WRITE(6,6)IDEG,IMIN,SEC 6 FORMAT('0',20X,'AZIMUTH =',2X,I4,2X,I3,2X,F7.3) WRITE(6,7)DIS 7 FORMAT('0',20X,'ELLIPSOID DISTANCE=',F15.4) IF(DABS(DIFF).LT.1.0D-9)GO TO 10 AAA=ALP 1 CONTINUE 10 CONTINUE PM=(PHI1+PHI2)/2.0D0 DP=DP*RO DL=DL*RO DA1=DL/RO*DSIN(PM)/(DCOS(DP/(2.0D0*RO))) DA2=((DL/RO)**3)/12.0D0*(DSIN(PM)/(DCOS(DP/(2.0D0*RO))) N-((DSIN(PM))**3)/((DCOS(DP/(2.0D0*RO)))**3)) DA=(DA1+DA2)*RO ALP2=ALP+PI+DA/RO IF(ALP.GT.PI)ALP2=ALP-PI+DA/RO IF(JJ.GT.9)WRITE(6,100) 100 FORMAT('0','WARNING---MAXIMUM # OF ITERATIONS PERFORMED IN ELLINV CSUBROUTINE---CHECK FOR ERRORS') RETURN END SUBROUTINE INVCOR(H1,H2,SD,PHI1,PHI2,ALP12,EDIST,A,B) C C SUBROUTINE TO COMPUTE SPATIAL DISTANCE GIVEN THE ELLIPSOID C DISTANCE AND NECESSARY PARAMETERS C C ***INPUT*** C C H1----HEIGHT OF PT. 1 ABOVE ELLIPSOID SURFACE(METRES) C H2----HEIGHT OF PT. 2 ABOVE ELLIPSOID SURFACE(METRES) C PHI1----GEODETIC LATITUDE OF PT.1 (RADIANS) C PHI2----GEODETIC LATITUDE OF PT.2 (RADIANS) C ALP12----GEODETIC AZIMUTH FROM PT.1 TO PT.2(RADIANS) C EDIST----ELLIPSOID DISTANCE(METRES) C A----SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID(METRES) C B----SEMI-MINOR AXIS OF REFERENCE ELLIPSOID(METRES) C C ***OUTPUT*** C C SD----SPATIAL DISTANCE(METRES) C IMPLICIT REAL*8(A-H,O-Z) PI=DARCOS(-1.0D0) ALP21=ALP12+PI E=(A**2-B**2)/(A**2) RN=A/((1.0D0-E*((DSIN(PHI1))**2))**.5D0) RM=A*(1.0D0-E)/((1.0D0-E*((DSIN(PHI1))**2))**1.5D0) R1=(RN*RM)/(RM*((DSIN(ALP12))**2)+RN*((DCOS(ALP12))**2)) RN=A/((1.0D0-E*((DSIN(PHI2))**2))**.5D0) RM=A*(1.0D0-E)/((1.0D0-E*((DSIN(PHI2))**2))**1.5D0) R2=(RN*RM)/(RM*((DSIN(ALP21))**2)+RN*((DCOS(ALP21))**2)) R=(R1+R2)/2.0D0 DD=2.0D0*R*DSIN(EDIST/(2.0D0*R)) SD=((DD**2)*(1.0D0+H1/R)*(1.0D0+H2/R)+(H1-H2)**2)**.5D0 RETURN END SUBROUTINE ERDML(PHI1,PHI2,A12,A21,S,A,BB,SPI,SLI,SPL,SA,SD,C, *RLAM1,RLAM2,B,C3,B4,SL) C ***************************************************************** C * THIS SUBROUTINE COMPUTES THE ERROR PROPAGATION IN THE * C * * C * DIRECT CASE FOR GUASSIAN MEAN LATITUDE * C * * C * * C * INPUT PARAMETERS * C * * C * PHI1- ELLIPSOIDAL LATITUDE POINT 1 - RADIANS * C * PHI2- ELLIPSOIDAL LATITUDE POINT 2 - RADIANS * C * A12- ELLIPSOIDAL AZIMUTH FROM POINT 1 TO 2 - RADIANS * C * A21- ELLIPSOIDAL AZIMUTH FROM POINT 2 TO 1 - RADIANS * C * S - ELLIPSOIDAL DISTANCE BETWEEN TWO POINTS - METERS * C * A,B - SEMI-MAJOR AND SEMI-MINOR AXIS OF THE REFERENCE * C * ELLIPSOID - METERS * C * SPI - VARIANCE OF PHI1 * C * SLI - VARIANCE OF RLAM1 * C * SPL - COVARIANCE BETWEEN PHI1 AND RLAM1 * C * SA - VARIANCE OF AZIMUTH A12 * C * SD- STANDARD DEVIATION OF ELLIPSOIDAL DISTANCE S * C * * C * * C * OUTPUT * C * * C * B- FIRST DESIGN MATRIX - 4*4 * C * B4- SECOND DESIGN MATRIX - 5*4 * C * SL- WEIGHT MATRIX - 4*4 SECONDS AND METERS SQUARED * C * C- FIRST V-C MATRIX -4*4 SECONDS SQUARED * C * C3 - SECOND V-C MATRIX - 5*5 SECONDS SQUARED * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C * * C ***************************************************************** C IMPLICIT REAL*8(A-H,K-Z) DIMENSION B(4,4),BT(4,4),SL(4,4),SS(4,4),C(4,4),B4(5,4),BS(5,4), *C3(5,5),BT4(4,5) DO 10 I=1,4 DO 10 J=1,4 B(I,J)=0.D0 10 CONTINUE DO 20 I=1,4 DO 20 J=1,4 SL(I,J)=0.D0 20 CONTINUE PI=DARCOS(-1.D0) RHO=6.48D5/PI RHO2=RHO*RHO C C FORMATION OF INPUT WEIGHT MATIX C SL(1,1)=SPI/RHO2 SL(1,2)=SPL/RHO2 SL(2,1)=SPL/RHO2 SL(2,2)=SLI /RHO2 SL(3,3)=SA /RHO2 SL(4,4)=SD**2.D0 PHIM=(PHI1+PHI2)/2.D0 DLAM=RLAM2-RLAM1 ALM=A12+(DLAM*DSIN(PHIM)/2.D0) A2=A*A B2=BB*BB E2=(A2-B2)/A2 CAM=DCOS(ALM) CAM2=CAM*CAM CP1=DCOS(PHI1) CP2=DCOS(PHI2) SP1=DSIN(PHI1) SP12=SP1*SP1 SP2=DSIN(PHI2) SP22=SP2*SP2 PHIM=(PHI1+PHI2)/2.D0 DPH=PHI2-PHI1 DLAM=RLAM2-RLAM1 ES1=(1.D0-E2*SP12) ES2=1.D0-E2*SP22 AN1=A/(ES1**.5D0) AN2=A/(ES2**.5D0) ANM=(AN1+AN2)/2.D0 AM1=A*(1.D0-E2)/(ES1**1.5D0) AM2=A*(1.D0-E2)/(ES2**1.5D0) AMM=(AM1+AM2)/2.D0 CM=DCOS(PHIM) SM=DSIN(PHIM) AN12=AN1*AN1 AM12=AM1*AM1 ANM2=ANM*ANM P1=1.D0/(1.D0+DABS((DLAM/DPH)*(ANM*CM/AMM))**2.D0) Q1=DLAM/(DPH*AMM) R1=ANM*CM/DPH R2=(AM1*E2*SP1*CP1*CM)/(2.D0*(1.D0-E2)) R3=(1.5D0*ANM*CM*AM1*E2*SP1*CP1)/(AMM*ES1) R4=-.5D0*ANM*SM R5=DLAM*CM/4.D0 PP1=P1*(Q1*(R1+R2-R3+R4))-R5 R22=(AM2*E2*SP2*CP2*CM)/(2.D0*(1.D0-E2)) R33=(1.5D0*ANM*CM*AM2*E2*SP2*CP2)/(AMM*ES2) PP3=P1*(Q1*(-1.D0*R1+R22-R33+R4))-R5 Q2=(ANM*CM)/(DPH*AMM) PP2=P1*(-1.D0*Q2)+SM/2.D0 PP4=P1*Q2-SM/2.D0 SAM=DSIN(ALM) SM2=SM*SM CM=DCOS(PHIM) CM2=CM*CM AMM2=AMM*AMM CC=S*SAM*DLAM*CM/(AMM*4.D0) B(1,1)=1.D0 B(2,2)=1.D0 B(3,1)=1.-(3.*S*CAM*AM1*AN12*E2*SP1*CP1)/(A2*2.D0*AMM2) -CC B(3,2)=(S*SAM*SM)/(2.D0*AMM) B(3,3)=-1.D0*S*SAM/( AMM) B(3,4)=CAM/AMM AME=(AM1*E2)/(1.-E2) B(4,1)=(S*SAM*(ANM*SM-AME*SP1*CP1*CM))/(2.D0*ANM2*CM2)+(S*CAM*DLAM */(4.D0*ANM)) B(4,2)=1.D0-(S*CAM*SM/(ANM*CM*2.D0)) B(4,3)=S*CAM/( ANM*CM) B(4,4)=SAM/(ANM*CM) CALL TRNSD(BT,4,B,4,4,4) CALL MMULD(SS,4,B,4,SL,4,4,4,4) CALL MMULD(C,4,SS,4,BT,4,4,4,4) TM=DTAN(PHIM) DO 35 I=1,5 DO 35 J=1,4 B4(I,J)=0.D0 35 CONTINUE DO 36 I=1,4 B4(I,I)=1.D0 36 CONTINUE B4(5,1)=PP1+(DLAM*CM/2.D0) B4(5,2)=PP2-SM B4(5,3)=PP3+(DLAM*CM/2.D0) B4(5,4)=PP4+SM CALL TRNSD(BT4,4,B4,5,5,4) CALL MMULD(BS,5,B4,5,C,4,5,4,4) CALL MMULD(C3,5,BS,5,BT4,4,5,4,5) DO 40 I=1,4 DO 40 J=1,4 C(I,J)=C(I,J)*RHO2 40 CONTINUE DO 30 I=1,5 DO 30 J=1,5 C3(I,J)=C3(I,J)*RHO2 30 CONTINUE DO 45 I=1,3 DO 45 J=1,3 SL(I,J)=SL(I,J)*RHO2 45 CONTINUE RETURN END SUBROUTINE ERIML(PHI1,RLAM1,PHI2,RLAM2,A12,A21,A,BB,B ,C2,SL,B2,C1 *) C ***************************************************************** C * * C * * C * THIS SUBROUTINE COMPUTES THE ERROR PROPAGATION IN THE * C * * C * INVERSE CASE FOR GUASSIAN MEAN LATITUDE * C * * C * THE ERROR PROPAGATION IS DONE TWO STEPS . THE FIRST * C * * C * V-C MATRIX IS COMPUTED FOR PHI1,RLAM1,PHI2,RLAM2,A12,A21 * C * * C * THIS V-C MATRIX IS USED AS A SIGMA L MATRIX IN THE NEXT * C * * C * STEP. THE V-C MATRIX COMPUTED IN THE SECOND STEP IS THE * C * * C * SAME AS THE FIRST EXCEPT IT CONTAINS A TERM FOR THE * C * * C * DISTANCE. * C * * C * * C * INPUT PARAMETERS * C * * C * PHI1- ELLIPSOIDAL LATITUDE POINT 1 - RADIANS * C * PHI2- ELLIPSOIDAL LATITUDE POINT 2 - RADIANS * C * A12- ELLIPSOIDAL AZIMUTH FROM POINT 1 TO 2 - RADIANS * C * A21- ELLIPSOIDAL AZIMUTH FROM POINT 2 TO 1 - RADIANS * C * S - ELLIPSOIDAL DISTANCE BETWEEN TWO POINTS - METERS * C * A,B - SEMI-MAJOR AND SEMI-MINOR AXIS OF THE REFERENCE * C * ELLIPSOID - METERS * C * SL- INPUT COVARIANCE MATRIX FOR PHI1,RLAM1,PHI2,RLAM2 * C * - SECONDS SQUARED C * * C * * C * OUTPUT * C * * C * B- FIRST DESIGN MATRIX-6*4 * C * B2-SECOND DESIGN MATRIX - 7*6 C * C1- V-C MATRIX FOR PHI1,RLAM1,PHI2,RLAM2,A12,A21 - * C * SECONDS SQUARED * C * C2- VARIANCE - COVARIANCE MATRIX FOR PHI1,RLAM1,PHI2, * C * RLAM2,A12,A21,S - SECONDS SQUARED,METERS SQUARED C * * C * * C * * C * * C * * C * * C * * C * * C ***************************************************************** C IMPLICIT REAL*8(A-H,O-Z) DIMENSION SL(4,4),B(6,4),BT(4,6),BS(6,4),C1(6,6),B2(7,6),B2T(6,7) *,BS2(7,6),C2(7,7) PI=DARCOS(-1.D0) RHO=6.48D5/PI RHO2=RHO*RHO A2=A*A BB2=BB*BB E2=(A2-BB2)/A2 CP1=DCOS(PHI1) CP2=DCOS(PHI2) SP1=DSIN(PHI1) SP12=SP1*SP1 SP2=DSIN(PHI2) SP22=SP2*SP2 PHIM=(PHI1+PHI2)/2.D0 DPH=PHI2-PHI1 DLAM=RLAM2-RLAM1 ES1=(1.D0-E2*SP12) ES2=1.D0-E2*SP22 AN1=A/(ES1**.5D0) AN2=A/(ES2**.5D0) ANM=(AN1+AN2)/2.D0 AM1=A*(1.D0-E2)/(ES1**1.5D0) AM2=A*(1.D0-E2)/(ES2**1.5D0) AMM=(AM1+AM2)/2.D0 CM=DCOS(PHIM) SM=DSIN(PHIM) P1=1.D0/(1.D0+DABS((DLAM/DPH)*(ANM*CM/AMM))**2.D0) Q1=DLAM/(DPH*AMM) R1=ANM*CM/DPH R2=(AM1*E2*SP1*CP1*CM)/(2.D0*(1.D0-E2)) R3=(1.5D0*ANM*CM*AM1*E2*SP1*CP1)/(AMM*ES1) R4=-.5D0*ANM*SM R5=DLAM*CM/4.D0 PP1=P1*(Q1*(R1+R2-R3+R4))-R5 R22=(AM2*E2*SP2*CP2*CM)/(2.D0*(1.D0-E2)) R33=(1.5D0*ANM*CM*AM2*E2*SP2*CP2)/(AMM*ES2) PP3=P1*(Q1*(-1.D0*R1+R22-R33+R4))-R5 Q2=(ANM*CM)/(DPH*AMM) PP2=P1*(-1.D0*Q2)+SM/2.D0 PP4=P1*Q2-SM/2.D0 DO 10 I=1,6 DO 10 J=1,4 B(I,J)=0.D0 10 CONTINUE C C COMPUTING THE FIRST DESIGN C B(1,1)=1.D0 B(2,2)=1.D0 B(3,3)=1.D0 B(4,4)=1.D0 B(5,1)=PP1 B(5,2)=PP2 B(5,3)=PP3 B(5,4)=PP4 B(6,1)=PP1+(DLAM*CM/2.D0) B(6,2)=PP2-SM B(6,3)=PP3+(DLAM*CM/2.D0) B(6,4)=PP4+SM DO 20 I=1,4 DO 20 J=1,4 SL(I,J)=SL(I,J)/RHO2 20 CONTINUE C C COMPUTING THE FIRST V-C MATRIX-B*SLL*BT C CALL TRNSD(BT,4,B,6,6,4) CALL MMULD(BS,6,B,6,SL,4,6,4,4) CALL MMULD(C1,6,BS,6,BT,4,6,4,6) C C COMPUTING THE SECOND DESIGN MATRIX C DO 31 I=1,7 DO 31 J=1,6 B2(I,J)=0.D0 31 CONTINUE DO 32 I=1,6 B2(I,I)=1.D0 32 CONTINUE CA=DCOS((A12+A21-PI )/2.D0) EA=(1.5D0*AM1*E2*SP1*CP1)/((ES1 )*CA)*DPH B2(7,1)=-1.D0*AMM/CA+EA B2(7,3)=AMM/CA+((1.5D0*AM2*E2*SP2*CP2)/(ES2*CA))*DPH B2(7,2)=0.D0 B2(7,4)=0.D0 B2(7,5)=(DPH*AMM/(2.D0*(DABS(CA)**2.D0)))*DSIN((A12+A21-PI )/2 *.D0) B2(7,6)=B2(7,5) C C COMPUTING THE SECOND V-C MATRIX - B2*C1*BT2 C CALL TRNSD(B2T,6,B2,7,7,6) CALL MMULD(BS2,7,B2,7,C1,6,7,6,6) CALL MMULD(C2,7,BS2,7,B2T,6,7,6,7) C C CHANGING MATRIX INTO SECONDS DO 30 I=1,6 DO 30 J=1,6 C1(I,J)=C1(I,J)*RHO2 30 CONTINUE DO 33 I=1,7 DO 33 J=1,7 C2(I,J)=C2(I,J)*RHO 33 CONTINUE DO 34 I=1,6 DO 34 J=1,6 C2(I,J)=C2(I,J)*RHO 34 CONTINUE C2(7,7)=C2(7,7)/RHO RETURN END