C C C*************************************************************************** C* * C* PROGRAM 'L E V A G R A V' : (BY: M.M. NASSAR, JANUARY 1975.) * C* ------------------------- * 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). * 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 PURPOSE OF THE MAIN PROGRAM 'LEVAGRAV' IS TO READ THE LATITUDE * C* AND LONGITUDE OF THE CENTRE OF THE ONE-DEGREE-SQUARE UNDER INVESTIGATION* 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 INVERSION PROCESS. * C* * C* THE LEAST-SQUARES PREDICTION OF THE ANOMALY AND HEIGHT FIELDS IS * C* PERFORMED VIA THE SUBROUTINE 'APPROX' * 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* ONE-DEGREE-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-ONE DATA CARD PER EACH ONE-DEGREE-SQUARE 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* 2-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* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* THE OUTPUT WILL BE :- * C* ------------------ * C* * C* 1-IDENTIFICATION OF THE LOCATION OF THE DEGREE-SQUARE IN QUESTION. * C* 2-RESULTS OF THE LEAST-SQUARES APPROXIMATION TO THE FREE-AIR ANOMALY * C* FIELD - AS OBTAINED FROM THE SUBROUTINE 'APPROX' * 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-DEGREE- * C* SQUARE, AS OBTAINED FROM THE SUBROUTINE 'GCPLOT'. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* IT SHOULD BE NOTED THAT THIS PROGRAM CAN HANDLE ONE-DEGREE-SQUARE ONLY * C* AT A TIME, WHICH IMPLIES THAT THE NUMBER OF DEGREE-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* THE FOLLOWING PARAMETERS HAVE TO BE SPECIFIED AT THE BEGINNING OF THE * C* MAIN PROGRAM DECK, NAMELY: THE MAXIMUM ALLOWABLE NUMBER OF DATA POINTS, * C* THE ORDER OF THE APPROXIMATING POLYNOMIAL, THE SCALE FACTOR FOR THE * C* COORDINATE SYSTEM,THE SCALE FACTOR FOR THE NORMAL EQUATIONS MATRIX, * C* AND THE DESIRED INTEGRATION DISTANCE (IN KM.). ALSO, A COMMON BLOCK * C* CONTAINING THE X,Y LOCAL COORDINATES FOR ALL THE DATA POINTS, AS WELL * C* AS THE ORDER OF THE POLYNOMIAL, MUST BE DECLARED BETWEEN THE MAIN * C* PROGRAM AND THE FUNCTION-SUBPROGRAM 'A'. * 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 ; OR OF THE S&NGULA-&TY OF THE NO-MAL EQUAT&ONS MATRIX. * 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 IQUADR: THE QUADRANT INDICATOR. C U : THE ABSOLUTE TERMS OF THE NORMAL EQUATIONS. C DIMENSION XFY(2000),YAL(2000) DIMENSION DGFA(2000),WFA(2000) DIMENSION HT(2000),WHT(2000) DIMENSION CIJ(16),SIGCIJ(16,16),TIJ(16),SIGTIJ(16,16) DIMENSION IQUADR(4) DIMENSION U(16) INTEGER POLDEG DOUBLE PRECISION DSIN,DCOS,DSQRT,DABS,DSIGN 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 SET-UP VALUES FOR THE DEGREE OF THE POLYNOMIAL 'POLDEG' (FOR THE FREE-AIR C ANOMALY AND THE HEIGHT), AND THE MAXIMUM DIMENSIONAL SPACE ALLOCATED C FOR EACH ARRAY 'MAXDIM', AND THE DESIRED INTEGRATION DISTANCE 'DIST' C (IN KILOMETERS) ALONG THE LEVELLING LINE, AND THE REASONABLE SCALE-FACTOR C 'SCALXY' TO BE MULTIPLIED BY THE LINEAR UNITS; SPECIFICALLY: THE X & Y C PLANE-COORDINATES, THE INTEGRATION LENGTH 'DIST' OF THE LEVELLING LINE, C AND THE MEAN RADIUS OF CURVATURE 'R0M' OF THE ELLIPSOID- IN ORDER TO GET C A WELL-CONDITIONED MATRIX-Q OF THE NORMAL EQUATIONS FOR THE INVERSION C PROCESS. C POLDEG=3 MAXDIM=2000 DIST=40.0D 0 SCALXY=0.03D 0 C C SET-UP THE NECESSARY VALUE OF THE SCALE-FACTOR 'SCALE', WHICH WILL BE C MULTIPLIED (IF REQUIRED) BY THE NORMAL-EQUATIONS MATRIX-Q, TO OVERCOME C THE OVERFLOW (OR UNDERFLOW) PROBLEM DURING THE INVERSION PROCESS. C SCALE=1.0D 0 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 DEGREE-SQUARE -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(2F8.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 DEGREE- 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 ONE 3-DEGREE-SQUARE CENTRE'//) WRITE(6,33) PHI,ALONG 33 FORMAT(10X,'THE CENTRE OF THE ONE-DEGREE-SQUARE UNDER INVESTIGATIO 1N IS LOCATED AT :'/,15X,'LATITUDE :',F6.1,2X,'DEGREES (NORTH)'/15 2X,'LONGITUDE :',F6.1,2X,'DEGREES (WEST)'//,40X,'================== 3============'//) 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-DEGREE-SQUARE 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 FYMAX=PHI+0.5D 0 FYMIN=PHI-0.5D 0 ALMAX=ALONG+0.5D 0 ALMIN=ALONG-0.5D 0 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 FOUR DEGREE-SQUARES 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.LT.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 'APPROX' 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***'/) CALL APPROX(NP,MPOLY,DGFA,WFA,SCALE,CIJ,SIGCIJ,EVFCIJ,U,DETQ) IF(DABS(DETQ).LT.1.0D-60) GO TO 960 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-DEGREE- C SQUARE UNDER INVESTIGATION. 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-DEGREE-SQUARE IS =',F10.4,' MGALS'/,10X,'WITH AN ESTIMAT 2ED STANDARD DEVIATION =',F10.4,' MGALS'//) C C******************************************************************************* C C CALL ALSO THE SAME SUBROUTINE 'APPROX' 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 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 APPROX(NP,MPOLY,HT,WHT,SCALE,TIJ,SIGTIJ,EVFTIJ,U,DETQ) IF(DABS(DETQ).LT.1.0D-60) GO TO 970 C C EXTRACT AND PRINT-OUT THE PREDICTED VALUE OF THE HEIGHT, ALONG WITH C ITS STANDARD DEVIATION, AT THE CENTRE OF THE ONE-DEGREE- SQUARE 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-D 1EGREE-SQUARE 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,DIST,SCALXY) 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-DEGREE-SQUARE AREA C SPECIFIED AROUND THE POINT OF INTEREST, FOR THE INTERPOLATION PURPOSES. C PLOTFY=0.5D 0*FYFACT*SCALXY PLOTAL=0.5D 0*ALFACT*SCALXY XFYMIN=-PLOTFY YALMIN=-PLOTAL XFYMAX=+PLOTFY YALMAX=+PLOTAL C CALL GCPLOT(YAL,XFY,NP,YALMIN,XFYMIN,YALMAX,XFYMAX) C CMFLSC=4.0D 0 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-DEGREE 1-SQUARE*****'/) WRITE(6,33) PHI,ALONG GO TO 9999 C**************R**************************************************************** C 960 WRITE(6,961) 961 FORMAT(///,10X,'ERROR MESSAGE : THE GRAMS MATRIX OF NORMAL EQUATIO 1NS (Q) FOR ESTIMATING THE FREE-AIR GRAVITY'/,15X,'ANOMALY SURFACE 2TURNED-OUT TO BE SINGULAR, SOLUTION IS IMPOSSIBLE'/,10X,'PLEASE CH 3ECK YOUR SCALING OF THE (X,Y) LOCAL COORDINATES OF DATA POINTS AND 4/OR'/15X,'THE SCALING OF THE MATRIX-Q ITSELF') GO TO 9999 970 WRITE(6,971) 971 FORMAT(///,10X,'ERROR MESSAGE : THE GRAMS MATRIX OF NORMAL EQUATIO 1NS (Q) FOR ESTIMATING THE HEIGHT (TOPOGRAPHY) '/,15X,'SURFACE 2TURNED-OUT TO BE SINGULAR, SOLUTION IS IMPOSSIBLE'/,10X,'PLEASE CH 3ECK YOUR SCALING OF THE (X,Y) LOCAL COORDINATES OF DATA POINTS AND 4/OR'/15X,'THE SCALING OF THE MATRIX-Q ITSELF') GO TO 9999 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 UNKNOWN COEFFI 2CIENTS (-VE DEGREES OF FREEDOM)') 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') 777 FORMAT('1') C 9999 GO TO 999 C C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C 99 STOP END C C C*************************************************************************** C* * C* SUBROUTINE 'A P P R O X' :- (BY: M.M. NASSAR, JANUARY 1975.) * C* ------------------------ * C* * C* THIS SUBROUTINE HAS BEEN DESIGNED AS A * C* GENERAL SUBROUTINE FOR PERFORMING A LEAST-SQUARES APPROXIMATION TO ANY * C* DATA SET, USING A MIXED-ALGEBRAIC POLYNOMIAL OF ANY ARBITRARY ORDER * C* IN TWO DIMENSIONS, WHICH IMPLIES THAT THIS SUBROUTINE HAS THE ADVANTAGE * C* OF EXECUTION-TIME-DIMENSIONING. * 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 INDIVIDUAL ELEMENTS OF THE DESIGN MATRIX-A ARE GENERATED (VIA THE * C* FUNCTION-SUBPROGRAM 'A') ONE AT A TIME WHENEVER NEEDED IN THE COMPUTAT- * C* IONS WITHOUT THE NECESSITY OF STORING THEM, TO SAVE COMPUTER STORAGE. * C* (SEE THE LISTINGS OF THE FUNCTION-SUBPROGRAM 'A' CONCERNING THE REQUIRED* C* COMMON BLOCK WITH THE MAIN PROGRAM). * C* * C* THE INVERSION OF THE NORMAL EQUATIONS MATRIX IS DONE VIA THE SUBROUTINE * C* 'CHOLD', WHICH IS BASED ON THE CHOLESKY DECOMPOSITION METHOD. IN THIS * C* RESPECT, APPROX CONTAINS THE POSSIBILITY OF SCALING THE NORMAL EQUATIONS* C* MATRIX (AS A WHOLE) BY THE DESIRED SCALE FACTOR, IF REQUIRED BEFORE * C* AND AFTER THE INVERSION PROCESS. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* USAGE:- CALL APPROX(N,M,F,W,SCALE,C,Q,EVF,U,DETQ) * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* THE REQUIRED INPUT DATA :- (VIA THE CALLING ROUTINE) * C* ----------------------- * C C* N : TOTAL NUMBER OF AVAILABLE DATA POINTS. * C* M : NUMBER OF UNKNOWN COEFFICIENTS OF THE APPROXIMATING POLYNOMIAL. * C* F : VECTOR OF OBSERVED DATA AT THE N DISCRETE POINTS. * C* W : VECTOR OF WEIGHTS FOR THE N OBSEVED QUANTITIES. * C* SCALE : DESIRED SCALE FACTOR FOR THE NORMAL EQUATIONS MATRIX. * C* * C* IT SHOULD BE NOTED HERE THAT THE VECTORS OF X & Y LOCAL COORDINATES , * C* AS WELL AS THE SPECIFIED ORDER OF THE POLYNOMIAL ARE PART OF THE INPUT * C* DATA ,HOWEVER THEY MUST BE DECLARED AS A COMMON BLOCK BETWEEN THE MAIN * C* PROGRAM AND THE FUNCTION 'A'. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* THE MAIN OUTPUT ARE :- * C* ------------------- * C* * C* C : VECTOR OF ESTIMATED COEFFICIENTS OF THE APPROXIMATING POLYNOMIAL * C* Q : ESTIMATED VARIANCE-COVARIANCE MATRIX OF THE COEFFICIENTS 'C'. * C* EVF : THE A POSTERIORI VARIANCE FACTOR. * C* U : VECTOR OF ABSOLUTE TERMS OF THE NORMAL EQUATIONS. * C* DETQ : DETERMINANT OF THE NORMAL EQUATIONS MATRIX. * C* * C* NOTE HERE THAT THE REASON FOR INCLUDING THE VECTOR-U IN THE ARGUMENT * C* LIST IS SOLELY TO ACHIEVE THE ADVANTAGE OF ITS EXECUTION TIME DIMENSIONS* C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* ON THE OTHER HAND THE PRINTED OUTPUT FROM THE LEAST-SQUARES APPROXIMAT- * C* ION TECHNIQUE ,BY THIS SUBROUTINE, INCLUDES :- * C* * C* 1-THE NORMAL EQUATIONS MATRIX. * C* 2-THE VECTOR-U OF ABSOLUTE TERMS OF THE NORMAL EQUATIONS. * C* 3-THE DETERMINANT OF THE NORMAL EQUATIONS MATRIX. * C* 4-THE VECTOR OF ESTIMATED COEFFICIENTS OF THE APPROXIMATING POLYNOMIAL * C* 5-THE VECTOR-V OF RESIDUALS. * C* 6-THE NUMBER OF DEGREES OF FREEDOM. * C* 7-THE ESTIMATED SQUARE-SUM OF THE WEIGHTED RESIDUALS. * C* 8-THE A POSTERIORI (ESTIMATED) VARIANCE FACTOR. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* FINALLY, IT SHOULD BE NOTED THAT THE ORDER OF THE APPROXIMATING POLYNOM-* C* IAL IS LEFT ARBITRARY, HOWEVER, THE NUMBER OF DATA POINTS IS ALLOWED * C* ONLY UP TO 2000, AND IF EXCEEDED THE CORRESPONDING DIMENSION STATEMENTS * C* IN THE MAIN PROGRAM AND THE FUNCTION 'A' MUST BE MODIFIED. * C* * C*************************************************************************** C SUBROUTINE APPROX(N,M,F,W,SCALE,C,Q,EVF,U,DETQ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION F(N),W(N) DIMENSION Q(M,M),U(M),C(M) C C INITIALIZE ZERO-VALUES FOR THE ELEMENTS OF SOME RESULTANT VECTORS & MATRICES C DO 15 I=1,M DO 15 J=1,M 15 Q(I,J)=0.0D 0 DO 25 I=1,M U(I)=0.0D 0 25 C(I)=0.0D 0 C C******************************************************************************* 888 FORMAT('1') C******************************************************************************* C C FORMULATE THE NORMAL-EQUATIONS MATRIX (WHICH IS KNOWN AS 'GRAM'S-MATRIX C IN THE LEAST-SQUARES TERMINOLOGY, AND WHOSE ELEMENTS ARE THE COMBINATIONS C OF THE SCALAR-PRODUCT OF THE BASE-FUNCTIONS). THIS MATRIX 'Q' IS GIVEN BY: C Q = AT * P * A , C WHERE 'A' IS THE DESIGN MATRIX, AND 'P' IS THE WEIGHT MATRIX. 'A' IS KNOWN C AS VANDERMONDE'S-MATRIX IN THE LEAST-SQUARES APPROXIMATION TERMINOLOGY, C THE ELEMENTS IN EACH ROW OF 'A' PERTAIN TO ONLY ONE OBSERVED DATA POINT, C AND ARE SIMPLY THE FUNCTIONAL VALUES OF THE M CHOSEN BASE-FUNCTIONS IN C THE MIXED (ALGEBRAIC, NON-LINEAR, THIRD-ORDER, IN TWO-DIMENSIONS) C APPROXIMATING POLYNOMIAL. I.E. ONE ROW IN A-MATRIX FOR EACH DATA POINT C C IT SHOULD BE NOTED HERE THAT THE INDIVIDUAL ELEMENTS 'A(ROW,COL.)' OF THE C A-MATRIX WHICH ARE NEEDED FOR COMPUTING THE Q-MATRIX, WILL BE EVALUATED C ONE AT A TIME THROUGH THE FUNCTION-SUBPROGRAM ''A'' WHICH IS AUTOMATICALLY C CALLED AS SOON AS THE EXECUTION HITS ANY EXPRESSION INVOLVING ANY ELEMENT C OF THE DESIGN MATRIX-A. C C THE 'P' MATRIX IS A DIAGONAL WEIGHT MATRIX, CONTAINING THE WEIGHT FUNCTIONS C W'S ON ITS MAIN DIAGONAL. C DO 20 I=1,M DO 20 J=1,M Q(I,J)=0.0D 0 DO 20 K=1,N Q(I,J)=Q(I,J)+A(K,I)*A(K,J)*W(K) 20 CONTINUE C C PRINT-OUT THE NORMAL-EQUATIONS MATRIX C WRITE(6,81) M,M 81 FORMAT(/,5X,'THE NORMAL EQUATIONS MATRIX-Q'/,10X,'DIMENSION : (',I 13,' ,',I3,' )') WRITE(6,82) 82 FORMAT(/,10X,'1',12X,'2',12X,'3',12X,'4',12X,'5',12X,'6',12X,'7',1 12X,'8'/) MHALF=M/2 MHALF1=MHALF+1 DO 110 I=1,M 110 WRITE(6,83) I,(Q(I,J),J=1,MHALF) 83 FORMAT(I4,8D13.4) WRITE(6,84) 84 FORMAT(/,10X,'9',11X,'10',11X,'11',11X,'12',11X,'13',11X,'14',11X, 1'15',11X,'16'/) DO 120 I=1,M 120 WRITE(6,83) I,(Q(I,J),J=MHALF1,M) C******************************************************************************* C C FORMULATE THE VECTOR-U OF ABSOLUTE TERMS OF THE NORMAL EQUATIONS (WHOSE C ELEMENTS ARE THE SCALAR PRODUCTS OF THE BASE-FUNCTIONS AND THE OBSERVED- C FUNCTION (F-VECTOR), I.E. U=AT*P*F WHERE A & P ARE AS DEFINED C ABOVE, AND 'F' IS THE VECTOR OF OBSERVED QUANTITIES -TO BE USED IN COMPUTING C THE BEST-FITTING POLYNOMIAL. C DO 30 I=1,M U(I)=0.0D 0 DO 30 J=1,N U(I)=U(I)+A(J,I)*W(J)*F(J) 30 CONTINUE C C PRINT-OUT THE VECTOR 'U' C WRITE(6,888) WRITE(6,85) M 85 FORMAT(//,5X,'THE VECTOR-U OF ABSOLUTE TERMS OF THE NORMAL EQUATIO 1NS'/,10X,'DIMENSIONS : (',I3,' , 1 )') WRITE(6,86) 86 FORMAT(/,12X,'1'/) DO 130 I=1,M 130 WRITE(6,87) I,U(I) 87 FORMAT(I4,D15.5) C******************************************************************************* IF(SCALE.EQ.1.0D 0) GO TO 215 C C SCALE THE NORMAL-EQUATIONS MATRIX-Q, BEFORE INVERSION, TO OVERCOME THE C OVERFLOW PROBLEM WHICH IS LIKELY TO OCCUR DURING THE INVERSION PROCESS. C DO 5 I=1,M DO 5 J=1,M 5 Q(I,J)=Q(I,J)*SCALE 215 CONTINUE C C CALL THE SUBROUTINE 'CHOLD' TO INVERT THE NORMAL-EQUATIONS MATRIX C 'Q' (I.E. GRAM'S-MATRIX), UNDER THE SAME NAME 'Q' AND IN THE SAME LOCATION C THIS INVERSE CAN BE REFERRED TO AS THE WEIGHT-COEFFICIENT MATRIX (IN THE C LEAST-SQUARES SENSE) OF THE COEFFICIENTS OF THE BEST-FITTING POLYNOMIAL. C ***FOR DETAILS REGARDING THE DESCRIPTION OF THE INPUT/OUTPUT DATA FOR C 'CHOLD', SEE ITS DOCUMENTATION.*** C CALL CHOLD(Q,M,M,DETQ,&97) IF(SCALE.EQ.1.0D 0) GO TO 225 C C RE-SCALE THE INVERSE MATIX-Q (OF THE NORMAL EQUATIONS), AFTER THE INVERSION C PROCESS, IN ORDER TO GET THE TRUE (CORRECT) WEIGHT-COEFFICIENT MATRIX C OF THE ESTIMATED POLYNOMIAL-COEFFICIENTS. C DO 210 I=1,M DO 210 J=1,M 210 Q(I,J)=Q(I,J)*SCALE 225 CONTINUE WRITE(6,372) DETQ 372 FORMAT(/10X,'DETERMINANT OF NORMAL-EQUATIONS MATRIX-Q =',D15.5/) C******************************************************************************* C C COMPUTE THE ESTIMATED VALUES OF THE COEFFICIENTS (C'S) OF THE BEST-FITTING C POLYNOMIAL. C=Q(INVERSE)*U C DO 40 I=1,M C(I)=0.0D 0 DO 40 J=1,M C(I)=C(I)+Q(I,J)*U(J) 40 CONTINUE C C PRINT-OUT THE VECTOR-C OF ESTIMATED COEFFICIENTS C WRITE(6,88) M 88 FORMAT(//,5X,'THE VECTOR-C OF ESTIMATED POLYNOMIAL-COEFFICIENTS'/, 110X,'DIMENSIONS : (',I3,' , 1 )') WRITE(6,86) DO 140 I=1,M 140 WRITE(6,87) I,C(I) C******************************************************************************* C C COMPUTE THE SQUARE-SUM OF THE WEIGHTED RESIDUALS VTPV=VT*P*V, WHERE P C IS AS DEFINED EARLIER, AND 'V' IS THE VECTOR OF RESIDUALS (OR DISCREPENCIES) C WHICH ARE GIVEN BY THE DIFFERENCE BETWEEN THE OBSERVED VALUES (VECTOR-F) C AND THE CORRESPONDING PREDICTED VALUES 'A*C' (AS COMPUTED FROM THE BEST- C FITTING POLYNOMIAL), I.E. V=A*C-F, IT SHOULD BE NOTED HERE THAT C THESE RESIDUALS HAVE THE SAME PHYSICAL UNITS AS THE OBSERVED VALUES 'F' C C IT SHOULD BE REMEMBERED HERE AGAIN THAT THE INDIVIDUAL ELEMENTS OF THE C DESIGN MATRIX-A ARE EVALUATED VIA THE FUNCTION-SUBPROGRAM ''A'' C VTPV=0.0D 0 WRITE(6,888) WRITE(6,122) N 122 FORMAT(//,5X,'THE VECTOR-V OF RESIDUALS'/,5X,'DIMENSIONS : (',I4,' 1 , 1 )') WRITE(6,86) DO 60 I=1,N AC=0.0D 0 DO 50 J=1,M 50 AC=AC+A(I,J)*C(J) C C COMPUTE AND PRINT-OUT THE VECTOR OF RESIDUALS -DISCREPANCIES BETWEEN C THE OBSERVED QUANTITIES AND THEIR CORRESPONDING PREDICTED VALUES. C V=F(I)-AC WRITE(6,87) I,V VTPV=VTPV+V*W(I)*V 60 CONTINUE C C COMPUTE THE NUMBER OF DEGREES OF FREEDOM, WHICH IS THE DIFFERENCE BETWEEN C THE TOTAL NUMBER OF USED OBSERVATION POINTS AND THE NUMBER OF UNKNOWN C POLYNOMIAL-COEFFICIENTS C NUMDF=N-M C C COMPUTE THE ESTIMATED (IN THE LEAST-SQUARES SENCE) VARIANCE FACTOR 'EVF' C -DENOTED BY SIGMA-ZERO-SQUARE-CAP. IT SHOULD BE NOTED HERE, THAT THE EVF C IS UNITLESS, DUE TO OUR STIPULATION OF THE WEIGHTS (INVERSE OF VARIANCES) C EVF=VTPV/FLOAT(NUMDF) C C COMPUTE THE ESTIMATED VARIANCE-COVARIANCE MATRIX 'EC' (SIGMA-C) OF THE C ESTIMATED POLYNOMIAL-COEFFICIENTS (C'S), SIMPLY BY MULTIPLYING THE C ESTIMATED VARIANCE-FACTOR BY THE INVERSE OF THE NOMAL-EQUATIONS MATRIX C I.E. EC=EVF*Q(INVERSE) C DO 70 I=1,M DO 70 J=1,M Q(I,J)=EVF*Q(I,J) 70 CONTINUE C C PRINT-OUT THE 'NUMDF' , 'VTPV' AND 'EVF' SCALARS, AND THE 'EC' MATRIX C WRITE(6,89) NUMDF,VTPV,EVF 89 FORMAT(/,/,5X,'THE NUMBER OF DEGREES OF FREEDOM =',I5//,5X,'THE ES 1TIMATED SQUARE-SUM OF THE WEIGHTED RESIDUALS =',D15.5//,5X,'THE ES 2TIMATED VARIANCE FACTOR =',D15.5/) WRITE(6,888) WRITE(6,91) M,M 91 FORMAT(/,5X,'THE ESTIMATED VARIANCE-COVARIANCE MATRIX-EC OF THE ES 1TIMATED POLYNOMIAL-COEFFICIENTS'/,15X,'DIMENSIONS : (',I3,' ,',I3, 2' )') WRITE(6,82) DO 150 I=1,M 150 WRITE(6,83) I,(Q(I,J),J=1,MHALF) WRITE(6,84) DO 160 I=1,M 160 WRITE(6,83) I,(Q(I,J),J=MHALF1,M) 97 RETURN END C C*************************************************************************** C* * C* FUNCTION-SUBPROGRAM ' A ' :- * C* ------------------------- * C* * C* THIS FUNCTION-SUBPROGRAM IS DESIGNED AS A * C* GENERAL FUNCTION, TO COMPUTE (IN DOUBLE PRECISION) SINGLE ELEMENTS * C* A(I,J) OF THE DESIGN MATRIX-A (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 (IDEG), WHERE I IS THE ROW * C* POSITION AND J IS THE COLUMN POSITION OF THE SOUGHT ELEMENT A(I,J). * 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-A 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 A(I,J)) * C* AS SOON AS THE EXECUTION IN THE CALLING ROUTINE HITS ANY * C* EXPRESSION THAT INVOLVES ANY INDIVIDUAL ELEMENT A(I,J) OF A. * 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 A(I,J) 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 2000 , OTHERWISE THE CORRESPONDING DIMENSION * C* STATEMENTS SHOULD BE MODIFIED,.. * C* * C*************************************************************************** C DOUBLE PRECISION FUNCTION A(I,J) DOUBLE PRECISION X,Y,XA,YA DIMENSION X(2000),Y(2000) C C DECLARE EXPLICITLY THE LABELLED COMMON BLOCK BETWEEN THE MAIN-PROGRAM C AND THIS FUNCTION-SUBPROGRAM C THIS COMMON BLOCK SHOULD CONTAIN THE VECTORS OF X & Y PLANE LOCAL C COORDINATES, AND THE DEGREE OF THE APPROXIMATING MIXED-ALGEBRAIC POLYNOMIAL C COMMON /XYMAIN/ X,Y,IDEG N=IDEG+1 K=(J-1)/N L=J-K-(N-1)*K-1 XA=X(I) YA=Y(I) C C CHECK AND RESPOND TO OVERCOME THE PROBLEM OF TERMINATING THE PROGRAM DUE C TO THE OCCURENCE OF A BASE EQUALS ZERO WITH AN EXPONENT ALSO EQUALS ZERO, C WHICH WILL RESULT IN AN ERROR-MESSAGE AND STOPPING THE PROGRAM-EXECUTION C IF(XA.EQ.0.0D 0.AND.K.EQ.0) XA=1.0D 0 IF(YA.EQ.0.0D 0.AND.L.EQ.0) YA=1.0D 0 A=XA**K*YA**L RETURN END C C C*************************************************************************** C* * C* SUBROUTINE 'CHOLD' : * C* ------------------ * C* THE PURPOSE OF THIS SUBROUTINE IS TO INVERT A * C* POSITIVE DEFINITE SYMMETRIC MATRIX - USING THE 'CHOLESKY' * C* DECOMPOSITION METHOD. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* THE SUBROUTINE WAS WRITTEN BY 'D.E. WELLS', IN 1973, (REFERENCE : * C* 'MATRICES',BY D.E. WELLS - UNB DEPT. OF SURVEYING ENG. LECTURE NOTES * C* # 15, THE SECOND PRINTING, AUGUST 1973). * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* USAGE : CALL CHOLD(A,IRDA,NA,DETA,&INT) * C* * C* REQUIRED INPUT DATA ARE :- (THROUGH THE CALLING ROUTINE) * C* ----------------------- * C* A : IS THE MATRIX WHOSE ROW DIMENSION IS (IRDA), AND CONTAINING * C* THE POSITIVE DEFINITE SYMMETRIC MATRIX (NA BY NA)-REQUIRED * C* TO BE INVERTED- AS ITS UPPER LEFT MOST SUBMATRIX, * C* IRDA : IS THE ROW DIMENSION OF THE MATRIX A, * C* NA : THE SIZE OF THE UPPER-LEFT MOST SUBMATRIX OF A, WHICH IS * C* REQUIRED TO BE INVERTED, * C* INT : IS AN INTEGER SPECIFYING THE STAEMENT NUMBER IN THE MAIN * C* CALLING ROUTINE, TO WHICH THE CONTROL SHOULD BE TRANSFERRED * C* IF THE MATRIX TO BE INVERTED IS EITHER SINGULAR (DETA<1D-10) * C* OR HAS NO SIZE (NA<1). * C* * C* OUTPUT DATA ARE :- (TO THE MAIN CALLING ROUTINE) * C* --------------- * C* A : THE ORIGINAL MATRIX (A) WHOSE UPPER-LEFT MOST SUBMATRIX * C* (NA BY NA) IS REPLACED BY ITS REQUIRED INVERSE (I.E. THE * C* ORIGINAL (NA BY NA) SUBMATRIX WILL BE DESTROYED), * C* DETA : IS THE DETERMINANT OF THE (NA BY NA) SUBMATRIX TO BE INVERTED. * C* * C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* C* * C* NOTE THAT NORMALLY, WE HAVE THE ORIGINAL MATRIX (A) IS THE SAME ONE * C* WANTED TO BE INVERTED, IN SUCH CASE THE ARGUMENTS 'NA' AND 'IRDA' * C* IN THE CALLING STATEMENT WILL BE EQUAL. * C* * C* NOTE ALSO, THAT THERE IS NO LIMITATION ON THE SIZE OF THE MATRIX TO * C* BE INVERTED. HOWEVER, THIS WILL BE RESTRICTED BY THE AVAILABLE * C* COMPUTER TIME. * C* * C*************************************************************************** C SUBROUTINE CHOLD(A,IRDA,NA,DETA,*) DOUBLE PRECISION A,DETA,SUM,SQRT,DSQRT,ABS,DABS,SING DIMENSION A(IRDA,NA) SQRT(SUM) = DSQRT(SUM) ABS(DETA) = DABS(DETA) DATA SING/1D-60/ C C CHOLESKI DECOMPOSITION OF INPUT MATRIX INTO TRIANGULAR MATRIX C IF(NA .LT. 1) GO TO 18 DETA = A(1,1) A(1,1) = SQRT(A(1,1)) IF(NA .EQ. 1) GO TO 6 DO 1 I = 2,NA 1 A(I,1) = A(I,1) / A(1,1) DO 5 J = 2,NA SUM = 0. J1 = J - 1 DO 2 K = 1,J1 2 SUM = SUM + A(J,K) ** 2 DETA = DETA * (A(J,J) - SUM) A(J,J) = SQRT(A(J,J) - SUM) IF(J .EQ. NA) GO TO 5 J2 = J + 1 DO 4 I = J2,NA SUM = 0. DO 3 K = 1,J1 3 SUM = SUM + A(I,K) * A(J,K) 4 A(I,J) = (A(I,J) - SUM) / A(J,J) 5 CONTINUE 6 IF(ABS(DETA) .LT. SING) GO TO 16 C C INVERSION OF LOWER TRIANGULAR MATRIX C DO 7 I = 1,NA 7 A(I,I) = 1. / A(I,I) IF(NA .EQ. 1) GO TO 10 N1 = NA - 1 DO 9 J = 1,N1 J2 = J + 1 DO 9 I = J2,NA SUM = 0. I1 = I - 1 DO 8 K = J,I1 8 SUM = SUM + A(I,K) * A(K,J) 9 A(I,J) = - A(I,I) * SUM C C CONSTRUCTION OF INVERSE OF INPUT MATRIX C 10 DO 15 J = 1,NA IF(J .EQ. 1) GO TO 12 J1 = J - 1 DO 11 I = 1,J1 11 A(I,J) = A(J,I) 12 DO 14 I = J,NA SUM = 0. DO 13 K = I,NA 13 SUM = SUM + A(K,I) * A(K,J) 14 A(I,J) = SUM 15 CONTINUE RETURN 16 WRITE(6,17) DETA 17 FORMAT(10X, 'SINGULAR MATRIX IN CHOLD. DET =',E20.5) RETURN 1 18 WRITE(6,19) 19 FORMAT(10X,'MATRIX OF DIMENSION ZERO IN CHOLD') RETURN 1 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 (40 KM. LONG), RADIATING FROM THE CENTRE OF * C* A ONE-DEGREE-SQUARE , 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) * 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-DEGREE-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-DEGREE- * C* SQUARE, RESPECTIVELY. * C* S : THE DESIRED INTEGRATION DISTANCE (IN KM.) ALONG THE SIMULAT-* C* ED LEVELLING LINE (E.G., 40 KM.). * 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* * 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-DEGREE- SQUARE AT A TIME. * C* * C*************************************************************************** C SUBROUTINE GCAFAZ(CIJ,TIJ,PHI0,R0M,S,SC) IMPLICIT REAL*8(A-H,O-Z) DIMENSION CIJ(16),TIJ(16),C(5),T(5),ODASH(4),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=10.905D 0 A1=13.244D 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-DEGREE- 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-DEGREE- 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 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-DEGREE-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