C C C*************************************************************************** C* * C* PROGRAM 'H V A D G C' : ( BY : M.M. NASSAR, OCT.,1974 ) * C* --------------------- * C* * C* THE MAIN PURPOSE OF THIS PROGRAM IS TO COMPUTE THE* C* ACTUAL GRAVITY-CORRECTION TO THE HEIGHT-DIFFERENCE (HELMERT, VIGNAL OR * C* DYNAMIC) BASED ON NORMAL GRAVITY ONLY AND BEING USED IN CANADA. * C* * C* FOR DETAILS CONCERNING THE MATHEMATICAL MODELS AND RELATED DEVELOPM- * C* ENTS, SEE : NASSAR AND VANICEK (1975).'LEVELLING AND GRAVITY'-TECHN- * C* ICAL REPORT # 33, SURVEYING ENGINEERING DEPT., U.N.B., FREDERICTON. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* THE REQUIRED INFORMATION (WHICH HAVE TO BE READ FROM DATA CARDS) AT THE * C* TWO END-POINTS OF THE LEVELLING SECTION IN QUESTION ARE AS FOLLOWS-FOR * C* EACH POINT: LATITUDE AND LONGITUDE (IN DECIMAL DEGREES), ELEVATION * C* ABOVE SEA-LEVEL AND ITS ASSIGNED STANDARD DEVIATION (IN FEET), FREE-AIR * C* GRAVITY ANOMALY AND ITS STANDARD DEVIATION (IN MILLIGALS). IT SHOULD * C* BE NOTED HERE THAT IF THE STANDARD DEVIATION OF THE FREE-AIR ANOMALY * C* IS NOT AVAILABLE, THEN IT SHOULD BE INPUT ON THE DATA CARDS AS '0.0' * C* (ZERO) ,WHICH MAKES THE PROGRAM (OR ITS SUBROUTINES) COMPUTE ITS VALUE * C* INTERNALLY VIA THE RIGOROUS ERROR PROPAGATIONS. IT HAS TO BE UNDERSTOOD* C* FROM THE ABOVE THAT, THREE HEIGHT-SYSTEMS (HELMERT, VIGNAL AND DYNAMIC) * C* ARE CONSIDERED UNDER INVESTIGATION BY THIS PROGRAM. * C* * C* THE ABOVE INFORMATION HAS TO BE INPUT ON ONE DATA CARD ONLY, PER EACH * C* BENCH MARK, ACCORDING TO THE FOLLOWING FORMAT :- * C* COLS 1- 2 : BLANK * C* COLS 3-10 : LATITUDE (+VE N.), IN DECIMAL DEGREES,...............F8.5 * C* COL 11 : BLANK * C* COLS 12-20 : LONGITUDE (+VE W.), IN DECIMAL DEGREES,..............F9.5 * C* COLS 21-23 : BLANK * C* COLS 24-30 : ELEVATION (ABOVE M.S.L.), IN FEET, ..................F7.1 * C* COLS 31-35 : BLANK * C* COLS 36-40 : ELEVATION STANDARD DEVIATION, IN FEET, ..............F5.1 * C* COLS 41-43 : BLANK * C* COLS 44-50 : FREE-AIR ANOMALY, IN MILLIGALS, .....................F7.2 * C* COLS 51-54 : BLANK * C* COLS 55-60 : ANOMALY STANDARD DEVIATION, IN MILLIGALS, ...........F6.2 * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* IN ADDITION TO THE MAIN PURPOSE MENTIONED ABOVE, THE PROGRAM COMPUTES * C* AND PRINTS-OUT SOME OTHER AUXILIARY INFORMATION ABOUT THE LEVELLING * C* SECTION UNDER INVESTIGATION. THE FOLLOWING IS A SUMMARY OF ITS OUTPUT: * C* * C* 1-TABLE SUMMARISING THE INPUT DATA. * C* 2-THE COMPUTED ELLIPSOIDAL LENGTH (IN KM.) AND GEODETIC AZIMUTH (IN DEG,* C* MIN & SEC) OF THE GIVEN LEVELLING SECTION. * C* 3-THE COMPUTED HELMERT-GRAVITY-CORRECTION (ALONG WITH ITS STANDARD * C* DEVIATION)-IN MILLIMETERS, FOR THE WHOLE GIVEN LEVELLING SECTION, THEN* C* THE NORMALIZED HELMERT-GRAVITY-CORRECTION-IN MILLIMETERS-FOR ONE * C* KILOMETER DISTANCE ONLY (ALONG WITH ITS CORRESPONDING STANDARD DEV.) * C* 4-SAME AS (3) ABOVE BUT FOR THE VIGNAL HEIGHT-SYSTEM. * C* 5-SAME AS (3) ABOVE BUT FOR THE DYNAMIC HEIGHT-SYSTEM. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* NOTE THAT THIS PROGRAM (AND HIS SUBROUTINES) IS WORKING IN DOUBLE PRECI-* C* SION, AND THE NUMBER OF THE LEVELLING SECTIONS THAT CAN BE HANDLED BY * C* IT ARE UNLIMITTED. FURTHER, THE READ STATEMENT OF THE INPUT DATA CAN * C* BE EASILY MODIFIED SUCH THAT CONSECUTIVE LEVELLING SECTIONS THAT * C* CONSTITUTE AN ENTIRE LEVELLING LINE (OR LOOP) CAN BE HANDLED IN ONE * C* RUN OF THE PROGRAM. * C* * C*************************************************************************** C IMPLICIT REAL*8(A-H,O-Z) DIMENSION PHI(2),ALONG(2),HT(2),HTACC(2),DGFA(2),DGACC(2) DOUBLE PRECISION DSIN,DCOS,DABS,DSQRT PI=3.141592653589793D 0 RHO=206264.8062471D 0 RHODEG=180.D 0/PI FEETMS=0.3048D 0 C 999 CONTINUE DO 10 I=1,2 READ(5,5,END=99)PHI(I),ALONG(I),HT(I),HTACC(I),DGFA(I),DGACC(I) 5 FORMAT(2X,F8.5,1X,F9.5,3X,F7.1,5X,F5.1,3X,F7.2,4X,F6.2) PHI(I)=PHI(I)/RHODEG ALONG(I)=ALONG(I)/RHODEG HT(I)=HT(I)*FEETMS HTACC(I)=HTACC(I)*FEETMS 10 CONTINUE C WRITE(6,7) 7 FORMAT('1',//,30X,'TABLE OF GIVEN INFORMATION AT BENCH MARKS 1 AND 1 2'/,30X,'================================================='/) WRITE(6,9) 9 FORMAT(5X,'STATION',7X,'LATITUDE',10X,'LONGITUDE',11X,'HEIGHT',2X, 1'HT.ACCURACY',2X,'FREE-AIR ANOM',2X,'FA-AN ACCURACY',/,6X,'NUMBER' 2,6X,'(+VE NORTH)',7X,'(+VE WEST)',9X,'(METERS)',4X,'(METERS)',7X,' 3(MGALS)',8X,'(MGALS)',/,16X,'DEG',2X,'MIN',2X,'SEC',6X,'DEG',2X,'M 4IN',2X,'SEC') C DO 20 J=1,2 LATDI=IDEG(PHI(J),MLATI,SLATI) LOGDI=IDEG(ALONG(J),MLOGI,SLOGI) WRITE(6,11) J,LATDI,MLATI,SLATI,LOGDI,MLOGI,SLOGI,HT(J),HTACC(J), 1DGFA(J),DGACC(J) 11 FORMAT(5X,I4,7X,I3,2X,I2,2X,F7.4,3X,I3,2X,I2,2X,F7.4,5X,F7.2,5X,F5 1.2,8X,F7.2,9X,F6.2) 20 CONTINUE C CALL INDPRO(PHI(1),ALONG(1),PHI(2),ALONG(2),AZ12,AZ21,DIST) DIST=DIST/1000.D 0 DISTN=DSQRT(DIST) IDAZ=IDEG(AZ12,MAZ,SAZ) WRITE(6,13) DIST 13 FORMAT(//,5X,'THE CALCULATED DISTANCE BETWEEN POINTS 1 AND 2 =',F1 12.6,3X,'KILOMETERS') WRITE(6,15) IDAZ,MAZ,SAZ 15 FORMAT(/42X,'DEG',2X,'MIN',2X,'SEC',/,5X,'THE CALCULATED AZIMUTH F 1ROM 1 TO 2 =',I4,2X,I2,2X,F6.3) C DO 40 JJJ=1,3 ISYSTM=JJJ C CALL DHCORG(PHI(1),PHI(2),HT(1),HTACC(1),HT(2),HTACC(2),DGFA(1),DG 1ACC(1),DGFA(2),DGACC(2),ISYSTM,DELTA,SIGMA) C C NORMALIZATION OF THE COMPUTED GRAVITY-CORRECTIONS (DELTA), AND ITS ST. C DEVIATION (SIGMA)-THAT ARE EVALUATED IN(MM. ) FOR THE ENTIRE SECTION C BETWEEN POINTS 1 AND 2- TO CORRESPOND TO ONLY 1 KM IN THE DIRECTION C (AZIMUTH) 1 TO 2, AND TO BE GIVEN IN MILLIMETER-UNITS. C NOTE THAT THE NORMALIZATION PROCESS HERE MEANS DIVIDING BY THE SQUARE-ROOT C OF THE LENGTH (IN KM.) OF THE GIVEN LEVELLING SECTION. THIS IS CONSISTENT C WITH THE INTERNATIONALLY ADOPTED NOTATION FOR THE FORMULA EXPRESSING THE C ACCUMULATED ERRORS IN PRECISE LEVELLING . C DELTAN=DELTA/DISTN SIGMAN=SIGMA/DISTN WRITE(6,17) DELTAN,SIGMAN 17 FORMAT(5X,'THE CORRESPONDING *NORMALIZED*-GRAVITY-CORRECTION FOR 1 1 KM IN THE DIRECTION 1-2 =',F10.4,2X,' MM / KM ',/,47X,'WITH COR 2RESPONDING STANDARD DEVIATION =',F10.4,2X,' MM / KM '/) 40 CONTINUE C WRITE(6,50) 50 FORMAT(//,8X,'NOTE :- **NORMALIZED** IN THE ABOVE TABLE MEANS DIVI 1DED BY THE SQUARE-ROOT'/,16X,'OF THE DISTANCE (IN KM.) BETWEEN THE 2 GIVEN POINTS 1 & 2') GO TO 999 99 STOP END C C*********************************************************************** C* * C* SUBROUTINE 'DHCORG' : (BY M.M. NASSAR, SEPT.,1974) * C* ------------------- * C* THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE * C* AND PRINT-OUT THE 'GRAVITY CORRECTION' (DUE TO THE NEGLECT OF THE * C* EARTH'S ACTUAL GRAVITY FIELD AND USING THE THEORITICAL NORMAL * C* GRAVITY INSTEAD-IN LEVELLING) TO BE ADDED TO THE HEIGHT DIFFERENCE * C* (H2-H1) BETWEEN BENCH-MARKS 1&2, AS BASED ON NORMAL GRAVITY ONLY * C* AND BEING USED BY THE NGS AND THE GSC. THE ESTIMATED STANDARD * C* DEVIATION OF THE COMPUTED GRAVITY CORRECTION IS ALSO COMPUTED AND * C* PRINTED-OUT BY THIS SUBROUTINE. * C* * C* THE FORMULAE CONTAINED HEREIN ARE BASED ON THE 1930 INTERNATIONAL * C* FORMULA FOR NORMAL GRAVITY, AND THE FORMULA FOR NORMAL GRAVITY USED * C* BY THE NGS (AND SUBSEQUENTLY BY THE GSC)-FOR DETAILS CONCERNING THE * C* LATTER FORMULA SEE THE USC&GS SPECIAL PUBLICATION NO. 18 OR 240. * C* * C* FOR DETAILS CONCERNING THE MATHEMATICAL MODELS AND RELATED DEVELOPM-* C* ENTS, SEE : NASSAR AND VANICEK (1975).'LEVELLING AND GRAVITY'-TECHN-* C* ICAL REPORT # 33, SURVEYING ENGINEERING DEPT., U.N.B., FREDERICTON. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* NOTE THAT SUCH GRAVITY CORRECTION TO (H2-H1) WILL BE REFERRED TO * C* (IN THIS SUBROUTINE) AS DIFFERENCE (HELMERT'S, VIGNAL'S, DYNAMIC), * C* SINCE IT IS ACTUALLY THE DIFFERENCE BETWEEN THE RIGOROUS CORRECTION * C* (HELMERT, VIGNAL, DYNAMIC) BASED ON THE ACTUAL GRAVITY, AND THE * C* CORRESPONDING APPROXIMATE CORRECTION BASED ON NORMAL GRAVITY. * C* * C* THIS SUBROUTINE COMPUTES THE ABOVE MENTIONED DIFFERENCE (ALONG WITH * C* ITS CORRESPONDING ST. DEV.) BY CONSIDERING EITHER ONE OF THE FOLLOW-* C* ING THREE HEIGHT-SYSTEMS: HELMERT'S (RIGOROUS)-ORTHOMETRIC HEIGHTS; * C* VIGNAL'S (APPROXIMATE)-ORTHOMETRIC HEIGHTS; OR DYNAMIC HEIGHTS .... * C* DEPENDING ON THE APPROPRIATE INPUTTED VALUE FOR THE CODE (ISYSTM). * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* USAGE:- CALL DHCORG(PHI1,PHI2,H1,SIGH1,H2,SIGH2,DGF1,SIGDG1, * C* DGF2,SIGDG2,ISYSTM,DELTA,SIGMA) * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-***** C* * C* THE REQUIRED INPUT DATA ARE :- (THROUGH THE CALLING ROUTINE) * C* ----------------------------- * C* PHI1 ,PHI2 : LATITUDES (+VE NORTH) OF POINTS 1&2 RESPECT. .IN RADIANS * C* H1 ,H2 : ELEVATIONS (ABOVE MSL) OF POINTS 1&2 RESPECT...IN METERS * C* SIGH1 ,SIGH2 : STANDARD DEVIATIONS ASSIGNED TO H1&H2 RESPECT..IN METERS * C* DGF1 ,DGF2 : FREE-AIR ANOMALIES AT POINTS 1&2 RESPECT. (EITHER DEDUCED* C* FROM THE OBSERVED GRAVITIES-IF ANY-AT THESE POINTS, OR * C* COMPUTED VIA AN INTERPOLATION TECHNIQUE).......IN MGALS * C* SIGDG1,SIGDG2: STANDARD DEVIATIONS OF DGF1&DGF2 RESPECT., AS OBTAINED * C* FROM THE INTERPOLATION TECHNIQUE-USED TO OBTAIN THE FREE-* C* AIR ANOMALIES .................................IN MGALS * C* HOWEVER, IF EITHER DGF1 OR DGF2 HAPPENS TO BE DEDUCED * C* FROM THE OBSERVED GRAVITY AT ITS CORRESPONDING BENCH MARK* C* ,THEN ITS ACCURACY WILL BE ESTIMATED FROM THE ACCURACY * C* ASSIGNED TO THE HEIGHT OF THE SAME BENCH MARK. THIS WILL * C* BE AUTOMATICALLY DONE IN THIS SUBROUTINE ONCE WE INPUT * C* ZERO VALUE FOR THE ACCURACY OF THE FREE-AIR ANOMALY IN * C* QUESTION, I.E., SIGDG1=0. AND/OR SIGDG2=0., WHICHEVER * C* THE CASE IS APPLICABLE. * C* ISYSTM : AN INTEGER CODE INDICATING THE DESIRED HEIGHT-SYSTEM * C* TO BE USED AS FOLLOWS :- * C* ISYSTM=1 ,FOR HELMERT'S-ORTHOMETRIC HEIGHTS, * C* ISYSTM=2 ,FOR VIGNAL'S-ORTHOMETRIC HEIGHTS, * C* ISYSTM=3 ,FOR DYNAMIC HEIGHTS. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-***** C* * C* THE OUTPUT RESULTS WILL BE :- * C* -------------------------- * C* DELTA : THE COMPUTED GRAVITY-CORRECTION (HELMERT, VIGNAL OR DYNAMIC) * C* DEPENDING ON THE SPECIFIED VALUE FOR THE CODE (ISYSTM)-IN MM. * C* SIGMA : THE ESTIMATED ACCURACY (STANDARD DEVIATION)) OF THE CORR- * C* ESPONDING COMPUTED 'DELTA' ...........................IN MM. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-***** C* * C* IT SHOULD BE NOTED HERE THAT THIS SUBROUTINE IS WORKING IN DOUBLE * C* PRECISION, AND FOR ONE LEVELLING SECTION ONLY AT A TIME. IN OTHER * C* WORDS, IT CAN COMPUTE AND PRINT-OUT THE GRAVITY CORRECTION (ALONG * C* WITH ITS STANDARD DEVIATION) FOR THE HEIGHT DIFFERENCE BETWEEN TWO * C* POINTS ONLY AT A TIME. * C* CONSEQUENTLY, IT CAN BE CALLED WITHIN A DO-LOOP (IN THE CALLING ROUTINE * C* WITH APPROPRIATE VALUES FOR ITS ARGUMENTS, TO PERFORM THE GRAVITY- * C* CORRECTIONS FOR ALL THE LEVELLING SECTIONS CONSTITUING THE ENTIRE LINE * C* * C*************************************************************************** C 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 1930 C GAMAEQ=978049.0D 00 C0=1.0D 00 C1=0.0052884D 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=10.905D 0 A1=13.244D 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 SIGMA=DSQRT(SIGSQ) C C PRINT-OUT THE FINAL RESULTS, NAMELY: THE COMPUTED GRAVITY-CORRECTION (HELMERT, C VIGNAL'S OR DYNAMIC) ALONG WITH ITS CORRESPONDING STANDARD DEVIATION. C -BOTH ARE PRINTED IN MILLIMETERS. C DELTA=DELTA*1000.0D 0 SIGMA=SIGMA*1000.0D 0 IF(ISYSTM.EQ.1) WRITE(6,60)DELTA,SIGMA IF(ISYSTM.EQ.2) WRITE(6,70)DELTA,SIGMA IF(ISYSTM.EQ.3) WRITE(6,80)DELTA,SIGMA 60 FORMAT(//,4X,'THE HELMERT-GRAVITY-CORRECTION FOR THE WHOLE SECTION 1 BETWEEN THE GIVEN POINTS 1&2 =',F10.4,2X,'MILLIMETERS'/,60X,'WITH 2 A STANDARD DEVIATION =',F10.4,2X,'MILLIMETERS'//) 70 FORMAT(//,4X,' THE VIGNAL-GRAVITY-CORRECTION FOR THE WHOLE SECTION 1 BETWEEN THE GIVEN POINTS 1&2 =',F10.4,2X,'MILLIMETERS'/,60X,'WITH 2 A STANDARD DEVIATION =',F10.4,2X,'MILLIMETERS'//) 80 FORMAT(//,4X,'THE DYNAMIC-GRAVITY-CORRECTION FOR THE WHOLE SECTION 1 BETWEEN THE GIVEN POINTS 1&2 =',F10.4,2X,'MILLIMETERS'/,60X,'WITH 2 A STANDARD DEVIATION =',F10.4,2X,'MILLIMETERS'//) 200 RETURN END 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*************************************************************************** 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 C C*************************************************************************** C* * C* SUBROUTINE 'INDPRO' : * C* ------------------- * C* THE PURPOSE OF THIS SUBROUTINE IS TO SOLVE THE * C* INVERSE GEODETIC PROBLEM (ON THE NAD REFERENCE ELLIPSOID), NAMELY : * C* GIVEN THE LATITUDES AND LONGITUDES OF TWO POINTS, THEN SOLVE FOR THE * C* TWO AZIMUTHS (FOREWARD & BACKWARD) AND THE DISTANCE BETWEEN THE TWO * C* POINTS. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* THIS SUBROUTINE WAS WRITTEN BY 'DON THOMSON' IN 1969, BASED ON THE * C* 'BESSEL'S APPROACH' FOR LONG LINES ON THE ELLIPSOID-(REFERENCE : * C* HANDBUCH DER VERMESSENGSKUNDE, BD.III-BY 'JORDON & EGGERT', 1962 - * C* ENGLISH TRANSLATION, ARMY MAP SERVICE, WASHINGTON.) * C* *- C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* USAGE : CALL INDPRO(PHI1,ONG1,PHI2,ONG2,AZ12,AZ21,DIST) * C* * C* REQUIRED INPUT DATA ARE :- (THROUGH THE CALLING ROUTINE) * C* ----------------------- * C* * C* PHI1 : LATITUDE OF FIRST POINT '1' (+VE NORTH)- IN RADIANS, * C* ONG1 : LONGITUDE OF FIRST POINT '1' (+VE WEST)- IN RADIANS, * C* PHI2 : LATITUDE OF SECOND POINT '2' (+VE NORTH)- IN RADIANS, * C* ONG2 : LONGITUDE OF SECOND POINT '2' (+VE WEST)- IN RADIANS. * C* * C* THE OUTPUT WILL BE :- * C* ------------------ * C* * C* AZ12 : THE AZIMUTH OF THE GEODESIC FROM '1' TO '2' - IN RADIANS, * C* AZ21 : THE AZIMUTH OF THE GEODESIC FROM '2' TO '1' - IN RADIANS, * C* DIST : THE GEODESIC DISTANCE BETWEEN '1' AND '2' - IN METERS. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* NOTE : THIS SUBROUTINE IS VALID (IN NORTH AMERICA) FOR SOLVING THE * C* INVERSE GEODETIC PROBLEM BETWEEN ONLY TWO POINTS AT A TIME. * C* * C*************************************************************************** C SUBROUTINE INDPRO (PHI1,ONG1,PHI2,ONG2,AZ12,AZ21,DIST) IMPLICIT REAL*8(A-H,O-Z) REAL*8 DIFF(10) 15 A=6378206.4D 0 B=6356583.8D 0 F=(A-B)/A ESQ1=2.*F-(F*F) ESQ2=ESQ1/(1.-ESQ1) RHO=206264.8062471D 0 PI=3.14159265359D 0 ONG1=2.D+0*PI-ONG1 ONG2=2.D+0*PI-ONG2 IF(PHI1.GE.1.305D+0)GO TO 200 B1=DATAN((DSQRT(1.-ESQ1))*DTAN(PHI1)) GO TO 201 200 ANG1=DATAN((DSQRT(1.-ESQ1))*DTAN(PI-(DABS(PHI1)))) B1=PI+ANG1 IF(PHI1.LT.1.D-15)B1=-PI-ANG1 201 IF(PHI2.GE.1.305D+0)GO TO 202 B2=DATAN((DSQRT(1.-ESQ1))*DTAN(PHI2)) GO TO 203 202 ANG2=DATAN((DSQRT(1.-ESQ1))*DTAN(PI-(DABS(PHI2)))) B2=PI+ANG2 IF(PHI2.LT.1.D-15)B2=-PI-ANG2 203 BL=DABS(ONG2-ONG1) IF(BL.GT.PI)BL=2.*PI-BL AL=BL I=1 DIFF(I)=0.0D 0 C 4 T=DARCOS((DSIN(B1)*DSIN(B2))+(DCOS(B1)*DCOS(B2)*DCOS(AL))) W=(DCOS(B1)*DCOS(B2)*DSIN(AL))/DSIN(T) B0=DARCOS(W) ALPHA=DARSIN(W) X=DCOS(T)-((2.*DSIN(B1)*DSIN(B2))/(DCOS(ALPHA)*DCOS(ALPHA))) Y=(2.*(X*X))-1. Z=(4.*(X*X*X))-(3.*X) A0=1.D+0-((F/4.D+0)*((1.D+0+F)*(DCOS(ALPHA)**2)))+((3.D+0/16.D+0)* 1(F**2)*(DCOS(ALPHA)**4)) A2=((F/4.D+0)*(1.D+0+F)*(DCOS(ALPHA)**2))-((F**2/4.D+0)*(DCOS(ALPH 1A)**4)) A4=(F**2/32.D+0)*(DCOS(ALPHA)**4) U=2.*DSIN(T)*DCOS(T) V=(3.*DSIN(T))-(4.*(DSIN(T)**3)) DIFF(I+1)=(RHO*F*DSIN(ALPHA))*((A0*T)+(A2*DSIN(T)*X)+(A4*U*Y)) DELTA=DIFF(I+1)-DIFF(I) AL=BL+(DIFF(I+1)/RHO) IF(DABS(DELTA).LE.1.D-06) GO TO 5 I=I+1 GO TO 4 C C AZIMUTH COMPUTATION C 5 E1=DSIN(AL)*DCOS(B2) D1=(DSIN(B2)*DCOS(B1))-(DCOS(AL)*DSIN(B1)*DCOS(B2)) E2=DSIN(AL)*DCOS(B1) D2=(DSIN(B2)*DCOS(B1)*DCOS(AL))-(DSIN(B1)*DCOS(B2)) IF(ONG2.GT.ONG1.AND.DABS(D1).LT.1.D-15)GO TO 101 IF(ONG2.LT.ONG1.AND.DABS(D1).LT.1.D-15)GO TO 102 IF(ONG1.GT.ONG2)GO TO 103 IF(D1.GT.1.D-15)AZ12=DATAN(E1/D1) IF(D1.LT.1.D-15)AZ12=PI+DATAN(E1/D1) IF(D2.GT.1.D-15)AZ21=PI+DATAN(E2/D2) IF(D2.LT.1.D-15)AZ21=2.*PI+DATAN(E2/D2) GO TO 85 103 IF(D1.LT.1.D-15)AZ12=PI-DATAN(E1/D1) IF(D1.GT.1.D-15)AZ12=2.*PI-DATAN(E1/D1) IF(D2.LT.1.D-15)AZ21=-DATAN(E2/D2) IF(D2.GT.1.D-15)AZ21=PI-DATAN(E2/D2) GO TO 85 101 AZ12=PI/2. AZ21=3.*PI/2. GO TO 85 102 AZ12=3.*PI/2. AZ21=PI/2. C C DISTANCE COMPUTATION C 85 AKSQ=ESQ2*DSIN(B0)*DSIN(B0) AX=1.+(AKSQ/4.)-((3.*AKSQ**2)/64.)+((5.*AKSQ**3)/256.)-((175.*AKSQ 1**4)/16384.) BX=-AKSQ/4.+((AKSQ**2)/16.)-((15.*AKSQ**3)/512.)+((35.*AKSQ**4)/20 148.) CX=(-AKSQ**2/128.)+((3.*AKSQ**3)/512.)-((35.*AKSQ**4)/8192.) DX=(-AKSQ**3/1536.)+((5.*AKSQ**4)/6144.) DIST=B*((AX*T)+(BX*X*DSIN(T))+(CX*Y*DSIN(2.*T))+(DX*Z*DSIN(3.*T))) ONG1=2.D+0*PI-ONG1 ONG2=2.D+0*PI-ONG2 RETURN END