C C**********E******************************************************************** C* * C*** PROGRAM 'L O O P G C'*** BY: M. M. NASSAR, MAY 1976. -***8 C**********E******************************************************************** C* * C* THE PURPOSE OF THIS COMPUTER PACKAGE IS TO COMPUTE THE ACCUMULATED GRAVITY * C* CORRECTION ALONG LEVELLING LINES OR LOOPS, COMPARE THEM NUMERICALLY WITH THE* C* THE CORRESPONDING STANDARD DEVIATION OF PRECISE SPIRIT LEVELLING IN A * C* TABULAR FORM , FOR THE PURPOSE OF COMPARISON ,ANALYSIS AND PLOTTING OF THE * C* OBTAINED RESULTS-FOR STUDYING THE INFLUNCE OF NEGLECTING THE GRAVITY * C* ANOMALIES ON HEIGHT COMPUTATIONS FROM SPIRIT LEVELLING. * C* * C* THREE SYSTEMS OF HEIGHTS ARE CONSIDERED * C* H E L M E R T , D Y N A M I C & V I G N A L * C* * C* FOR DETAILS ABOUT THE MATHEMATICAL BASIS AND PURPOSE OF SUCH AN INVESTIGATION C* -SEE: NASSAR'S PH.D. THESIS, 'GRAVITY FIELD AND LEVELLED HEIGHTS IN CANADA' * C* ,SURVEYING ENGINEERING DEPT., U.N.B., FEB. 1977. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-***-*-* C* NOTE:- THE DOCUMENTATION OF THE ENTIRE PROGRAM PACKAGE IS DONE WITHIN * C* THE PROGRAM SEGMENTS VIA SELF-EXPLAINED COMMENT CARDS. C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-***-*-* C**********E******************************************************************** C**********E******************************************************************** C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION PHI(100),ALONG(100),BMN(100) DIMENSION HT(100),HTACC(100),DGFA(100),DGACC(100) DIMENSION ACCDIS(100),ACCSTD(100),ACCHGC(100),ACCHSD(100) DIMENSION ACCVGC(100),ACCVSD(100),ACCDGC(100),ACCDSD(100) DIMENSION TITLE(20) DOUBLE PRECISION DSIN,DCOS,DABS,DSQRT IPAGE=49 ISKIP=IPAGE-1 PI=3.141592653589793D 0 RHO=206264.8062471D 0 RHODEG=180.D 0/PI STDEVN=1.33D 0 C C******************************************************************************* 555 CONTINUE C READ-IN THE IDENTIFICATION TITLE DESCRIBING THE INPUTTED LDVELLING-ROUTE: C THIS WILL BE PRINTED ON EACH PAGE OF THE CORRESPONDING OUTPUT. READ(5,86,END=777) (TITLE(I),I=1,20) 86 FORMAT(20A4) C C READ-IN THE TOTAL NUMBER OF BENCH-MARKS FORMING THE INPUTTED LEVELLING LINE C NOTE:- IF A CLOSED LOOP IS INPUTTED ,THE FIRST (LAST) BENCH-MARK HAS TO C BE COUNTED TWICE.- BECAUSE IT WILL BE TREATED AS TWO BENCH-MARKS. READ(5,41) NBM 41 FORMAT(I4) DO 10 I=1,NBM READ(5,5) BMN(I),PHI(I),ALONG(I),HT(I),HTACC(I),DGFA(I),DGACC(I) 5 FORMAT(6X,A4,2X,F8.5,1X,F9.5,2X,F8.2,4X,F6.2,3X,F7.2,4X,F6.2) PHIDEG=PHI(I) PHI(I)=PHI(I)/RHODEG ALONG(I)=ALONG(I)/RHODEG C IF(DGACC(I).EQ.0.0) DGACC(I)=DSQRT(0.0025D0+(0.3086*HTACC(I))**2) 10 CONTINUE C C PRINT-OUT THE TABLE OF INPUT INFORMATION AT EACH BENCH MARK WRITE(6,73) WRITE(6,87) (TITLE(I),I=1,20) 87 FORMAT(15X,20A4,/) WRITE(6,7) 7 FORMAT(30X,'TABLE OF GIVEN INFORMATION AT BENCH MARKS ',/,30 1X,'========================================='/) WRITE(6,9) 9 FORMAT(5X,' B.MARK',7X,'LATITUDE',10X,'LONGITUDE',10X,'HEIGHT',3X, 1'HT.ACCURACY',2X,'FREE-AIR ANOM',2X,'F.A.A.ACCURACY',/,6X,'NUMBER' 2,6X,'(+VE NORTH)',7X,'(+VE WEST)',9X,'(METERS)',4X,'(METERS)',6X,' 3(MGALS)',9X,'(MGALS)',/,16X,'DEG',2X,'MIN',2X,'SEC',6X,'DEG',2X,'M 4IN',2X,'SEC') ICOUNT=0 DO 20 J=1,NBM ICOUNT=ICOUNT+1 IF(ICOUNT.EQ.IPAGE) WRITE(6,73) IF(ICOUNT.EQ.IPAGE) WRITE(6,87) (TITLE(I),I=1,20) IF(ICOUNT.EQ.IPAGE) WRITE(6,9) LATDI=IDEG(PHI(J),MLATI,SLATI) LOGDI=IDEG(ALONG(J),MLOGI,SLOGI) BMA=BMN(J) WRITE(6,11)BMA,LATDI,MLATI,SLATI,LOGDI,MLOGI,SLOGI,HT(J),HTACC(J), 1DGFA(J),DGACC(J) 11 FORMAT(7X,A4,5X,I3,2X,I2,2X,F5.2,5X,I3,2X,I2,2X,F5.2,6X,F7.2,6X,F5 1.2,8X,F7.2,9X,F6.2) 20 CONTINUE C C PRINT-OUT HEAHING AND SUBTITLES FOR THE TABLE OF COMPUTED GRAVITY CORRECTIONS C (HELMERT, VIGNAL, AND DYNAMIC), ALONG THE INPUTTED LEVELLING LINE (OR LOOP). WRITE(6,73) WRITE(6,87) (TITLE(I),I=1,20) WRITE(6,63) 63 FORMAT(34X,'***TABLE OF COMPUTED GRAVITY CORRECTIONS***',/20X,'FOR 1 EACH INDIVIDUAL LEVELLING-SECTION ALONG THE GIVEN LEVELLING-ROUTE 2',/20X,'========================================================== 3==========='/) C C COMPUTE AND PRINT-OUT THE TOTAL NUMBER OF LEVELLING SECTIONS FORMING THE LINE NUMSEC=NBM-1 WRITE(6,93) NUMSEC 93 FORMAT(10X,'TOTAL NUMBER OF LEVELLING-SECTIONS ALONG THE LINE =',I 14) WRITE(6,1) 1 FORMAT(/,53X,'G R A V I T Y - C O R R E C T I O N S',/,3X,'LEVELLI 1NG-SECTION',2X,'LENGTH',3X,'AZIMUTH',9X,'HELMERT',14X,'VIGNAL',15X 2,'DYNAMIC',/,6X,'FROM',4X,'TO',7X,'(KM)',3X,'(DEGREES)',2X,'MAGNIT 2UDE',2X,'STANDARD',2X,'MAGNITUDE',2X,'STANDARD',2X,'MAGNITUDE',2X, 4'STANDARD',/,44X,'(MM)',4X,'(MM/KM.)',5X,'(MM)',4X,'(MM/KM.)',5X,' 5(MM)',4X,'(MM/KM.)') C SUMDIS=0.0D 0 SUMDH=0.0D 0 SUMDV=0.0D 0 SUMDD=0.0D 0 SUMHVR=0.0D 0 SUMVVR=0.0D 0 SUMDVR=0.0D 0 ACCDIS(1)=0.0D 0 ACCSTD(1)=0.0D 0 ACCHGC(1)=0.0D 0 ACCVGC(1)=0.0D 0 ACCDGC(1)=0.0D 0 ACCHSD(1)=0.0D 0 ACCVSD(1)=0.0D 0 ACCDSD(1)=0.0D 0 ICOUNT=0 C DO 80 K=2,NBM ICOUNT=ICOUNT+1 IF(ICOUNT.EQ.ISKIP) WRITE(6,73) IF(ICOUNT.EQ.ISKIP) WRITE(6,87) (TITLE(I),I=1,20) IF(ICOUNT.EQ.ISKIP) WRITE(6,1) PHI1=PHI(K-1) PHI2=PHI(K) ALONG1=ALONG(K-1) ALONG2=ALONG(K) HT1=HT(K-1) HT2=HT(K) HTACC1=HTACC(K-1) HTACC2=HTACC(K) DGFA1=DGFA(K-1) DGFA2=DGFA(K) DGACC1=DGACC(K-1) DGACC2=DGACC(K) C C COMPUTE AND PRINT-OUT THE LENGTH (KM), AZIMUTH (DEGREES), HELMERT & C VIGNAL & DYNAMIC GRAVITY CORRECTIONS- ALONG WITH THEIR STANDARDIZED VALUES. C (I.E. NORMALIZED : IN MM./KM.).. C FOR EACH LEVELLING SECTION IN THE INPUTTED LINE (OR LOOP). C ANAD=6378206.4D 0 BNAD=6356583.8D 0 ALEST1=2.0D 0*PI-ALONG1 ALEST2=2.0D 0*PI-ALONG2 CALL VININV(ANAD,BNAD,PHI1,ALEST1,PHI2,ALEST2,AZ12,AZ21,DIST12) DIST12=DIST12/1000.0D 0 AZIM12=AZ12*RHODEG SUMDIS=SUMDIS+DIST12 ACCDIS(K)=SUMDIS ACCSTD(K)=STDEVN*DSQRT(ACCDIS(K)) C CALL DHCORG(PHI1,PHI2,HT1,HTACC1,HT2,HTACC2,DGFA1,DGACC1,DGFA2,DGA 1CC2,1,DELTAH,SIGMAH) SUMDH=SUMDH+DELTAH SUMHVR=SUMHVR+SIGMAH**2 ACCHGC(K)=SUMDH ACCHSD(K)=DSQRT(SUMHVR) C CALL DHCORG(PHI1,PHI2,HT1,HTACC1,HT2,HTACC2,DGFA1,DGACC1,DGFA2,DGA 1CC2,2,DELTAV,SIGMAV) SUMDV=SUMDV+DELTAV SUMVVR=SUMVVR+SIGMAV**2 ACCVGC(K)=SUMDV ACCVSD(K)=DSQRT(SUMVVR) C CALL DHCORG(PHI1,PHI2,HT1,HTACC1,HT2,HTACC2,DGFA1,DGACC1,DGFA2,DGA 1CC2,3,DELTAD,SIGMAD) SUMDD=SUMDD+DELTAD SUMDVR=SUMDVR+SIGMAD**2 ACCDGC(K)=SUMDD ACCDSD(K)=DSQRT(SUMDVR) C C STANDARDIZE THE COMPUTED GRAVITY-CORRECTIONS ***FOR EACH 'SINGLE' LEVELLING- C SECTION*** TO CORRESPOND TO 1 KM. LENGTH. C NOTE THAT THE NORMALIZATION PROCESS HERE MEANS DIVIDING BY THE SQUARE-ROOT C OF THE LENGTH (IN KM.) OF THE GIVEN LEVELLING LINE . THIS IS CONSISTENT C WITH THE INTERNATIONALLY ADOPTED NOTATION FOR THE FORMULA EXPRESSING THE C ACCUMULATED ERRORS IN PRECISE LEVELLING . C SQRS12=DSQRT(DIST12) HGC1KM=DELTAH/SQRS12 VGC1KM=DELTAV/SQRS12 DGC1KM=DELTAD/SQRS12 C C PRINT-OUT THE TABLE OF GRAVITY-CORRECTIONS FOR ALL THREE-SYSTEMS OF HEIGHTS BMN1=BMN(K-1) BMN2=BMN(K) WRITE(6,45) BMN1,BMN2,DIST12,AZIM12,DELTAH,HGC1KM,DELTAV,VGC1KM,DE 1LTAD,DGC1KM 45 FORMAT(6X,A4,3X,A4,3X,F6.2,5X,F7.2,2X,F9.4,2X,F8.4,2X,F9.4,2X,F8.4 1,2X,F9.4,2X,F8.4) 80 CONTINUE C C PRINT-OUT A SUMMARY OF RESULTS FOR THE ACCUMULATED GRAVITY-CORRECTIONS C (ALONG WITH THEIR CORRESPONDING ACCUMULATED STANDARD DEVIATIONS)- AT EACH C BENCH MARK ALONG THE GIVEN LEVELLING-ROUTE- FOR ALL THREE HEIGHT-SYSTEMS WRITE(6,73) WRITE(6,87) (TITLE(I),I=1,20) WRITE(6,47) 47 FORMAT(35X,'*** S U M M A R Y O F R E S U L T S ***'/,14X,'TABLE 1 OF ACCUMULATED GRAVITY-CORRECTIONS AND THEIR ACCUMULATED STANDARD 2 DEVIATIONS',/,27X,'***AT EACH BENCH MARK ALONG THE GIVEN LEVELLIN 3G ROUTE***',/,30X,'=============================================== 4==='/) WRITE(6,49) 49 FORMAT(41X,'*A C C U M U L A T E D* G R A V I T Y - C O R R E C T 1I O N S',/,3X,'BENCH-MARK',2X,'ACCUMULATED',2X,'ACCUMULATED',08X,' 2HELMERT',14X,'VIGNAL',15X,'DYNAMIC',/,5X,'NUMBER',7X,'LENGTH',5X,' 4ST. ERROR',3X,'MAGNITUDE',2X,'ST. DEV.',2X,'MAGNITUDE',2X,'ST. DEV 5.',2X,'MAGNITUDE',2X,'ST. DEV.',/,19X,'(KM)',8X,'(MM)',09X,'(MM)', 66X,'(MM)',7X,'(MM)',6X,'(MM)',7X,'(MM)',6X,'(MM)') ICOUNT=0 DO 50 J=1,NBM ICOUNT=ICOUNT+1 IF(ICOUNT.EQ.IPAGE) WRITE(6,73) IF(ICOUNT.EQ.IPAGE) WRITE(6,87) (TITLE(I),I=1,20) IF(ICOUNT.EQ.IPAGE) WRITE(6,49) WRITE(6,53) BMN(J),ACCDIS(J),ACCSTD(J),ACCHGC(J),ACCHSD(J),ACCVGC( 1J),ACCVSD(J),ACCDGC(J),ACCDSD(J) 53 FORMAT(6X,A4,8X,F6.2,4X,F9.4,3X,F9.4,2X,F8.4,2X,F9.4,2X,F8.4,2X,F9 1.4,2X,F8.4) 50 CONTINUE GO TO 555 777 CONTINUE WRITE(6,73) 73 FORMAT('1',////) 99 STOP END C C*********************************************************************** C C*************************************************************************** C* * C* FUNCTION 'RADIAN' : * C* ----------------- * C* THE PURPOSE OF THIS FUNCTION SUBROUTINE IS TO * C* CONVERT THE GIVEN ANGLE IN (DEG,MIN & SEC) INTO ITS EQUIVALENT VALUE * C* IN RADIANS TO BE USED IN OTHER COMPUTATION. * C* * C* THIS FUNCTION SUBROUTINE HAS AN ENTRY-POINT 'IDEG' TO PERFORM THE * C* INVERSE OF THE ABOVE CASE , NAMELY TO CONVERT THE GIVEN ANGLE FROM * C* RADIANS TO ITS EQUIVALENT VALUE IN (DEG,MIN & SEC). * C* * C* NOTICE : IF THE INPUTTED ANGLE IN RADIANS IS NEGATIVE ,THE -VE SIGN * C* WILL APPEAR ONLY BESIDE THE DEGREES AFTER CONVERSION ,AND THIS WILL * C* IMPLY THAT THE WHOLE ANGLE (DEG,MIN & SEC) IS -VE. * C* * C* THIS FUNCTION SUBROUTINE AND ITS ENTRY-POINT ARE VALID FOR ONLY ONE * C* CONVERSION AT A TIME. * C* * C*************************************************************************** FUNCTION RADIAN(DEG,XMIN,SEC) IMPLICIT REAL*8(A-H,O-Z) PI=3.141592653589793 RADIAN=(DEG+XMIN/60.+SEC/3600.)*PI/180. RETURN C C*************************************************************************** C ENTRY IDEG(RAD,MIN,SEC) PI=3.141592653589793 SIGN=RAD RAD=DABS(RAD) DEG=RAD*180./PI+0.005/3600. IDEG=DEG MIN=(DEG-IDEG)*60. SEC=((DEG-IDEG)*60.-MIN)*60.-0.005 IF(SIGN.LT.0.) IDEG=-IDEG AMIN=MIN AMIN=DABS(AMIN) MIN=AMIN SEC=DABS(SEC) RAD=SIGN RETURN END SUBROUTINE VININV(A,B1,PHI1,ONG1,PHI2,ONG2,A12,A21,S12) 00000010 C 00000020 C SUBROUTINE 'VININV' INVERSE SOLUTION OF GEODESICS ON THE ELLIPSOID 00000030 C BETWEEN POINTS 1 AND 2 00000040 C 00000050 C REFERENCE: 'DIRECT AND INVERSE SOLUTIONS OF GEODESICS ON THE 00000060 C ELLIPSOID WITH APPLICATION OF NESTED EQUATIONS' BY T. VINCENTY 00000070 C SURVEY REVIEW # 176 APRIL 1975 PAGES 88 TO 93 00000080 C 00000090 C INPUT: A=SEMI-MAJOR AXIS OF REFERENCE ELLIPSOID 00000100 C B1=SEMI-MINOR AXIS OF REFERENCE ELLIPSOID *OR* RECIPROCAL O00000110 C THE FLATTENING 00000120 C PHI1=GEODETIC LATITUDE OF POINT 1 00000130 C ONG1=GEODETIC LONGITUDE OF POINT 1 00000140 C PHI2=GEODETIC LATITUDE OF POINT 2 00000150 C ONG2=GEODETIC LONGITUDE OF POINT 2 00000160 C OUTPUT: A12=GEODESIC AZIMUTH POINT 1 TO 2 00000170 C A21=GEODESIC AZIMUTH POINT 2 TO 1 00000180 C S12=GEODESIC DISTANCE POINT 1 TO 2 00000190 C 00000200 C NOTES: ALL LENGTH MEASURES ARE IN METRES 00000210 C ALL ANGULAR QUANTITIES ARE IN RADIANS 00000220 C LATITUDES POSITIVE NORTH 00000230 C LONGITUDES POSITIVE EAST 00000240 C AZIMUTHS POSITIVE CLOCKWISE FROM NORTH 00000250 C ITERATION TO 1*10-12 RADIANS 00000260 C ACCURACY OF DISTANCES GREATER THAN 0.01 MM AT ANY DISTANCE 00000270 C ACCURACY OF AZIMUTHS GREATER THAN 1*10-5 ARC-SECONDS 00000280 C ROUTINE CANNOT HANDLE ANTIPODAL SOLUTION - MESSAGE PRINTED 00000290 C AND ROUTINE STOPS IF ANTIPODAL SOLUTION FOUND 00000300 C FOR EQUIVALENT DIRECT SOLUTION SEE SUBROUTINE 'VINDIR' 00000310 C 00000320 C C. CHAMBERLAIN MARCH 1976 00000330 C 00000340 IMPLICIT REAL*8(A-H,O-Z) 00000350 DATA PI/3.141592653589793D0/,SING/1.0D-12/ 00000360 C 00000370 C ELLIPSOID PARAMETERS (ALSO VALID FOR A SPHERE OF RADIUS A) 00000380 C 00000390 IF(B1.EQ.0.0)B1=A 00000400 IF(B1.NE.0.0.AND.B1.LT.1000.0)F=1.0D0/B1 00000410 IF(B1.NE.0.0.AND.B1.LT.1000.0)B=A-A*F 00000420 IF(B1.GT.1000.0)B=B1 00000430 IF(B1.GT.1000.0)F=(A-B)/A 00000440 IF(B.EQ.A)B1=0.0D0 00000450 B2=B*B 00000460 OMF=1.0D0-F 00000470 C 00000480 C REDUCED LATITUDE OF POINTS 1 AND 2 AND THEIR TRIG FUNCTIONS AND 00000490 C COMBINATIONS 00000500 C 00000510 TPI=2.0D0*PI 00000520 POT=PI/2.0D0 00000530 U1=POT 00000540 IF(PHI1.LT.0.0)U1=-U1 00000550 IF(DABS(DABS(PHI1)-POT).GT.SING)U1=DATAN(OMF*DTAN(PHI1)) 00000560 SU1=DSIN(U1) 00000570 CU1=DCOS(U1) 00000580 U2=POT 00000590 IF(PHI2.LT.0.0)U2=-U2 00000600 IF(DABS(DABS(PHI2)-POT).GT.SING)U2=DATAN(OMF*DTAN(PHI2)) 00000610 SU2=DSIN(U2) 00000620 CU2=DCOS(U2) 00000630 CU1SU2=CU1*SU2 00000640 SU1CU2=SU1*CU2 00000650 SU1SU2=SU1*SU2 00000660 CU1CU2=CU1*CU2 00000670 C 00000680 C CHANGE IN LONGITUDE ON AUXILLARY SPHERE ITERATION TO 1*10-12 RADIANS00000690 C (EQUATIONS 13,14,15,16,17,18,10,11) 00000700 C 00000710 ITER=0 00000720 DLON=ONG2-ONG1 00000730 IF(DABS(DLON-PI).LT.SING.AND.DABS(PHI1).LT.SING.AND.DABS(PHI2).LT.00000740 1SING)GO TO 30 00000750 DLAS=DLON 00000760 10 SDLAS=DSIN(DLAS) 00000770 CDLAS=DCOS(DLAS) 00000780 SSIG=DSQRT((CU2*SDLAS)**2+(CU1SU2-SU1CU2*CDLAS)**2) 00000790 CSIG=SU1SU2+CU1CU2*CDLAS 00000800 TSIG=SSIG/CSIG 00000810 SIG=DARCOS(CSIG) 00000820 SALPHA=CU1CU2*SDLAS/SSIG 00000830 CALPH2=1.0D0-SALPHA**2 00000840 CTSM=0.0D0 00000850 IF(DABS(CALPH2).GT.SING )CTSM=CSIG-2.0D0*SU1SU2/CALPH2 00000860 CTSM2=CTSM**2 00000870 C=F/16.0D0*CALPH2*(4.0D0+F*(4.0D0-3.0D0*CALPH2)) 00000880 DLASUP=DLON+(1.0D0-C)*F*SALPHA*(SIG+C*SSIG*(CTSM+C*CSIG*00000890 1 (-1.0D0+2.0D0*CTSM2))) 00000900 IF(DABS(DLASUP-DLAS).LT.SING )GO TO 20 00000910 DLAS=DLASUP 00000920 ITER=ITER+1 00000930 IF(ITER.GT.50)GO TO 30 00000940 GO TO 10 00000950 C 00000960 C GEODESIC DISTANCE (EQUATIONS 3,4,6,19) 00000970 C 00000980 20 USQ=CALPH2*(A*A-B2)/B2 00000990 A3=1.0D0+USQ/16384.0D0*(4096.0D0+USQ*(-768.0D0+USQ*(320.0D0- 00001000 1 175.0D0*USQ))) 00001010 B4=USQ/1024.0D0*(256.0D0+USQ*(-128.0D0+USQ*(74.0D0-47.0D0*USQ)))00001020 DELSIG=B4*SSIG*(CTSM+B4/4.0D0*(CSIG*(-1.0D0+2.0D0*CTSM2)-B4/ 00001030 1 6.0D0*CTSM*(-3.0D0+4.0D0*SSIG**2)*(-3.0D0+4.0D0*CTSM2))) 00001040 S12=B*A3*(SIG-DELSIG) 00001050 C 00001060 C GEODESIC AZIMUTHS 1 TO 2 AND 2 TO 1 (EQUATIONS 20 AND 21) 00001070 C 00001080 A12=DATAN2(CU2*SDLAS,CU1SU2-SU1CU2*CDLAS) 00001090 A21=DATAN2(-CU1*SDLAS,SU1CU2-CU1SU2*CDLAS) 00001100 IF(A12.GT.TPI)A12=A12-TPI 00001110 IF(A12.LT.0.0)A12=A12+TPI 00001120 IF(A21.GT.TPI)A21=A21-TPI 00001130 IF(A21.LT.0.0)A21=A21+TPI 00001140 RETURN 00001150 C 00001160 C ERROR EXIT FOR ANTIPODAL SOLUTION 00001170 C 00001180 30 PRINT 1000 00001190 STOP 00001200 1000 FORMAT('1'//////////15X,'SUBROUTINE VININV WILL NOT HANDLE ANTIPOD00001210 1AL SOLUTION'/'1') 00001220 END 00001230 SUBROUTINE DHCORG(PHI1,PHI2,H1,SIGH1,H2,SIGH2,DGF1,SIGDG1,DGF2,SIG 1DG2,ISYSTM,DELTA,SIGMA) IMPLICIT REAL*8(A-H,O-Z) DOUBLE PRECISION DSIN,DCOS,DSQRT C C COMPUTATION OF THE AVERAGE HEIGHT, HEIGHT DIFFERENCE, AVERAGE FREE-AIR C ANOMALY, AND THE FREE-AIR ANOMALY DIFFERENCE BETWEEN THE GIVEN TWO C POINTS 1 & 2 C H12BAR=(H1+H2)/2.D 0 DH12=H2-H1 DGFBAR=(DGF1+DGF2)/2.D 0 DDGF12=DGF1-DGF2 C C COMPUTATION OF THE NORMAL GRAVITY AT POINTS 1 & 2 ,BASED ON THE C INTERNATIONAL FORMULA OF 1967. C GAMAEQ=978031.8D 00 C0=1.0D 00 C1=0.0053024D 00 C2=0.0000059D 00 T11=(DSIN(PHI1))**2 T21=(DSIN(2.*PHI1))**2 T12=(DSIN(PHI2))**2 T22=(DSIN(2.*PHI2))**2 GAMA1=GAMAEQ*(C0+C1*T11-C2*T21) GAMA2=GAMAEQ*(C0+C1*T12-C2*T22) C C COMPUTATION OF THE NORMAL GRAVITY AT POINTS 1 & 2 ,BASED ON THE C USC&GS USED FORMULA (THEIR SPECIAL PUBLICATION # 18) ,IN 1914. C GREF45=980624.0D 00 C0STAR=1.0D 00 C1STAR=0.002644D 00 C2STAR=0.000007D 00 T1STAR=DCOS(2.*PHI1) T2STAR=DCOS(2.*PHI2) G1STAR=GREF45*(C0STAR-C1STAR*T1STAR+C2STAR*(T1STAR**2)) G2STAR=GREF45*(C0STAR-C1STAR*T2STAR+C2STAR*(T2STAR**2)) C C COMPUTATION OF THE DIFFERENCES BETWEEN THE TWO ABOVE FORMULAS OF C NORMAL GRAVITY (INTERNATIONAL 1930 - USC&GS 1914) AT POINTS 1 & 2 C DGAMA1=GAMA1-G1STAR DGAMA2=GAMA2-G2STAR DDGAMA=DGAMA1-DGAMA2 DDGBAR=(DGAMA1+DGAMA2)/2.D 0 C C ESTABLISHING THE UPPER-TRIANGLE ELEMENTS OF THE VARIANCE-COVARIANCE C MATRIX OF THE: AVERAGE HEIGHT BETWEEN POINTS 1&2 ; HEIGHT DIFFERENCE C BETWEEN THEM ; FREE-AIR ANOMALY AT POINT 1 ; AND FREE-AIR ANOMALY AT C POINT 2. THIS V-C MATRIX WILL BE COMMON TO ALL THE THREE HEIGHT C SYSTEMS IN QUESTION. (IT SHOULD BE NOTED HERE THAT IN THIS SUBROUTINE C WE ARE NOT GOING TO DEAL WITH MATRICES AT ALL, BUT RATHER WE WORK WITH C THEIR APPROPRIATE NEEDED ELEMENTS, TO SAVE STORAGE. C S11=0.25D 0*(SIGH1**2+SIGH2**2) S12=0.50D 0*(SIGH2**2-SIGH1**2) S13=0.1543D 0*(SIGH1**2) S14=0.1543D 0*(SIGH2**2) S22=4.D 0*S11 S23=-2.D 0*S13 S24=2.D 0*S14 S33=-0.3086D 0*S23+0.0025D 0 IF(SIGDG1.NE.0.0D 0) S33=SIGDG1**2 S34=0.0 S44=0.3086D 0*S24+0.0025D 0 IF(SIGDG2.NE.0.0D 0) S44=SIGDG2**2 C C COMPUTATION OF THE HELMERT'S ORTHOMETRIC DIFFERENCE - DELTA ORTHOM. C EQUALS TO (RIGOROUS - NONRIGOROUS) ORTHOMETRIC CORRECTION. - IN METERS. C IF(ISYSTM.NE.1) GO TO 10 DELTAO=H12BAR*(DDGF12+DDGAMA+0.2238D 00*DH12)/GREF45 DELTA=DELTAO C C COMPUTATION OF THE LINEAR COEFFICIENT-MATRIX OF THE HELMERT'S- C ORTHOMETRIC-DIFFERENCE FUNCTION, WHICH IS NECESSARY FOR PROPAGATING C THE ERRORS AND DEDUCING ITS ACCURACY. C Q11=DELTA/H12BAR Q12=0.2238D 0*H12BAR/GREF45 Q13=H12BAR/GREF45 Q14=-Q13 GO TO 100 C C COMPUTATION OF THE VIGNIAL'S ORTHOMETRIC DIFFERENCE - DELTA VIGNIAL C EQUALS TO (RIGOROUS- NONRIGOROUS) VIGNAL ORTHOMETRIC-CORRECTION. C -IN METERS. C 10 IF(ISYSTM.NE.2) GO TO 20 DELTAV=(H12BAR*DDGAMA+DH12*DGFBAR)/GREF45 DELTA=DELTAV C C COMPUTATION OF THE LINEAR-COEFFICIENT MATRIX OF THE VIGNAL'S- C ORTHOMETRIC-DIFFERENCE FUNCTION, WHICH IS NEEDED FOR PROPAGATING C THE ERRORS AND DEDUCING ITS ACCURACY. C Q11=0.0 Q12=DGFBAR/GREF45 Q13=DH12/(2.*GREF45) Q14=Q13 GO TO 100 C C COMPUTATION OF THE DYNAMIC DIFFERENCE- DELTA DYNAMIC C EQUALS TO (RIGOROUS - NONRIGOROUS) DYNAMIC CORRECTION C -IN METERS. C 20 IF(ISYSTM.NE.3) GO TO 30 DELTAD=(DGFBAR+DDGBAR)*DH12/GREF45 DELTA=DELTAD C C COMPUTATION OF THE LINEAR-COEFFICIENT-MATRIX OF THE DYNAMIC--DIFFERENCE C FUNCTION, WHICH IS NEEDED FOR PROPAGATING THE ERRORS AND DEDUCING ITS C ACCURACY. C A0=-6.295D 0 A1=-0.358D 0 A2=1.076D 0 PHIBAR=(PHI1+PHI2)/2. Q11=0.0D 0 Q12=(DGFBAR+A0-A1*DSIN(PHIBAR)**2+A2*DSIN(2.*PHIBAR)**2)/GREF45 Q13=DH12/(2.*GREF45) Q14=Q13 GO TO 100 30 WRITE(6,50) 50 FORMAT('1'////10X,'THE INPUT CODE (ISYSTM), WHICH SPECIFIES THE 1HEIGHT SYSTEM IN USE ,IS NOT CORRECT'//) GO TO 200 100 CONTINUE C C COMPUTATION OF THE ACCURACY-OR STANDARD DEVIATION-(SIGMA, IN METERS) C OF THE COMPUTED: HELMERT'S-ORTHOMETRIC-DIFFERENCE,/OR, VIGNAL'S C ORTHOMETRIC-DIFFERENCE,/, OR THE DYNAMIC-DIFFERENCE --- ACCORDING TO C THE APPROPRIATE VALUE OF THE INPUT CODE (ISYSTM), INDICATING WHICH C SYSTEM OF HEIGHTS IN USE. C R11=Q11*S11+Q12*S12+Q13*S13+Q14*S14 R12=Q11*S12+Q12*S22+Q13*S23+Q14*S24 R13=Q11*S13+Q12*S23+Q13*S33+Q14*S34 R14=Q11*S14+Q12*S24+Q13*S34+Q14*S44 SIGSQ=R11*Q11+R12*Q12+R13*Q13+R14*Q14 SIGSQ=DABS(SIGSQ) SIGMA=DSQRT(SIGSQ) C C COMPUTE THE GRAVITY-CORRECTION ALONG WITH ITS STANDARD DEVIATION - IN MM. C DELTA=DELTA*1000.0D 0 SIGMA=SIGMA*1000.0D 0 200 RETURN END