C *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-***-*-* C* NOTE:- THE DOCUMENTATION OF THE ENTIRE PROGRAM PACKAGE IS DONE WITHIN * C* THE PROGRAM SEGMENTS VIA SELF-EXPLAINED COMMENT CARDS. C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-***-*-* C C*************************************************************************** C* * C* PROGRAM ' A R E A G C ' : (BY: M.M. NASSAR, 1975). * C* ------------------------- REVISED:JANUARY 1976. * C* * C* THIS COMPUTER PROGRAM PACKAGE (MAIN PROGRAM * C* AND ITS ASSOCIATED SUBPROGRAMS) HAS BEEN DESIGNED MAINLY FOR THE INVEST-* C* IGATION OF THE VARIATION OF THE ACTUAL GRAVITY CORRECTIONS WITH AZIMUTH * C* , FOR A PARTICULAR HEIGHT SYSTEM (HELMERT, VIGNAL, DYNAMIC) AT ANY * C* DESIRED LOCATION IN CANADA- BASED ON THE TECHNIQUES PROPOSED BY * C* NASSAR AND VANICEK (1975). AND USING THE NEW GRAVITY DATA FILE OF THE * C* EARTH PHYSICS BRANCH- BASED ON THE 1967 IUGG SYSTEM. * 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* FOR DETAILS CONCERNING THE REVISION TO THE NEW INTERNATIONAL SYSTEM OF * C* 1967; SEE NASSAR'S PH.D. THESIS, SURVEYING ENG., UNB, OCT. 1976. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* THE PURPOSE OF THE MAIN PROGRAM 'LEVAGRAV' IS TO READ THE LATITUDE * C* AND LONGITUDE OF THE CENTRE OF THE ONE-SQUARED-AREA UNDER INVESTIGATION* C* (USUALLY ONE-DEGREE BY ONE-DEGREE OR LESS,IN BOTH LATITUDE & LONGITUDE) * C* , SEEKS THE APPROPRIATE DATA FILE AND STORES THE RELEVANT INFORMATION * C* NEEDED FOR THE PREDICTION OF BOTH THE FREE-AIR ANOMALY AND HEIGHT FIELDS* C* ALSO, IT COMPUTES THE (X,Y) LOCAL COORDINATES OF THE SELECTED DATA * C* POINTS FOR THE SAME REASON MENTIONED ABOVE. THE MAIN PROGRAM CONTAINS * C* THE POSSIBILITY OF SCALING THE COORDINATE SYSTEM XY BY A SUITABLE SCALE * C* FACTOR ,IN ORDER TO GET A WELL-CONDITIONED MATRIX OF NORMAL EQUATIONS * C* -FOR THE LEAST-SQUARES SOLUTION PROCESS, AND ALSO TO OVERCOME THE * C* PROBLEM OF OVERFLOW ON THE IBM/SYSTEMS. * C* * C* THE LEAST-SQUARES PREDICTION OF THE ANOMALY AND HEIGHT FIELDS IS * C* PERFORMED VIA THE SUBROUTINE 'O R T H O', WHICH UTILIZES THE GRAM- * C* SCHMIDT ORTHOGONALIZATION PROCESS WHEN SOLVING THE NORMAL EQUATIONS * C* SYSTEM FOR THE APPROXIMATING POLYNOMIAL COEFFICIENTS. THIS TECHNIQUE * C* HELPS IN OVERCOMMING THE PROBLEM OF UNDERFLOW AND FALSE SINGULARITY * C* DURING THE CLASSICAL INVERSION PROCESS OF NORMAL EQUATIONS MATRIX. * C* * C* THE PREDICTION OF THE ANOMALY AND HEIGHT PROFILES ,ALONG THE SIMULATED * C* LEVELLING LINES, AND THE PATTERN OF THE GRAVITY CORRECTIONS VS. AZIMUTH * C* IN CASE OF THE ABOVE MENTIONED THREE SYSTEMS, ARE OBTAINED VIA THE * C* SUBROUTINE 'GCAFAZ'. * C* * C* FINALLY, A PLOT SHOWING THE DISTRIBUTION OF DATA POINTS WITHIN THE * C* PARTICULAR SQUARE-AREA IN QUESTION, IS FURNISHED BY THE SUBROUTINE * C* 'GCPLOT'. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* REQUIRED INPUT DATA ARE :- (BOTH EXTERNAL CARDS AND INTERNAL FILES)* C* ----------------------- * C* * C* 1-SEVEN DISC DATA FILES (ALONG WITH THEIR APPROPRIATE DD-JCL CARDS) * C* ,COVERING THE CANADIAN TERRITORY AND GEOGRAPHICALLY OVERLAPPED, HAVE * C* TO BE SUPPLIED. EACH RECORD ON THESE FILES PERTAIN TO ONLY ONE * C* GRAVITY STATION AND CONTAINS ITS LATITUDE AND LONGITUDE (BOTH IN * C* DECIMAL DEGREES) , OBSERVED GRAVITY (IN GALS), FREE-AIR ANOMALY (IN * C* MGAL), ELEVATION ABOVE MSL AND ELEVATION STANDARD DEVIATION (IN FEET) * C* THE DATA REFERENCE NUMBERS (UNIT NUMBERS) OF THESE 7 DATA SETS SHOULD * C* BE: 21, 22, 23, 24, 25, 16 AND 17 RESPECTIVELY, SINCE THOSE SPECIFIC * C* NUMBERS HAVE BEEN USED IN THE MAIN PROGRAM OF THIS PACKAGE. * C* FOR MORE DETAILS ABOUT THE EXACT GEOGRAPHICAL PARTITIONING OF THE * C* NEW GRAVITY DATA, SEE: NASSAR'S PH.D. THESIS, 1976. * C* * C* 2-ONE CARD CONTAINS A CODE 'LSPRNT' FOR THE OPTIONAL PRINTING-OUT OF THE* C* INTERMEDIATE RESULTS OF THE LEAST-SQUARES APPROXIMATION TECHNIQUE. * C* LSPRNT=0 : NO PRINT-OUT OF INTERMEDIATE RESULTS IS REQUIRED. * C* LSPRNT>0 : PRINT-OUT OF THE INTERMEDIATE RESULTS IS NEEDED. * C* ACCORDING TO THE FORMAT: (I1), I.E. INTEGERS IN COLUMN ONE ONLY. * C* * C* 3-ONE CARD CONTAINING THE FOLLOWING INFORMATION ABOUT THE THE DESIRED * C* AREA UNDER INVESTIGATION. * C* POLDEG: SPECIFIED ORDER OF THE APPROXIMATING POLYNOMIAL; * C* EXTENT: THE SIDE (EXTENT) OF THE SQUARED-AREA IN LAT.&LONG. DIRECTIONS* C* MAXDIM: ALLOWABLE NUMBER OF DATA POINTS WITHIN THE AREA, FOR STORAGE. * C* DIST : INTEGRATION LENGTH ALONG SIMULATED LEVELLING LINES. (IN KM.) * C* ACCORDING TO THE FORMAT: (I1,1X,F4.2,2X,I4,3X,F5.2). * C* * C* 4-ONE CARD CONTAINING THE FOLLOWING SCALE FACTOR:- * C* SCALXY: TO SCALE THE LOCAL RECTANGULAR COORDINATES X & Y; * C* ACCORDING TO THE FORMAT: (D10.3). * C* * C* 5-ONE DATA CARD PER EACH ONE-BLOCK-SQUARED TO BE INVESTIGATED, * C* CONTAINING THE LATITUDE AND LONGITUDE OF ITS CENTRE, AS FOLLOWS:- * C* COLS 1- 8 : LATITUDE (+VE NORTH) ,IN DECIMAL DEGREES............F8.2 * C* COLS 9-16 : LONGITUDE (+VE WEST) ,IN DECIMAL DEGREES............F8.2 * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* THE OUTPUT WILL BE :- * C* ------------------ * C* * C* 1-IDENTIFICATION OF THE LOCATION OF THE SQUARED-AREA IN QUESTION. * C* 2-RESULTS OF THE LEAST-SQUARES APPROXIMATION TO THE FREE-AIR ANOMALY * C* FIELD - AS OBTAINED FROM THE SUBROUTINE 'ORTHO'. * C* 3-SAME AS (2) ABOVE, BUT FOR THE HEIGHT FIELD. * C* 4-RESULTS AND PLOTS OF THE GRAVITY CORRECTIONS VS. AZIMUTH, IN CASE * C* OF HELMERT, VIGNAL AND DYNAMIC SYSTEMS OF HEIGHTS - AS OBTAINED FROM * C* THE SUBROUTINE 'GCAFAZ'. * C* 5-A PLOT SHOWING THE DISTRIBUTION OF DATA POINTS WITHIN THE ONE-SQUARED * C* AREA , AS OBTAINED FROM THE SUBROUTINE 'GCPLOT'. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* IT SHOULD BE NOTED THAT THIS PROGRAM CAN HANDLE ONE-SQUARED-AREA ONLY * C* AT A TIME, WHICH IMPLIES THAT THE NUMBER OF AREA-SQUARES THAT CAN BE * C* PROCESSED BY THE PROGRAM IS, THEORITICALLY, UNLIMITED. ON THE OTHER HAND* C* ,THIS WILL BE RESTRICTED BY THE AVAILABLE COMPUTER TIME. * C* * C* IT SHOULD BE NOTED HERE THAT A COMMON BLOCK CONTAINING THE X, Y LOCAL * C* COORDINATES FOR ALL THE DATA POINTS WITHIN THE SPECIFIED AREA, AS WELL * C* AS THE ORDER OF THE APPROXIMATING POLYNOMIAL FOR ANOMALY AND HEIGHT * C* FIELDS, MUST BE DECLARED BETWEEN THE MAIN PROGRAM AND THE FUNCTION- * C* SUBPROGRAM 'PHI'. * C* * C* FINALLY, SEVERAL ERROR MESSAGES MAY BE PRINTED-OUT WHEN NECESSARY. THE * C* USER IS NOTIFIED OF THE NON-EXISTENCE OF GRAVITY DATA ON FILES FOR A * C* PARTICULAR CHOSEN BLOCK; OF ANY INSUFFIENCY OR POOR DISTRIBUTION OF * C* DATA POINTS ; * C* THESE MESSAGES WILL BE SELF-EXPLANATORY AND READY TO GUIDE THE USER. * C* C*************************************************************************** C IMPLICIT REAL*8(A-H,O-Z) C C SPECIFY THE DIMENSIONS OF THE INVOLVED ARRAYS C C XFY : THE X-PLANE-COORDINATE (+VE NORTH), IN KILOMETERS. C YAL : THE Y-PLANE-COORDINATE (+VE EAST) , IN KILOMETERS. C DGFA : THE FREE-AIR ANOMALY, IN MGALS. C WFA : THE WEIGHT FUNCTION OF THE FREE-AIR ANOMALY C HT : THE ELEVATION OF THE POINT ABOVE M.S.L., IN METERS. C WHT : THE WEIGHT FUNCTION OF THE ELEVATIONS (HEIGHTS) C CIJ : THE COEFFICIENTS OF THE APPROXIMATING POLYNOMIAL FOR THE C FREE-AIR ANOMALY SURFACE. C SIGCIJ: THE VARIANCE-COVARIANCE MATRIX OF 'CIJ' C TIJ : THE COEFFICIENTS OF THE APPROXIMATING POLYNOMIAL FOR THE C TOPOGRAPHY (HEIGHT) SURFACE. C SIGTIJ: THE VARIANCE-COVARIANCE MATRIX OF 'TIJ' . C ALPHAS: WORKING SQUARE MATRIX FOR THE ORTHOGONALIZATION PROCESS. C CFOUR : FOURIER COEFFICIENTS (FROM ORTHO SUBROUTINE). C SC2FOR: VECTOR OF SQUARES OF NORMS OF FOURIER COEFFICIENTS. C IQUADR: THE QUADRANT INDICATOR. C DIMENSION XFY(3000),YAL(3000) DIMENSION DGFA(3000),WFA(3000) DIMENSION HT(3000),WHT(3000) DIMENSION CIJ(16),SIGCIJ(16,16),TIJ(16),SIGTIJ(16,16) DIMENSION ALPHAS(16,16),CFOUR(16),SC2FOR(16) DIMENSION IQUADR(4) INTEGER POLDEG DOUBLE PRECISION DSIN,DCOS,DSQRT,DABS,DSIGN MROWS=16 C C DECLARE EXPLICITLY THE LABELLED COMMON BLOCK BETWEEN THIS MAIN-PROGRAM C AND THE FUNCTION-SUBPROGRAM 'A' (WHICH WILL EVALUATE THE INDIVIDUAL C ELEMENTS OF THE DESIGN MATRIX-A OF THE BASE-FUNCTIONS OF THE POLYNOMIAL). C COMMON /XYMAIN/ XFY,YAL,POLDEG C C READ-IN THE PRINTING CODE 'LSPRNT' ,FOR THE INTERMEDIATE RESULS OF THE C LEAST-SQUARES APPROXIMATION TECHNIQUE. C LSPRNT=0 : NO PRINTING IS REQUIRED. C LSPRNT>0 : PRINTING IS NEEDED. C READ(5,85) LSPRNT 85 FORMAT(I1) C C READ-IN THE FOLLOWING INFORMATION CONCERNING THE AREA OF INTEREST: C POLDEG: ORDER OF APPROXIMATING POLYNOMIAL FOR ANOMALY & HEIGHT FIELDS. C EXTENT: THE SIZE OF THE AREA (IN DEGREES) IN BOTH LAT. & LONG. DIRECTIONS. C MAXDIM: THE MAXIMUM DIMENSIONAL SPACE ALLOCATED FOR THE INVOLVED ARRAYS. C DIST : THE DESIRED INTEGRATION DISTANCE (IN KM.) ALONG LEVELLING LINES. C READ(5,86) POLDEG,EXTENT,MAXDIM,DIST 86 FORMAT(I1,1X,F4.2,2X,I4,3X,F5.2) C C READ-IN THE NECESSARY SCALE FACTOR , AS FOLLOWS:- C SCALXY: TO BE MULTIPLIED BY THE LINEAR UNITS: SPECIFICALLY: THE X & Y PLANE C COORDINATES, THE INTEGRATION LENGTH 'DIST' OF THE LEVELLING LINE C AND THE MEAN RADIUS OF CURVATURE 'R0M' OF THE ELLIPSOID, IN ORDER TO C GET A WELL-CONDITIONED MATRIX-Q OF THE NORMAL EQUATIONS.FOR INVERSION C READ(5,87) SCALXY 87 FORMAT(D10.3) C C COMPUTE THE TOTAL NUMBER OF UNKNOWN COEFFICIENTS (WHICH EQUALS THE C NUMBER OF BASE-FUNCTIONS) FOR THE APPROXIMATING POLYNOMIAL . C MPOLY=(POLDEG+1)*(POLDEG+1) C C SET-UP VALUES FOR SOME CONVERSION CONSTANTS. C PI=3.141592653589793D 0 RHODEG=180.D 0/PI RHOSEC=RHODEG*3600.D 0 FEETMS=0.3048D 0 C C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C 999 CONTINUE C C READ-IN THE LATITUDE (+VE N) AND THE LONGITUDE (+VE W ) IN DECIMAL DEGREES C OF THE CENTRE OF THE DESIRED AREA-SQUAREA -UNDER INVESTIGATION C NOTE THAT THIS CENTRE WILL BE REFERRED TO AS THE POINT OF INTEREST C READ(5,31,END=99) PHI,ALONG 31 FORMAT(2F6.2) C C CONVERT THE LATITUDE OF THE GIVEN POINT OF INTEREST TO RADIANS C PHI0=PHI/RHODEG C C PRINT-OUT SOME HEADINGS FOR THE JOB AND IDENTIFICATION FOR THE AREA C SQUARE UNDER INVESTIGATION. C WRITE(6,32) 32 FORMAT('1',//,10X,'INVESTIGATION OF THE ACTUAL GRAVITY CORRECTION 1TO THE HEIGHT DIFFERENCES -BASED ON NORMAL GRAVITY-'/,10X,'VERSUS 2THE AZIMUTH OF THE LEVELLING LINE -MEASURED CLOCKWISE FROM THE CEN 3TRE OF THE ONE-BLOCK-SQUARED'//) WRITE(6,33) PHI,ALONG,EXTENT 33 FORMAT(10X,'THE CENTRE OF THE ONE-SQUARED-BLOCK UNDER INVESTIGATIO 1N IS LOCATED AT :'/,15X,'LATITUDE :',F7.2,2X,'DEGREES (NORTH)'/15 2X,'LONGITUDE :',F7.2,2X,'DEGREES (WEST)'/,10X,'THE SIDES OF THIS S 3QUARED-AREA ARE',F5.2,' DEGREES IN BOTH LATITUDE & LONGITUDE'//,40 4X,'=============================='//) C C DEFINE THE PARAMETERS (SEMI-MAJOR & SEMI-MINOR AXES) OF THE REFERENCE C ELLIPSOID IN CANADA, THE 1927-NAD OF CLARK 1866 ELLIPSOID.***IN KILOMETERS**** C ANAD=6378.2064D 0 BNAD=6356.5838D 0 ECC1SQ=(ANAD**2-BNAD**2)/ANAD**2 C C COMPUTATION OF THE MEAN RADIUS OF CURVATURE (IN KILOMETERS) OF THE NAD- C ELLIPSOID, AT THE LATITUDE (PHI-DEGREES) OF THE GIVEN POINT OF INTEREST. C DENO=1.0D 0-ECC1SQ*DSIN(PHI0)**2 RCN=ANAD/DSQRT(DENO) RCM=RCN*(1.D 0-ECC1SQ)/DENO R0M=DSQRT(RCN*RCM) C C COMPUTE THE CONVERSION FACTORS FROM CURVILINEAR DISTANCES ON THE ELLIPSOID C IN THE LATITUDE AND LONGITUDE DIRECTIONS, TO THE CORRESPONDING LINEAR C DISTANCES ON THE PLANE IN X AND Y DIRECTIONS, RESPECTIVELY. C FYFACT=R0M/RHODEG ALFACT=DCOS(PHI0)*FYFACT C C SET-UP THE ALLOWABLE INTERPOLATION BOUNDARIES (IN BOTH LATITUDE AND C LONGITUDE) WITHIN EACH OF THE SEVEN PARTITIONED DATA SETS. C PHIL21=41.5D 0 PHIU21=63.5D 0 ALGL21=52.5D 0 ALGU21=73.5D 0 PHIL22=43.5D 0 PHIU22=63.5D 0 ALGL22=73.5D 0 ALGU22=91.5D 0 PHIL23=48.5D 0 PHIU23=63.5D 0 ALGL23=91.5D 0 ALGU23=118.5D 0 PHIL24=48.5D 0 PHIU24=63.5D 0 ALGL24=118.5D 0 ALGU24=140.5D 0 PHIL25=63.5D 0 PHIU25=70.5D 0 ALGL25=61.5D 0 ALGU25=142.5D 0 PHIL16=70.5D 0 PHIU16=80.5D 0 ALGL16=66.5D 0 ALGU16=142.5D 0 PHIL17=80.5D 0 PHIU17=82.5D 0 ALGL17=51.5D 0 ALGU17=78.5D 0 C C TESTING THE GIVEN LOCATION (LONGITUDE & LATITUDE) OF THE BENCH MARK C AGAINST THE ABOVE SPECIFIED INTERPOLATION BOUNDARIES , TO DETERMINE THE C APPROPRIATE FILE (WHOSE DATA REFERENCE NUMBER IS DENOTED BY 'IUNIT'), C THAT CONTAINS THE NECESSARY DATA REQUIRED FOR INTERPOLATION. C IF(ALONG.GE.ALGL21.AND.ALONG.LE.ALGU21) GO TO 1 IF(ALONG.GE.ALGL22.AND.ALONG.LE.ALGU22) GO TO 2 IF(ALONG.GE.ALGL23.AND.ALONG.LE.ALGU23) GO TO 3 IF(ALONG.GE.ALGL24.AND.ALONG.LE.ALGU24) GO TO 4 C C PRINT-OUT A MESSAGE TO THE USER IF THE INPUTTED BENCH MARK IS LOCATED C OUTSIDE THE RANGE OF OBSERVED GRAVITY DATA EXISTING ON THE AVAILABLE C DATA FILES. C WRITE(6,10) 10 FORMAT('1'///,5X,'SORRY: THE GIVEN POINT IS LOCATED IN AN AREA WHE 1RE NO OBSERVED GRAVITY DATA EXIST ON THE AVAILABLE FILES'//) GO TO 9999 1 IF(PHI.GE.PHIL21.AND.PHI.LE.PHIU21) GO TO 11 IF(PHI.GE.PHIL25.AND.PHI.LE.PHIU25) GO TO 15 IF(PHI.GE.PHIL16.AND.PHI.LE.PHIU16) GO TO 6 IF(PHI.GE.PHIL17.AND.PHI.LE.PHIU17) GO TO 7 WRITE(6,10) GO TO 9999 2 IF(PHI.GE.PHIL22.AND.PHI.LE.PHIU22) GO TO 12 IF(PHI.GE.PHIL25.AND.PHI.LE.PHIU25) GO TO 15 IF(PHI.GE.PHIL16.AND.PHI.LE.PHIU16) GO TO 6 IF(PHI.GE.PHIL17.AND.PHI.LE.PHIU17) GO TO 7 WRITE(6,10) GO TO 9999 3 IF(PHI.GE.PHIL23.AND.PHI.LE.PHIU23) GO TO 13 IF(PHI.GE.PHIL25.AND.PHI.LE.PHIU25) GO TO 15 IF(PHI.GE.PHIL16.AND.PHI.LE.PHIU16) GO TO 6 WRITE(6,10) GO TO 9999 4 IF(PHI.GE.PHIL24.AND.PHI.LE.PHIU24) GO TO 14 IF(PHI.GE.PHIL25.AND.PHI.LE.PHIU25) GO TO 15 IF(PHI.GE.PHIL16.AND.PHI.LE.PHIU16) GO TO 6 WRITE(6,10) GO TO 9999 C C SPECIFY THE APPROPRIATE DATA-SET REFERENCE-NUMBER ,WHICH CONTAINS THE C NECESSARY GRAVITY DATA AROUND THE POINT OF INTEREST. C 11 IUNIT=21 GO TO 5 12 IUNIT=22 GO TO 5 13 IUNIT=23 GO TO 5 14 IUNIT=24 GO TO 5 15 IUNIT=25 GO TO 5 6 IUNIT=16 GO TO 5 7 IUNIT=17 5 CONTINUE C C SET-UP THE MAX. & MIN. LIMITS (IN BOTH LAT. & LONG.) -ALL IN DEGREES- FOR C A CHOSEN ONE-AREA-SQUARED AREA SURROUNDING THE POSITION OF THE GIVEN POINT C OF INTEREST,THE GRAVITY POINTS (STORED ON DISC FILES) WITHIN THIS AREA C WILL BE USED IN PERFORMING THE INTERPOLATION TECHNIQUE C EXTD=EXTENT DISD=DIST C 960 CONTINUE C AREA=EXTD/2. FYMAX=PHI+AREA FYMIN=PHI-AREA ALMAX=ALONG+AREA ALMIN=ALONG-AREA C C INITIALIZE ZERO VALUES FOR THE 4-QUADRANTS INDICATOR 'IQUADR' C DO 40 I=1,4 40 IQUADR(I)=0 C C READ-IN THE DATA RESIDING ON DISC AS DATA SET NO. 'IUNIT' CONTAINING C FOR EACH POINT:- LAT.+VE N AND LONG.+VE W ,BOTH IN DECIMAL DEGREES ; C OBSERVED GRAVITY IN GALS ; FREE-AIR GRAVITY ANOMAL IN MGALS ; ELEVATION C AND ITS ACCURACY IN FEET. C REWIND IUNIT NP=1 50 READ(IUNIT,END=60) ALAT,ALOG,GRAV,DG,H,HACC C C CHECK, SELECT, AND STORE THE INFORMATION FOR ONLY THOSE POINTS LOCATED C WITHIN THE SPECIFIED SQUARED-AREA AROUND THE POINT OF INTEREST. C IF(ALAT.LT.FYMIN.OR.ALAT.GT.FYMAX) GO TO 50 IF(ALOG.LT.ALMIN.OR.ALOG.GT.ALMAX) GO TO 50 C C COMPUTE AND STORE THE X (+VE NORTH) AND Y (+VE EAST) PLANE ORTHOGONAL C COORDINATES OF EACH SELECTED POINT-IN KILOMETERS- ,REFERRED TO THE POINT C OF INTEREST AS THE ORIGIN OF SUCH LOCAL PLjNE COORDINATE SYSTEM. C XFY(NP)=(ALAT-PHI)*FYFACT YAL(NP)=(ALONG-ALOG)*ALFACT C C SCALE THE X AND Y LOCAL PLANE-COORDINATES OF EACH SELECTED POINT, IN UNITS C OF (1/SCALXY)-KM TO HELP GETTING A WELL-CONDITIONED MATRIX-Q OF THE NORMAL- C EQUATIONS. IT SHOULD BE NOTED THAT THIS SCALING WILL BE EQUIVALENT TO C A VARIABLE SCALING OF EACH COLUMN IN THE DESIGN MATRIX-A (VANDERMONDE'S) C XFY(NP)=XFY(NP)*SCALXY YAL(NP)=YAL(NP)*SCALXY C C SET-UP AND STORE THE ELEVATION OF EACH SELECTED POINT ,ALONG WITH ITS C WEIGHT FUNCTION (INVERSE OF ITS VARIANCE) - IN METER UNITS C HT(NP)=H*FEETMS VARHT=(HACC*FEETMS+1.0D-5)**2 WHT(NP)=1.0D 0/VARHT C C SET-UP AND STORE THE FREE-AIR GRAVITY-ANOMALY ,AT EACH SELECTED POINT, C ALONG WITH ITS WEIGHT FUNCTION (INVERSE OF ITS VARIANCE, WHICH IS COMPUTED C FROM THE VARIANCE OF THE ELEVATION) -IN MGAL UNITS C DGFA(NP)=DG WFA(NP)=1.0D 0/(0.05D 0**2+0.3086D 0**2*VARHT) C C STORE ANY INTEGER NUMBER, SAY 1, IN THE POSITION 'JACK' OF THE QUADRANT- C INDICATOR ARRAY 'IQUADR', TO HELP TESTING FOR THE EXISTANCE OF AT LEAST C ONE DATA POINT IN EACH QUADRANT AROUND THE POINT OF INTEREST C C 'JACK' TAKES THE VALUES: 1,2,3 OR 4 ,ACCORDING TO THE QUADRANTS WITH C X,Y=(-,-),(-,+),(+,-) AND (+,+) RESPECTIVELY. C JACK=2.05D 0+DSIGN(1.0D 0,XFY(NP))+(1.0D 0+DSIGN(1.0D 0,YAL(NP)))/ 12.0D 0 IQUADR(JACK)=1 C C UPDATE THE TOTAL-NuMBER OF ACCUMULATED SELECTED DATA POrNTS WITHIN-THE AREA C NP=NP+1 C C TEST IF THE NUMBER OF DATA POINTS FOUND EXCEEDS THE MAXIMUM ALLOCATED C STORAGE FOR IT ,IN THE DIMENSION STATEMENT. C IF(NP.GT.MAXDIM) GO TO 60 GO TO 50 60 CONTINUE C C COMPUTE AND PRINT-OUT THE TOTAL NUMBER OF AVAILABLE DATA POINTS C NP=NP-1 WRITE(6,71) NP 71 FORMAT(10X,'THE TOTAL NUMBER OF AVAILABLE DATA POINTS IS =',I5,2X, 1'POINTS'/) C******************************************************************************* C C CHECK IF THE NUMBER OF DATA POINTS AVAILABLE IN THE CHOSEN AREA ,IS C SUFFICIENT FOR ESTIMATING THE COEFFICIENTS OF THE APPROXIMATING POLYNOMIAL C (WHOSE NUMBER IS = MPOLY), OTHERWISE PRINT-OUT A MESSAGE TO THE USER , C AND THEN STOP EXECUTION. C IF(NP.LE.MPOLY) GO TO 980 C C TEST FOR THE EXISTANCE OF AT LEAST ONE DATA POINT IN EACH OF THE FOUR C QUADRANTS AROUND THE POINT OF INTEREST, OTHERWISE PRINT-OUT A MESSAGE TO C THE USER, AND THEN STOP EXECUTION. C IQVAL=IQUADR(1)*IQUADR(2)*IQUADR(3)*IQUADR(4) IF(IQVAL.EQ.0) GO TO 990 C C TEST IF THE NUMBER OF DEGREES OF FREEDOM IS MORE THAN 50% OF THE TOTAL C NUMBER OF UNKNOWN POLYNOMIAL-COEFFICIENTS, OTHERWISE GIVE A WARNING TO C THE USER, AND THEN CONTINUE EXECUTION AND FIND THE RESULTS. C NUMDF=NP-MPOLY IF(NUMDF.LT.MPOLY/2) WRITE(6,65) 65 FORMAT(///,10X,'WARNING:- BE CAREFUL WHEN USING AND ANALYZING THE 1FOLLOWING RESULTS'/,10X,'SINCE THE NUMBER OF DEGREES OF FREEDOM IS 2 LESS THAN 50% OF THE UNKNOWNS'/) C******************************************************************************* C C CALL THE SUBROUTINE 'ORTHO' TO COMPUTE AND PRINt-OUT THE ESTIMATED C COEFFICIENTS OF THE BEST-FITTING POLYNOMIAL, APPROXIMATING THE FREE-AIR C GRAVITY SURFACE, ALONG WITH THEIR ASSOCIATED VARIANCE-COVARIANCE MATRIX. C WRITE(6,333) 333 FORMAT(3X,'***THE FOLLOWING RESULTS PERTAIN TO THE BEST-FITTING PO 1LYNOMIAL FOR THE FREE-AIR GRAVITY-ANOMALY SURFACE***'/) C CALL ORTHO(NP,MPOLY,DGFA,WFA,MROWS,CIJ,SIGCIJ,ALPHAS,CFOUR,SC2FOR, 1LSPRNT) C C EXTRACT AND PRINT-OUT THE PREDICTED VALUE OF THE FREE-AIR ANOMALY, ALONG C WITH ITS ESTIMATED STANDARD DEVIATION, AT THE CENTRE OF THE ONE-AREA- C SQUARED UNDER INVESTIGATION (POINT OF INTEREST), FOR VISUAL EXAMINATION. C DGPRED=CIJ(1) DGPSTD=DSQRT(SIGCIJ(1,1)) WRITE(6,512) DGPRED,DGPSTD 512 FORMAT(//5X,'THE PREDICTED FREE-AIR ANOMALY AT THE CENTRE OF THE G 1IVEN ONE-AREA-SQUARED IS =',F10.4,' MGALS'/,10X,'WITH AN ESTIMAT 2ED STANDARD DEVIATION =',F10.4,' MGALS'//) C C******************************************************************************* C C CALL ALSO THE SAME SUBROUTINE 'ORTHO' MENTIONED ABOVE ,TO COMPUTE C AND PRINT-OUT THE ESTIMATED COEFFICIENTS OF THE BEST-FITTING POLYNOMIAL C APPROXIMATING THE TOPOGRAPHY (HEIGHT) SURFACE, ALONG WITH THEIR C ASSOCIATED VARIANCE-COVARIANCE MATRIX. C C C WRITE(6,777) WRITE(6,444) 444 FORMAT(3X,'***THE FOLLOWING RESULTS PERTAIN TO THE BEST-FITTING PO 1LYNOMIAL FOR THE TOPOGRAPHY (HEIGHT) SURFACE***'/) C CALL ORTHO(NP,MPOLY,HT,WHT,MROWS,TIJ,SIGTIJ,ALPHAS,CFOUR,SC2FOR,LS 1PRNT) C C EXTRACT AND PRINT-OUT THE PREDICTED VALUE OF THE HEIGHT, ALONG WITH C ITS STANDARD DEVIATION, AT THE CENTRE OF THE ONE-AREA-SQUARED IN QUESTION C HTPRED=TIJ(1) HTPSTD=DSQRT(SIGTIJ(1,1)) WRITE(6,513) HTPRED,HTPSTD 513 FORMAT(//5X,'THE PREDICTED HEIGHT AT THE CENTRE OF THE GIVEN ONE-A 1REA-SQUARED IS =',F10.4,' METERS'/,10X,'WITH AN ESTIMATED STANDA 2RD DEVIATION =',F10.4,' METERS'//) C C******************************************************************************* C C CALL THE SUBROUTINE 'GCAFAZ' TO COMPUTE, PRINT-OUT, AND PLOT THE GRAVITY- C CORRECTION (TO THE HEIGHT DIFFERENCE) AS A FUNCTION OF THE AZIMUTH OF A C HYPOTHETICAL RUNNING LEVELLING LINE OF 1 KM LONG (AVERAGE LENGTH OF A C SECTION FOR FIRST-ORDER WORK). C IT SHOULD BE NOTED HERE THAT THE PLOTTED GRAVITY-CORRECTIONS ARE C STANDARDIZED FOR 1 KM., BY DIVIDING THE CORRESPONDING ACCUMULATED C CORRECTIONS BY THE SQUARE-ROOT OF THE TOTAL ACCUMULATED DISTANCE IN KM. C THIS IS TO BE CONSISTENT WITH THE ADOPTED FORMULA FOR LEVELLING C ACCURACY SPECIFICATIONS C WRITE(6,777) WRITE(6,555) 555 FORMAT(/,3X,'***THE FOLLOWING INFORMATION PERTAIN TO THE GRAVITY-C 1ORRECTION TO THE HEIGHT-DIFFERENCE VS. AZIMUTH***'/) C CALL GCAFAZ(CIJ,TIJ,PHI0,R0M,DISD,SCALXY,MPOLY) C C******************************************************************************* C C CALL THE SUBROUTINE 'GCPLOT' TO PLOT A SKETCH SHOWING THE APPROXIMATE C DISTRIBUTION OF THE USED DATA POINTS WITHIN THE ONE-AREA-SQUARED C SPECIFIED AROUND THE POINT OF INTEREST, FOR THE INTERPOLATION PURPOSES. C PLOTFY= AREA *FYFACT*SCALXY PLOTAL= AREA *ALFACT*SCALXY XFYMIN=-PLOTFY YALMIN=-PLOTAL XFYMAX=+PLOTFY YALMAX=+PLOTAL C CALL GCPLOT(YAL,XFY,NP,YALMIN,XFYMIN,YALMAX,XFYMAX) C CMFLSC=( EXTD*60.)/15. WRITE(6,554) CMFLSC 554 FORMAT(5X,'SCALE OF THE ABOVE PLOT IS: 1.0 CENTIMETER REPRESENTS', 1F5.2,2X,'MINUTES (IN BOTH LAT.& LONG.)'/) WRITE(6,552) 552 FORMAT(10X,'*****DISTRIBUTION OF DATA POINTS WITHIN THE ONE-AREA S 1QUARED*****'/) WRITE(6,33) PHI,ALONG,EXTD GO TO 9999 C**************R**************************************************************** C 980 WRITE(6,981) 981 FORMAT(///,10X,'ERROR MESSAGE: THE LEAST SQUARES ESTIMATION FOR TH 1E POLYNOMIAL COEFFICIENTS IS IMPOSSIBLE'/,25X,'BECAUSE THE NUMBER 2OF AVAILABLE DATA POINTS IS LESS THAN THE NUMBER OF',/,25X,'UNKNOW 3N COEFFICIENTS (-VE DEGREES OF FREEDOM)') C C DOUBLE THE SIZE OF THE GIVEN ORIGINAL BLOCK, AS WELL AS THE INTEGRATION C DISTANCE FOR THE GRAVITY CORRECTIONS., AND PRINT-OUT THE NEW SIZE OF BLOCK C EXTD=EXTD*2. DISD=DISD*2. WRITE(6,961) EXTD 961 FORMAT(/,5X,'***THE SIDES OF THE CURRENT SQUARED-BLOCK ARE =',F5.2 1,' DEGREES***'/) IF(EXTD.LE.1.0) GO TO 960 GO TO 9999 990 WRITE(6,991) 991 FORMAT(///,10X,'ERROR MESSAGE: THE CHOSEN SIZE OF THE AREA AROUND 1THE ABOVE POINT OF INTEREST'/,25X,'DOES NOT HAVE A RELIABLE DISTRI 2BUTION OF AT LEAST ONE DATA POINT IN EACH QUADRANT') C C DOUBLE THE SIZE OF THE GIVEN ORIGINAL BLOCK, AS WELL AS THE INTEGRATION C DISTANCE FOR THE GRAVITY CORRECTIONS., AND PRINT-OUT THE NEW SIZE OF BLOCK C EXTD=EXTD*2. DISD=DISD*2. WRITE(6,961) EXTD IF(EXTD.LE.1.0) GO TO 960 777 FORMAT('1') C 9999 GO TO 999 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* 99 STOP END C C C*************************************************************************** C* * C* SUBROUTINE 'G C A F A Z' :- ( BY : M.M. NASSAR, JANUARY 1975.) * C* ------------------------ * C* * C* THE PURPOSE OF THIS SUBROUTINE IS TO COMPUTE* C* AND PRINT-OUT THE ACCUMULATED ACTUAL GRAVITY-CORRECTIONS (IN MM.) ALONG * C* A SIMULATED LEVELLING LINE ( S KM. LONG), RADIATING FROM THE CENTRE OF * C* A GIVEN SQUARED-BLOCK AS VARRYING WITH THE AZIMUTH OF THE LINE, AT * C* FIVE DEGREES INTERVALS. * C* * C* THE COMPUTED GRAVITY-CORRECTIONS ARE THEN STANDARDIZED PER 1 KM. OF THE * C* LEVELLING LINE (BY DIVIDING BY THE SQUARE-ROOT OF ITS LENGTH IN KM.). * C* THREE SYSTEMS OF HEIGHTS: HELMERT; VIGNAL; AND DYNAMIC, ARE CONSIDERED * C* UNDER INVESTIGATION BY THIS SUBROUTINE. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* NOTE THAT SUCH GRAVITY CORRECTIONS HAVE BEEN FORMULATED TO BE ADDED * C* THE CORRESPONDING HEIGHT DIFFERENCES BASED ON NORMAL GRAVITY ONLY , AND * C* CURRENTLY BEING USED IN CANADA. MOREOVER, THESE CORRECTIONS HERE , ARE * C* BASED ON A THIRD-ORDER TWO-DIMENSIONAL APPROXIMATING POLYNOMIAL TO THE * C* FREE-AIR ANOMALY AND HEIGHT FIELDS, AND FURTHER ON A FOURTH-ORDER ONE- * C* DIMENSIONAL APPROXIMATING PPOLYNOMIALS TO THE ANOMALY AND HEIGHT PROFILE* C* ALONG THE SIMULATED LEVELLING LINES. * 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* USAGE :- CALL GCAFAZ(CIJ,TIJ,PHI0,R0M,S,SC,MCT) * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* REQUIRED INPUT DATA :- (VIA THE CALLING ROUTINE) * C* ------------------ * C* * C* CIJ , TIJ : VECTORS OF COEFFICIENTS (16 EACH) OF APPROXIMATING POLYNOM- * C* IALS TO THE FREE-AIR GRAVITY ANOMALY AND HEIGHT FIELDS * C* WITHIN THE ONE BLOCK-SQUARE AREA, RESPECTIVELY. * C* PHI0, R0M : THE LATITUDE (IN RADIANS) OF, AND THE MEAN RADIUS OF CURVAT-* C* URE (IN KM.) OF THE EARTH AT, THE CENTRE OF THE ONE BLOCK- * C* SQUARE, RESPECTIVELY. * C* S : THE DESIRED INTEGRATION DISTANCE (IN KM.) ALONG THE SIMULAT-* C* ED LEVELLING LINE (E.G., 40 KM.IN CASE OF ONE-DEGREE-SQUARE)* C* SC : THE SUITABLE SCALE FACTOR APPLIED TO THE XY COORDINATE * C* SYSTEM TO GET A WELL-CONDITIONED NORMAL EQUATIONS MATRIX, * C* WHICH MUST BE APPLIED ALSO HERE ON R0M AND S. * C* MCT : THE TOTAL NUMBER OF COEFFICIENTS OF THE APPROXIMATING * C* POLYNOMIAL TO THE FREE-AIR ANOMALY AND HEIGHT SURFACES * C* (SIZE OF BOTH 'CIJ' & 'TIJ' ARRAYS). * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* THE OUTPUT WILL BE :- * C* ------------------ * C* * C* 1-TABLE CONTAINING THE ACCUMULATED AS WELL AS THE STANDARDIZED GRAVITY- * C* CORRECTIONS VS. AZIMUTH, IN CASE OF HELMERT, VIGNAL AND DYNAMIC * C* SYSTEMS OF HEIGHTS, FOR SIMULATED LEVELLING LINES AT 5-DEGREES * C* INTERVALS IN AZIMUTH. * C* 2-SUMMARY OF RESULTS CONTAINING THE MAXIMUM AND MINIMUM ABSOLUTE VALUES * C* OF THE STANDARDIZED GRAVITY-CORRECTIONS (MM/KM.) , AS WELL AS THE * C* RESPECTIVE DIRECTIONS (AZIMUTH) IN WHICH THEY OCCUR- AGAIN FOR THE * C* THREE SYSTEMS OF HEIGHTS UNDER INVESTIGATION. * C* 3-A SEPARATE PLOT FOR EACH OF THE THREE SYSTEMS CONSIDERED, SHOWING * C* THE PATTERN OF THE GRAVITY CORRECTION VS. AZIMUTH. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* NOTE THAT THIS SUBROUTINE CAN HANDLE ONLY ONE-BLOCK-SQUARED AT A TIME. * C* * C*************************************************************************** C SUBROUTINE GCAFAZ(CIJ,TIJ,PHI0,R0M,S,SC,MCT) IMPLICIT REAL*8(A-H,O-Z) DIMENSION CIJ(MCT),TIJ(MCT) DIMENSION C(5),T(5),ODASH(4) DIMENSION O(8),V(8),D(8) DIMENSION XO(72),YO(72),XV(72),YV(72),XD(72),YD(72) DOUBLE PRECISION DSIN,DCOS,DSQRT,DABS DOUBLE PRECISION DMAX1,DMIN1 C C SET-UP VALUES FOR SOME CONVERSION CONSTANTS. C PI=3.141592653589793D 0 RHODEG=180.D 0/PI RHOSEC=RHODEG*3600.D 0 CONST5=0.2238D 0 C C SCALING THE DISTANCE 'S', WHICH IS GIVEN IN KILOMETERS, TO THE SAME C SCALE-FACTOR 'SC' OF THE X & Y PLANE-COORDINATES, TO GET CONSISTENT C UNITS AND MEANINGFUL RESULTS. ALSO THE SAME SCALING SHOULD BE APPLIED C TO THE MEAN RADIUS OF CURVATURE OF THE ELLIPSOID 'R0M'. C SS=S*SC R0MS=R0M*SC C C COMPUTE THE SQUARE-ROOT OF THE GIVEN INTEGRATION DISTANCE 'S-KM.' (BEFORE C ANY SCALING). THIS WILL BE NEEDED FOR THE STANDARDIZATION PROCESS OF THE C COMPUTED GRAVITY-CORRECTIONS- FOR 1 KM. ONLY. C DISTN=DSQRT(S) C C COMPUTATION OF SOME CONSTANT COEFFICIENTS -COMMON TO ALL SYSTEMS OF C HEIGHTS IN QUESTION- WHICH ARE FUNCTIONS OF THE LATITUDE (PHI0) OF THE C CENTRE OF EACH ONE-DEGREE-SQUARE. THESE COEFFICIENTS WILL BE IN MILLIGALS. C GREF45=980624.0D 0 A0=-6.295D 0 A1=-0.358D 0 A2=1.076D 0 SPHI0=DSIN(PHI0) S2PHI0=DSIN(2.D 0*PHI0) S4PHI0=DSIN(4.D 0*PHI0) CONST1=A1*S2PHI0-2.D 0*A2*S4PHI0 CONST2=A0-A1*SPHI0**2+A2*S2PHI0**2 C C INITIALIZING THE AZIMUTH OF THE LINE RADIATING FROM THE BLOCK-CENTRE C ALPHA=0.D 0 C C COMPUTATIONS OF THE GRAVITY-CORRECTIONS TO ALL THE THREE HEIGHT SYSTEMS C UNDER INVESTIGATION, FOR DIFFERENT AZIMUTHS WITH FIVE-DEGREES INTERVALS. C C INITIALIZE SOME REASONABLE VALUES - TO HELP EXTRACTING THE ABSOLUTE- C MAXIMUM AND THE ABSOLUTE-MINIMUM VALUES OF THE GRAVITY-CORRECTION. C ZOMAX=0.0D 0 ZVMAX=0.0D 0 ZDMAX=0.0D 0 ZOMIN=1.0D 03 ZVMIN=1.0D 03 ZDMIN=1.0D 03 IORGIN=0 C C PRINT-OUT HEADINGS FOR THE TABLE OF COMPUTED GRAVITY-CORRECTIONS VS. C AZIMUTH, FOR THE THREE HEIGHT-SYSTEMS UNDER INVESTIGATION. C WRITE(6,1) 1 FORMAT(/,54X,'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',16X,'VIGNAL',15X 2,'DYNAMIC',/,7X,'FROM',3X,'TO',7X,'(KM)',3X,'(DEGREES)',5X,'(MM)', 36X,'(MM/KM)',5X,'(MM)',6X,'(MM/KM)',5X,'(MM)',6X,'(MM/KM)'/) C C***********S T A R T OF LOOP FOR GRAVITY-CORRECTION**************************** C DO 100 KAZ=1,72 AZ=ALPHA/RHODEG DS=DSIN(AZ) DC=DCOS(AZ) CONST4=DC/R0MS C C COMPUTATION OF THE COEFFICIENTS (C'S) OF THE BEST FITTING POLYNOMIAL TO THE C FREE-AIR GRAVITY ANOMALY PROFILE ALONG THE LEVELLING LINE ,AS FUNCTIONS OF C THE AZIMUTH OF THE LINE AND THE COEFFICIENTS (CIJ'S) OF THE BEST FITTING C POLYNOMIAL TO THE FREE-AIR GRAVITY-ANOMALY SURFACE WITHIN THE ONE-BLOCK- C SQUARE BLOCK IN QUESTION (OBTAINED FROM THE INTERPOLATION SUBROUTINE) C -'A P P R O X'. C C(1)=CIJ(1) C(2)=CIJ(2)*DS+CIJ(5)*DC C(3)=CIJ(3)*DS**2+CIJ(6)*DS*DC+CIJ(9)*DC**2 C(4)=CIJ(4)*DS**3+CIJ(7)*DC*DS**2+CIJ(10)*DS*DC**2+CIJ(13)*DC**3 C(5)=CIJ(8)*DC*DS**3+CIJ(11)*(DS*DC)**2+CIJ(14)*DS*DC**3 C C COMPUTATION OF THE COEFFICIENTS (T'S) OF THE BEST FITTING POLYNOMIAL TO THE C TOPOGRAPHY (HEIGHT) PROFILE ALONG THE LEVELLING LINE , AS FUNCTIONS OF THE C AZIMUTH OF THE LINE AND THE COEFFICIENTS (TIJ'S) OF THE BEST FITTING C POLYNOMIAL TO THE TOPOGRAPHY (HEIGHT) SURFACE WITHIN THE ONE-BLOCK- C SQUARE BLOCK IN QUESTION (OBTAINED FROM THE INTERPOLATION SUBROUTINE) C -'A P P R O X'. C T(1)=TIJ(1) T(2)=TIJ(2)*DS+TIJ(5)*DC T(3)=TIJ(3)*DS**2+TIJ(6)*DS*DC+TIJ(9)*DC**2 T(4)=TIJ(4)*DS**3+TIJ(7)*DC*DS**2+TIJ(10)*DS*DC**2+TIJ(13)*DC**3 T(5)=TIJ(8)*DC*DS**3+TIJ(11)*(DS*DC)**2+TIJ(14)*DS*DC**3 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C COMPUTATION OF HELMERT-GRAVITY-CORRECTION ( DELTA-HELMERT ) FOR THE C GIVEN LINE WITH LENGTH S KM, FROM THE INTEGRAL FORMULA DERIVED BY C SUPERIMPOSING BOTH THE FREE-AIR ANOMALY PROFILE AND THE TOPOGRAPHY PROFILE C ALONG THE LINE. C C 1-COMPUTE THE AUXILLIARY COEFFICIENTS(O- TELDA) C ODASH(1)=-(C(2)-CONST1*CONST4-CONST5*T(2))/GREF45 ODASH(2)=-(C(3)-CONST5*T(3))/GREF45 ODASH(3)=-(C(4)-CONST5*T(4))/GREF45 ODASH(4)=-(C(5)-CONST5*T(5))/GREF45 C C 2-COMPUTE THE EXPLICIT COEFFICIENTS (O'S) OF THE DIFFERENTIAL FORMULA C OF HELMERT GRAVITY-CORRECTION. C O(1)=T(1)*ODASH(1) O(2)=T(2)*ODASH(1)+2.D 0*T(1)*ODASH(2) O(3)=T(3)*ODASH(1)+2.D 0*T(2)*ODASH(2)+3.D 0*T(1)*ODASH(3) O(4)=T(4)*ODASH(1)+2.D 0*T(3)*ODASH(2)+3.D 0*T(2)*ODASH(3)+4.D 0*T 1(1)*ODASH(4) O(5)=T(5)*ODASH(1)+2.D 0*T(4)*ODASH(2)+3.D 0*T(3)*ODASH(3)+4.D 0*T 1(2)*ODASH(4) O(6)=2.D 0*T(5)*ODASH(2)+3.D 0*T(4)*ODASH(3)+4.D 0*T(3)*ODASH(4) O(7)=3.D 0*T(5)*ODASH(3)+4.D 0*T(4)*ODASH(4) O(8)=4.D 0*T(5)*ODASH(4) C C 3-SUBSTITUTE IN THE EXPLICIT INTEGRAL FORMULA FOR THE HELMERT-GRAVITY- C CORRECTION TO OBTAIN ITS VALUE (IN METERS) FOR A LINE OF S KM LENGTH C DELTAO=O(1)*SS+O(2)*SS**2/2.0D 0+O(3)*SS**3/3.0D 0+O(4)*SS**4/4.0D 1 0+O(5)*SS**5/5.0D 0+O(6)*SS**6/6.0D 0+O(7)*SS**7/7.0D 0+O(8)*SS** 28/8.0D 0 C C COMPUTE 'DELTA-HELMERT' ,IN MILLIMETERS. C DELTAO=DELTAO*1000.D 0 C C COMPUTE THE STANDARDIZED 'DELTA-HELMERT' ,IN MILLIMETERS. C NOTE THAT THE STANDARDIZATION (NORMALIZATION) IS DONE HERE FOR 1 KM. C BY DIVIDING THE ACCUMULATED GRAVITY-CORRECTION BY THE SQUARE-ROOT OF C THE INTEGRATION DISTANCE (S), I.E. BY DIVIDING BY 'DISTN' COMPUTED BEFORE. C DELTON=DELTAO/DISTN C C COMPUTE THE MAXIMUM AND MINIMUM ABSOLUTE VALUES OF 'DELTON' FOR THE C PURPOSE OF PRINTING-OUT A SUMMARY OF THE RESULT. C ALONG WITH THEIR RESPECTIVE AZIMUTHb. C ZOMAX=DMAX1(ZOMAX,DABS(DELTON)) IF(ZOMAX.EQ.DABS(DELTON)) AZOMAX=ALPHA ZOMIN=DMIN1(ZOMIN,DABS(DELTON)) IF(ZOMIN.EQ.DABS(DELTON)) AZOMIN=ALPHA C C COMPUTE AND STORE THE X (HORIZONTAL) AND Y (VERTICAL) PLANE COMPONENTS OF C NORMALIZED HELMERT-GRAVITY-CORRECTION (IN MILLIMETERS), AS A FUNCTION OF C THE AZIMUTH ONLY ,SINCE THEY CORRESPOND TO 1 KM LENGTH. C THIS IS FOR THE PURPOSE OF PLOTTING THESE GRAV. CORREC. VS. AZIMUTH. C XO(KAZ)=DABS(DELTON)*DS YO(KAZ)=DABS(DELTON)*DC C C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C COMPUTATION OF VIGNAL--GRAVITY-CORRECTION 'DELTA-VIGNAL' FOR THE C GIVEN LINE WITH LENGTH S KM, FROM THE INTEGRAL FORMULA DERIVED BY C SUPERIMPOSING BOTH THE FREE-AIR ANOMALY PROFILE AND THE TOPOGRAPHY PROFILE C ALONG THE LINE. C C 1-COMPUTE THE EXPLICIT COEFFICIENTS (V'S) OF THE DIFFERENTIAL FORMULA (AS A C FUNCTION OF DISTANCE) FOR THE VIGNAL-GRAVITY-CORRECTION C CONST6=CONST1*CONST4 V(1)=(T(2)*C(1)+CONST6*T(1))/GREF45 V(2)=(T(2)*(C(2)+CONST6)+T(3)*2.D 0*C(1))/GREF45 V(3)=(T(2)*C(3)+T(3)*(2.D 0*C(2)+CONST6)+T(4)*3.D 0*C(1))/GREF45 V(4)=(T(2)*C(4)+T(3)*2.D 0*C(3)+T(4)*(3.D 0*C(2)+CONST6)+T(5)*4.D 10*C(1))/GREF45 V(5)=(T(2)*C(5)+T(3)*2.D 0*C(4)+T(4)*3.D 0*C(3)+T(5)*(4.D 0*C(2)+C 1ONST6))/GREF45 V(6)=(T(3)*2.D 0*C(5)+T(4)*3.D 0*C(4)+T(5)*4.D 0*C(3))/GREF45 V(7)=(T(4)*3.D 0*C(5)+T(5)*4.D 0*C(4))/GREF45 V(8)=(T(5)*4.D 0*C(5))/GREF45 C C 2-SUBSTITUTE IN THE EXPLICIT INTEGRAL FORMULA FOR THE VIGNAL-GRAVITY- C CORRECTION TO OBTAIN ITS VALUE (IN METERS) FOR A LINE OF S KM LENGTH C DELTAV=V(1)*SS+V(2)*SS**2/2.0D 0+V(3)*SS**3/3.0D 0+V(4)*SS**4/4.0D 1 0+V(5)*SS**5/5.0D 0+V(6)*SS**6/6.0D 0+V(7)*SS**7/7.0D 0+V(8)*SS** 28/8.0D 0 C C COMPUTE 'DELTA-VIGNAL' ,IN MILLIMETERS. C DELTAV=DELTAV*1000.D 0 C C COMPUTE THE STANDARDIZED 'DELTA-VIGNAL' ,IN MILLIMETERS. C NOTE THAT THE STANDARDIZATION (NORMALIZATION) IS DONE HERE FOR 1 KM. C BY DIVIDING THE ACCUMULATED GRAVITY-CORRECTION BY THE SQUARE-ROOT OF C THE INTEGRATION DISTANCE (S), I.E. BY DIVIDING BY 'DISTN' COMPUTED BEFORE. C DELTVN=DELTAV/DISTN C C COMPUTE THE MAXIMUM AND MINIMUM ABSOLUTE VALUES OF 'DELTVN' FOR THE C PURPOSE OF PRINTING-OUT A SUMMARY OF THE RESULT. C ALONG WITH THEIR RESPECTIVE AZIMUTHS. C ZVMAX=DMAX1(ZVMAX,DABS(DELTVN)) IF(ZVMAX.EQ.DABS(DELTVN)) AZVMAX=ALPHA ZVMIN=DMIN1(ZVMIN,DABS(DELTVN)) IF(ZVMIN.EQ.DABS(DELTVN)) AZVMIN=ALPHA C C COMPUTE AND STORE THE X (HORIZONTAL) AND Y (VERTICAL) PLANE COMPONENTS OF THE C NORMALIZED VIGNAL-GRAVITY-CORRECTION (IN MILLIMETERS), AS A FUNCTION OF C THE AZIMUTH ONLY ,SINCE THEY CORRESPOND TO 1 KM LENGTH. C THIS IS FOR THE PURPOSE OF PLOTTING THESE GRAV. CORREC. VS. AZIMUTH. C XV(KAZ)=DABS(DELTVN)*DS YV(KAZ)=DABS(DELTVN)*DC C C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C C COMPUTATION OF THE DYNAMIC-GRAVITY-CORRECTION 'DELTA-DYNAMIC', FOR THE C GIVEN LINE WITH LENGTH S KM, FROM THE INTEGRAL FORMULA DERIVED BY C SUPERIMPOSING BOTH THE FREE-AIR ANOMALY PROFILE AND THE TOPOGRAPHY PROFILE C ALONG THE LINE. C C 1-COMPUTE THE EXPLICIT COEFFICIENTS (D'S) OF THE DIFFERENTIAL FORMULA (AS A C FUNCTION OF DISTANCE) FOR THE DYNAMIC-GRAVITY-CORRECTION. C CONST7=CONST2+C(1) D(1)=(T(2)*CONST7)/GREF45 D(2)=(T(2)*C(2)+T(3)*2.0D 0*CONST7)/GREF45 D(3)=(T(2)*C(3)+T(3)*2.0D 0*C(2)+T(4)*3.0D 0*CONST7)/GREF45 D(4)=(T(2)*C(4)+T(3)*2.0D 0*C(3)+T(4)*3.0D 0*C(2)+T(5)*4.0D 0*CONS 1T7)/GREF45 D(5)=(T(2)*C(5)+T(3)*2.0D 0*C(4)+T(4)*3.0D 0*C(3)+T(5)*4.0D 0*C(2) 1)/GREF45 D(6)=(T(3)*2.D 0*C(5)+T(4)*3.D 0*C(4)+T(5)*4.D 0*C(3))/GREF45 D(7)=(T(4)*3.D 0*C(5)+T(5)*4.D 0*C(4))/GREF45 D(8)=(T(5)*4.D 0*C(5))/GREF45 C C 2-SUBSTITUTE IN THE EXPLICIT INTEGRAL FORMULA FOR THE DYNAMIC-GRAVITY- C CORRECTION TO OBTAIN ITS VALUE (IN METERS) FOR A LINE OF S KM LENGTH C DELTAD=D(1)*SS+D(2)*SS**2/2.0D 0+D(3)*SS**3/3.0D 0+D(4)*SS**4/4.0D 1 0+D(5)*SS**5/5.0D 0+D(6)*SS**6/6.0D 0+D(7)*SS**7/7.0D 0+D(8)*SS** 28/8.0D 0 C C COMPUTE 'DELTA-DYNAMIC' ,IN MILLIMETERS. C DELTAD=DELTAD*1000.D 0 C C COMPUTE THE STANDARDIZED 'DELTA-DYNAMIC' ,IN MILLIMETERS. C NOTE THAT THE STANDARDIZATION (NORMALIZATION) IS DONE HERE FOR 1 KM. C BY DIVIDING THE ACCUMULATED GRAVITY-CORRECTION BY THE SQUARE-ROOT OF C THE INTEGRATION DISTANCE (S), I.E. BY DIVIDING BY 'DISTN' COMPUTED BEFORE. C DELTDN=DELTAD/DISTN C C COMPUTE THE MAXIMUM AND MINIMUM ABSOLUTE VALUES OF 'DELTDN' FOR THE C PURPOSE OF PRINTING-OUT A SUMMARY OF THE RESULT. C ALONG WITH THEIR RESPECTIVE AZIMUTHS. C ZDMAX=DMAX1(ZDMAX,DABS(DELTDN)) IF(ZDMAX.EQ.DABS(DELTDN)) AZDMAX=ALPHA ZDMIN=DMIN1(ZDMIN,DABS(DELTDN)) IF(ZDMIN.EQ.DABS(DELTDN)) AZDMIN=ALPHA C C COMPUTE AND STORE THE X (HORIZONTAL) AND Y (VERTICAL) PLANE COMPONENTS OF THE C NORMALIZED DYNAMIC-GRAVITY-CORRECTION (IN MILLIMETERS), AS A FUNCTION OF C THE AZIMUTH ONLY ,SINCE THEY CORRESPOND TO 1 KM LENGTH. C THIS IS FOR THE PURPOSE OF PLOTTING THESE GRAV. CORREC. VS. AZIMUTH. C XD(KAZ)=DABS(DELTDN)*DS YD(KAZ)=DABS(DELTDN)*DC C C PRINT-OUT THE FOLLOWING INFORMATION FOR EACH SIMULATED LEVELLING LINE C RADIATING FROM THE CENTRE OF THE CELL:- SERIAL NUMBER OF ITS END-POINTS C ; LENGTH (IN KM); AZIMUTH (IN DEGREES); HELMERT GRAVITY-CORRECTION (IN MM) C ; STANDARDIZED HELMERT GRAVITY CORRECTION (IN MM/KM ); VIGNAL GRAVITY- C CORRECTION (IN MM); STANDARDIZED VIGNAL GRAVITY CORRECTION (IN MM/KM ); C DYNAMIC GRAVITY-CORRECTION (IN MM); AND THE STADARDIZED DYNAMIC GRAVITY C CORRECTION ( IN MM/KM ). C IF(KAZ.EQ.37) WRITE(6,3) 3 FORMAT('1'///) IF(KAZ.EQ.37) WRITE(6,1) WRITE(6,2) IORGIN,KAZ,S,ALPHA,DELTAO,DELTON,DELTAV,DELTVN,DELTAD,D 1ELTDN 2 FORMAT(9X,I1,4X,I2,7X,F4.1,5X,F5.1,4X,F9.4,3X,F8.4,2X,F9.4,3X,F8.4 1,2X,F9.4,3X,F8.4) C C UPDATE THE AZIMUTH OF THE RUNNING LINE FROM THE CENTER OF THE BLOCK , BY C ADDING THE DESIRED INCREMENT (E.G. HERE WE ADD 5 DEGREES). C ALPHA=ALPHA+5.D 0 100 CONTINUE C C*******END OF LOOP FOR GRAVITY-CORRECTION******************** C C PRINT-OUT A SUMMARY OF THE MAXIMUM & MINIMUM ABSOLUTE VALUES OF THE C GRAVITY-CORRECTION FOR EACH OF THE THREE DIFFERENT HEIGHT-SYSTEMS IN C QUESTION (HELMERT'S, VIGNAL'S, & DYNAMIC) FOR THE ONE-BLOCK -SQUARE C UNDER INVESTIGATION. C WRITE(6,22) 22 FORMAT('1'///35X,'***S U M M A R Y O F R E S U L T S***'/,38X,'- 1--------------------------------'/) WRITE(6,24) ZOMAX,ZOMIN 24 FORMAT(06X,'THE MAXIMUM-ABSOLUTE VALUE OF THE NORMALIZED HELMERT-G 1RAVITY-CORRECTION IS =',F10.4,2X,' MM / KM '/,06X,'THE MINIMUM-ABS 2OLUTE VALUE OF THE NORMALIZED HELMERT-GRAVITY-CORRECTION IS =',F10 3.4,2X,' MM / KM '/) WRITE(6,10) AZOMAX,AZOMIN 10 FORMAT(15X,'AT CORRESPONDING AZIMUTHS OF ',F8.2,' AND',F8.2,' D 1EGREES, RESPECTIVELY'//) WRITE(6,26) ZVMAX,ZVMIN 26 FORMAT(06X,'THE MAXIMUM-ABSOLUTE VALUE OF THE NORMALIZED VIGNAL-G 1RAVITY-CORRECTION IS =',F10.4,2X,' MM / KM '/,06X,'THE MINIMUM-ABS 2OLUTE VALUE OF THE NORMALIZED VIGNAL-GRAVITY-CORRECTION IS =',F10 3.4,2X,' MM / KM '/) WRITE(6,10) AZVMAX,AZVMIN WRITE(6,28) ZDMAX,ZDMIN 28 FORMAT(06X,'THE MAXIMUM-ABSOLUTE VALUE OF THE NORMALIZED DYNAMIC-G 1RAVITY-CORRECTION IS =',F10.4,2X,' MM / KM '/,06X,'THE MINIMUM-ABS 2OLUTE VALUE OF THE NORMALIZED DYNAMIC-GRAVITY-CORRECTION IS =',F10 3.4,2X,' MM / KM '/) WRITE(6,10) AZDMAX,AZDMIN C*-*-*-*-*-*-*-*-*-*-*-*-**-*-*-**--*-*-*-*-*-*-*--*-*-*-*-*--*-*---*-*-*-*-*-*- C C CALL THE SUBROUTINE 'GCPLOT' FOR PLOTTING THE NORMALIZED GRAVITY- C CORRECTIONS, AS A FUNCTION OF THE AZIMUTH, FOR THE THREE DIFFERENT HEIGHT- C SYSTEMS UNDER INVESTIGATION. C XOMIN=-ZOMAX YOMIN=-ZOMAX XOMAX=+ZOMAX YOMAX=+ZOMAX C CALL GCPLOT(XO,YO,72,XOMIN,YOMIN,XOMAX,YOMAX) C CMOSC=XOMAX/7.5D 0 WRITE(6,21) CMOSC 21 FORMAT(5X,'SCALE OF THE ABOVE PLOT IS: 1.0 CENTEMETER REPRESENTS', 1F8.4,2X,'MM / KM ,OF GRAVITY-CORRECTION') WRITE(6,8) 8 FORMAT(//,5X,'***** VARIATION OF HELMERT-GRAVITY-CORRECTION WITH A 1ZIMUTH *****'/) C XVMIN=-ZVMAX YVMIN=-ZVMAX XVMAX=+ZVMAX YVMAX=+ZVMAX C CALL GCPLOT(XV,YV,72,XVMIN,YVMIN,XVMAX,YVMAX) C CMVSC=XVMAX/7.5D 0 WRITE(6,21) CMVSC WRITE(6,9) 9 FORMAT(//,5X,'***** VARIATION OF VIGNAL-GRAVITY-CORRECTION WITH A 1ZIMUTH *****'/) C XDMIN=-ZDMAX YDMIN=-ZDMAX XDMAX=+ZDMAX YDMAX=+ZDMAX C CALL GCPLOT(XD,YD,72,XDMIN,YDMIN,XDMAX,YDMAX) C CMDSC=XDMAX/7.5D 0 WRITE(6,21) CMDSC WRITE(6,11) 11 FORMAT(//,5X,'***** VARIATION OF DYNAMIC-GRAVITY-CORRECTION WITH A 1ZIMUTH *****'/) RETURN END C C C*************************************************************************** C* * C* SUBROUTINE 'G C P L O T' :- ( BY : M.M. NASSAR, JANUARY 1975.) * C* ------------------------ * C* * C* THE PURPOSE OF THIS SUBROUTINE IS TO PLOT * C* TWO PERPENDICULAR AXES: X (HORIZONTAL) AND Y (VERTICAL), INTERSECTING * C* AT THE CENTRE OF A SQUARE 6 BY 6 INCHES; AND THEN PLOT 'N' NUMBER OF * C* POINTS DEFINED BY THEIR X & Y COORDINATES RELATIVE TO THE CENTRE OF * C* THE SQUARE AND GIVEN IN THE USER'S UNITS. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* NOTE THAT THE TRANSFORMATION OF COORDINATES FROM THE USER UNITS TO THE * C* PLOTTING UNITS IS DONE AUTOMATICALLY VIA THE PLOTTING ROUTINES CALLED * C* BY THIS SUBROUTINE. FOR MORE DETAILS CONCERNING THE PLOTTING ROUTINES* C* AVAILABLE AT THE U.N.B. COMPUTING CENTRE, SEE: ''COMPUTER PLOTTING'', * C* A USER GUIDE MANUAL BY: U.G. GUJAR, CL # 278, COMPUTING CENTRE, U.N.B. * C* FREDERICTON, N.B., OCTOBER, 1972. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* USAGE :- CALL GCPLOT(X,Y,N,XMIN,YMIN,XMAX,YMAX) * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* REQUIRED INPUT DATA :- (VIA THE CALLING ROUTINE) * C* ------------------- * C* * C* X : VECTOR OF X COORDINATES OF THE N-POINTS, IN USER'S UNITS. * C* Y : VECTOR OF Y COORDINATES OF THE N-POINTS, IN USER'S UNITS. * C* N : TOTAL NUMBER OF POINTS TO BE PLOTTED. * C* XMIN : MINIMUM X COORDINATE IN THE USER'S UNITS. * C* YMIN : MINIMUM Y COORDINATE IN THE USER'S UNITS. * C* XMAX : MAXIMUM X COORDINATE IN THE USER'S UNITS. * C* YMAX : MAXIMUM Y COORDINATE IN THE USER'S UNITS. * C* * C* THE OUTPUT WILL BE THE PLOTTED SQUARE, THE TWO PERPENDICULAR AXES, AND * C* THE N-POINTS. FINALLY, THIS SUBROUTINE PRINTS-OUT THE TOTAL NUMBER 'N' * C* OF PLOTTED POINTS NO LIMITATIONS ON THE NUMBER OF POINTS. * C* * C*************************************************************************** C SUBROUTINE GCPLOT(X,Y,N,XMIN,YMIN,XMAX,YMAX) IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(N),Y(N) XLNGTH=6. YLNGTH=6. CALL DEVICE(1627) CALL AREA(XLNGTH,YLNGTH) CALL SETPLT(XMIN,YMIN,XMAX,YMAX) CALL PRNTCH('+') RECXL=2.0D 0*XMIN CALL RECT(XMIN,YMIN,RECXL) CALL AXES(XMIN,YMIN,XMAX,YMAX,0.,0.,-1.,-1.) CALL PRNTCH('.') DO 30 K=1,N CALL NOWPLT(0,X(K),Y(K)) 30 CONTINUE CALL ENDPLT WRITE(6,19) N 19 FORMAT(/,5X,'TOTAL NUMBER OF PLOTTED POINTS SURROUNDING THE CENTRA 1L-ORIGIN IS =',I6,2X,'POINTS') RETURN END C C*************************************************************************** C C THIS SUBROUTINE ORTHOGONALIZES THE MATRIX PHI USING THE GRAM-SCHMIDT C METHOD, COMPUTES THE FOURIER COEFFICIENTS OF THE ORTHOGONALIZED MATRIX 02 C DERIVES THE COEFFICIENTS OF PHI,COMPUTES THE VARIANCES OF THE FOURIER 0277 C COEFFICIENTS AND THE VARIANCE-COVARIANCE MATRIX OF THE COEFFICIENTS 0278 C INPUTS : C 1. PHI(OPTIONAL - COULD BE FUNCTION SUBPROGRAM INSTEAD) - AN N BY M 0280 C CONTAINING THE BASE FUNCTIONS EVALUATED FOR EACH OBSERVATION C 2. N - THE NUMBER OF OBSERVATIONS C 3. M - THE NUMBER OF BASE FUNCTIONS (>2) C 4. W - WEIGHT FUNCTIONS C 5. F - FUNCTIONAL VALUES C 6. NRD - THE MAXIMUM ROW DIMENSION OF PHI C 7. MRD- THE MAXIMUM COLUMN DIMENSION OF PHI MATRIX. C C OUTPUTS : C 1. ALPHA - AN MRD BY M MATRIX CONTAINING THE ALPHA'S USED IN COMPUTI C THE ORTHOGONALIZED MATRIX AND IN COMPUTING THE COEFFICIENTS OF PH C 2. C - THE M FOURIER COEFFICIENTS OF THE ORTHOGONALIZED MATRIX C 3. D - THE M COEFFICIENTS OF THE INPUT MATRIX PHI C 4. SUMD - THE VARIANCE-COVARIANCE MATRIX OF THE COEFFICIENTS C 5. SC2 - THE SQUARES OF THE NORMS OF THE ORTHOGONALIZED MATRIX C 6. V - THE N RESIDUALS C C*************************************************************************** C SUBROUTINE ORTHO(N,M,F,W,MRD,D,SUMD,ALPHA,C,SC2,LSPRNT) IMPLICIT REAL*8(A-H,O-Z) DIMENSION W(N),F(N),D(M),SUMD(MRD,M) DIMENSION ALPHA(MRD,M),C(M),SC2(M) DIMENSION IHEAD(8) C C TEST FOR NEGATIVE DEGREES OF FREEDOM C IF (N.LT.M) GO TO 100 K=1 ALPHA(M,M)=1.D0 C C DETERMINE THE ALPHA'S FOR COMPUTATION OF ORTHOGONALIZED MATRIX C 10 DO 3 J=K,M IF(J.NE.K) GO TO 6 ALPHA(K,K)=1.D0 GO TO 3 6 SC1=0.D0 SC2(K)=0.D0 SC3=0.D0 DO 2 I=1,N P=PHI(I,K) IF(K.EQ.1) GO TO 4 K1=K-1 DO 5 J1=1,K1 5 P=P+ALPHA(J1,K)*PHI(I,J1) 4 SC1=SC1+W(I)*PHI(I,J)*P SC3=SC3+F(I)*W(I)*P 2 SC2(K)=SC2(K)+W(I)*P**2 ALPHA(J,K)=-SC1/SC2(K) ALPHA(K,J)=ALPHA(J,K) 3 CONTINUE C C DETERMINE THE FOURIER COEFFICIENTS FOR THE ORTHOGONALIZED MATRIX C C(K)=SC3/SC2(K) K=K+1 IF(K.LT.3) GO TO 10 C C DETERMINE THE ALPHA'S USED IN COMPUTING THE COEFFICIENTS OF PHI C JK=K-1 9 JL=K JK=JK-1 JJ=K-JK-1 DO 8 LM=1,JJ JL=JL-1 8 ALPHA(JK,K)=ALPHA(JK,K)+ALPHA(JK,JL)*ALPHA(K,JL) IF(JK.NE.1) GO TO 9 IF(K.LT.M) GO TO 10 C C DETERMINE THE LAST FOURIER COEFFICIENT C SC2(K)=0.D0 SC3=0.D0 DO 7 I=1,N P=PHI(I,K) K1=K-1 DO 1 J=1,K1 1 P=P+ALPHA(J,K)*PHI(I,J) SC2(K)=SC2(K)+W(I)*P**2 7 SC3=SC3+F(I)*W(I)*P C(K)=SC3/SC2(K) C C DETERMINE THE COEFFICIENTS OF PHI C DO 13 I=1,M D(I)=C(I) IF(I.EQ.M) GO TO 13 II=I+1 DO 14 J=II,M 14 D(I)=D(I)+ALPHA(I,J)*C(J) 13 CONTINUE C C PRINT-OUT THE ESTIMATED COEFFICIENTS OF THE APPROXIMATING POLYNOMIAL. C IF(LSPRNT.EQ.0) GO TO 51 WRITE(6,52) M 52 FORMAT(//,5X,'VECTOR OF ESTIMATED POLYNOMIAL COEFFICIENTS'/15X,'DI 1MENSIONS : (',I3,' , 1 )') WRITE(6,53) 53 FORMAT(/,15X,'1'/) DO 40 JJ=1,M WRITE(6,41) JJ,D(JJ) 41 FORMAT(3X,I4,D15.5) 40 CONTINUE 51 CONTINUE C C COMPUTE THE VARIANCE OF THE FOURIER COEFFICIENTS AND THE VARIANCE-COVA 0 C MATRIX OF THE COEFFICIENTS C DO 15 I=1,M DO 15 J=1,M 15 SUMD(I,J)=0.D0 C C PRINT-OUT THE VECTOR OF RESIDUALS. C IF(LSPRNT.EQ.0) GO TO 54 WRITE(6,888) WRITE(6,55) N 55 FORMAT(//,5X,'VECTOR OF RESIDUALS'/,3X,'DIMENSIONS : (',I4,' , 1 ) 1') WRITE(6,53) 54 CONTINUE SC4=0.D0 DO 22 I=1,N PN=0.D0 DO 21 J=1,M 21 PN=PN+D(J)*PHI(I,J) C C COMPUTE THE RESIDUALS, SQUARE-SUM OF THE WEIGHTED RESIDUALS AND THE C A POSTERIORI VARIANCE-FACTOR. C V1=F(I)-PN IF(LSPRNT.EQ.0) GO TO 56 WRITE(6,41) I,V1 56 CONTINUE V2=V1**2 22 SC4=SC4+V2*W(I) SIGMA2=SC4/(N-M) C C PRINT-OUT THE NUMBER OF DEGREES OF FREEDOM, SQUARE-SUM OF THE WEIGHTED- C RESIDUALS AND THE ESTIMATED VARIANCE-FACTOR (A POSTERIORI). C NUMDF=N-M WRITE(6,43) NUMDF,SC4,SIGMA2 43 FORMAT(//,5X,'DEGREES OF FREEDOM =',I6/,5X,'SQUARE-SUM OF WEIGHTED 1-RESIDUALS =',D15.5,/,5X,'A POSTERIORI VARIANCE-FACTOR =',D15.5/) C C COMPUTE THE VARIANCE-COVARIANCE MATRIX OF THE ORIGINAL COEFFICIENTS. C DO 23 I=1,M SUMC=SIGMA2/SC2(I) DO 23 J=1,I DO 23 K=J,I 23 SUMD(J,K)=SUMD(J,K)+ALPHA(J,I)*ALPHA(K,I)*SUMC DO 24 I=1,M IT=I+1 IF(IT.GT.M) GO TO 30 DO 24 J=IT,M 24 SUMD(J,I)=SUMD(I,J) C 30 CONTINUE IF(LSPRNT.EQ.0) GO TO 57 C C SPECIFY AND INITIALIZE SOME OF THE REQUIRED PARAMETERS FOR PRINTING- C OUT THE VAR.-COV. MATRIX OF THE ESTIMATED PO YNOMIAL COEFFICIENTS. C 888 FORMAT('1') DO 50 I=1,8 50 IHEAD(I)=0 GROUP8=M/8. M8TH=GROUP8 MTOTAL=M8TH RESTM=GROUP8-M8TH IF(RESTM.GE.0.12) MTOTAL=M8TH+1 C C PRINT-OUT THE VARIANCE-COVARIANCE MATRIX OF THE POLYNOMIAL COEFFICIENTS. C WRITE(6,888) WRITE(6,91) M,M 91 FORMAT(/,5X,'THE ESTIMATED VARIANCE-COVARIANCE MATRIX OF THE ESTIM 1ATED POLYNOMIAL COEFFICIENTS'/,25X,'DIMENSIONS : (',I3,' ,',I3,' ) 2') MSTART=1 MEND=8 DO 151 JIM=1,MTOTAL IHEAD(1)=MSTART DO 160 J=2,8 160 IHEAD(J)=IHEAD(J-1)+1 WRITE(6,82) (IHEAD(J),J=1,8) 82 FORMAT(/,10X,8(I2,11X),/) DO 150 I=1,M 150 WRITE(6,83) I,(SUMD(I,J),J=MSTART,MEND) 83 FORMAT(I4,8D13.4) MSTART=MSTART+8 MEND=MEND+8 IF(MEND.GT.M) MEND=M 151 CONTINUE 57 CONTINUE RETURN C 100 WRITE(6,101) 101 FORMAT('0','*ERROR* NEGATIVE DEGREES OF FREEDOM') RETURN END C C*************************************************************************** C* * C* FUNCTION-SUBPROGRAM ' PHI' :- * C* --------------------------- * C* * C* THIS FUNCTION-SUBPROGRAM IS DESIGNED AS A * C* GENERAL FUNCTION, TO COMPUTE (IN DOUBLE PRECISION) SINGLE ELEMENTS * C* PHI(L,K) OF THE DESIGN MATRIX (WHICH IS ALSO KNOWN AS VANDERMONDE'S * C* MATRIX), NEEDED IN THE FORMATION OF THE NORMAL EQUATIONS SYSTEM IN THE * C* LEAST-SQUARES APPROXIMATION TECHNIQUE- USING A TWO-DIMENSIONAL MIXED- * C* ALGEBRAIC POLYNOMIAL OF ANY DESIRED ORDER (NK) , WHERE L IS THE ROW * C* POSITION AND K IS THE COLUMN POSITION OF THE SOUGHT ELEMENT PHI(L,K) * C* * C* NOTE THAT THE MAIN IDEA BEHIND THIS FUNCTION IS TO SAVE A LARGE AMOUNT * C* OF COMPUTER STORAGE, DUE TO THE FACT THAT MATRIX-PHI DOES NOT NEED TO BE* C* GENERATED IN FULL, SAY WITH DIMENSIONS (N BY M), AND INSTEAD WE COMPUTE * C* ONLY ONE OF ITS ELEMENTS AT A TIME WHEN REQUIRED. * C* * C* FOR MORE DETAILS CONCERNING THE MATHEMATICAL BASIS OF THIS FUNCTION, * C* SEE: 'INTRODUCTION TO APPROXIMATION THEORY' BY: E.W. CHENEY, 1966, * C* MCGRAW-HILL BOOK COMPANY, NEW YORK. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* USAGE:- THIS FUNCTION WILL BE AUTOMATICALLY CALLED (TO EVALUATE PHI(L,K)* C* AS SOON AS THE EXECUTION IN THE CALLING ROUTINE HITS ANY * C* EXPRESSION THAT INVOLVES ANY INDIVIDUAL ELEMENT PHI(L,K) OF PHI * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* THE REQUIRED INPUT DATA ARE THE VECTORS OF X & Y LOCAL ORTHOGONAL * C* COORDINATES OF ALL THE DATA POINTS AVAILABLE IN THE INTERPOLATION * C* AREA (E. G. ONE-DEGREE-SQUARE), AND THE DESIRED ORDER OF THE APPROXIM- * C* ATING POLYNOMIAL.. THIS DATA MUST BE SPECIFIED AS A LABELLED COMMON * C* BLOCK BETWEEN THIS FUNCTION-SUBPROGRAM AND THE MAIN ROUTINE THAT * C* COMPUTES THE X,Y COORDINATES.. * C* * C* THE OUTPUT WILL BE THE NUMERICAL VALUE OF PHI(L,K) TRANSFERRED TO THE * C* RESPECTIVE EXPRESSION IN THE CALLING ROUTINE . * C* * C* IT SHOULD BE NOTED THAT ALTHOUGH THERE ARE NO LIMITATIONS PUT ON THE * C* VALUE OF THE ORDER OF THE APPROXIMATING POLYNOMIAL, THE NUMBER OF DATA * C* POINTS IS NOT TO EXCEED 3000 , OTHERWISE THE CORRESPONDING DIMENSION * C* STATEMENTS SHOULD BE MODIFIED,.. * C* * C*************************************************************************** C C PHI COMPUTES THE ELEMENTS OF THE DESIGN MATRIX FOR ORTHO C DOUBLE PRECISION FUNCTION PHI(L,K) DOUBLE PRECISION X,Y DIMENSION X(3000),Y(3000) COMMON /XYMAIN/ X,Y,NK NN=NK+1 I=(K-1)/NN J=K-I-(NN-1)*I-1 IF(X(L).EQ.0.0.AND.I.EQ.0) X(L)=1.0 IF(Y(L).EQ.0.0.AND.J.EQ.0) Y(L)=1.0 PHI=X(L)**I*Y(L)**J RETURN END