C DATA SET PREDICT AT LEVEL 011 AS OF 09/30/80 C ******************************************************************* 00001**9 C * * 00002**9 C * * 00003**9 C * * 00004**9 C * PROGRAM 'P R E D I C T' : * 00005**9 C * ========================= * 00006**9 C * * 00007**9 C * * 00008**9 C * * 00009**9 C * * 00010**9 C * THIS COMPUTER PROGRAM PACKAGE HAS BEEN DESIGNED TO PREDIC * 00011**9 C * T FOR GRAVITY VALUES AT THE CENTRE OF A ONE DEGREE SQUARE BLOC * 00012**9 C * K USING (1) LEAST-SQUARES SURFACE FIT TECHNIQUE, (2) METHOD OF * 00013**9 C * WEIGHTED MEAN, AND (3) COMBINATION OF SURFACE FIT AND COLLOCAT * 00014**9 C * ION METHODS. IT IS A MODIFIED VERSION OF PART OF PROGRAM * 00015**9 C * 'LEVAGRAV' DEVELOPED BY M.M. NASSAR IN 1975. * 00016**9 C * * 00017**9 C * FOR DETAILS CONCERNING THE MATHEMATICAL MODELS, SEE: * 00018**9 C * KASSIM (1980). ''AN EVALUATION OF THREE TECHNIQUES FOR THE PR- * 00019**9 C * EDICTION OF GRAVITY ANOMALIES IN CANADA''. M.S.C THESIS, DEPT. * 00020**9 C * OF SURV. ENG.,U. N.B., FREDERICTON. * 00021**9 C * * 00022**9 C * THE MAIN PROGRAM 'PREDICT' READS THE LATITUDE AND LONGIT- * 00023**9 C * UDE OF THE CENTRE OF THE ONE DEGREE SQ. BLOCK, SEEKS THE APPR- * 00024**9 C * OPRIATE DATA FILE AND STORES ALL AVAILABLE DATA FOR THE BLOCK * 00025**9 C * COMPUTING AT THE SAME TIME THE X,Y COORDINATES OF EACH DATA * 00026**9 C * POINT WITH THE PREDICTION POINT AS ORIGIN. THE CLOSEST DATA TO * 00027**9 C * THE PREDICTION POINT ARE LATER SELECTED FOR THE PREDICTION. * 00028**9 C * * 00029**9 C * 'PREDICT' IS SET UP SUCH THAT ONLY ONE TECHNIQUE CAN BE * 00030**9 C * USED AT ATIME. PROVISION IS MADE TO HAVE DATA SCREENED FOR * 00031**9 C * CLUSTERS OF POINTS WHICH ARE REPLACED BY A WEIGHTED MEAN. * 00032**9 C * * 00033**9 C * PREDICTION BY THE SURFACE FIT METHOD IS DONE VIA SUBROUTINE * 00034**9 C * 'APPROX' * 00035**9 C * * 00036**9 C * PREDICTION BY THE WEIGHTED MEAN METHOD IS DONE VIA MAIN PROGR * 00037**9 C * AM * 00038**9 C * * 00039**9 C * PREDICTION BY THE COMBINATION METHOD IS DONE VIA SUBROUTINES * 00040**9 C * 'APPROX' AND 'GRAV'. THE COVARIANCE FUNCTION FOR GRAV MAY BE * 00041**9 C * OBTAINED BY LINKING THIS PACKAGE WITH THE LIBRARY PACKAGE * 00042**9 C * CONTAINING SUBROUTINE 'COVAX'; OR ALTERNATILY, A PLYNOMIAL * 00043**9 C * APPROXIMATION TO THE ANALYTIC FUNCTION MAY BE USED. * 00044**9 C * * 00045**9 C * SINCE THE PROGRAM IS DESIGNED TO PREDICT FOR ONE DEGREE * 00046**9 C * BLOCK AT A TIME, IT CAN PROCESS AN UNLIMITED NUMBER OF BLOCKS * 00047**9 C * PROVIDED ENOUGH TIME IS SPECIFIED. * 00048**9 C * * 00049**9 C * PROVISION IS MADE FOR RESULTS OF A SERIES OF PREDICTIONS * 00050**9 C * TO BE TESTED STATISTICALLY. THE TESTS INCLUDE CHI-SQUARE GOOD- * 00051**9 C * NESS OF FIT TEST, MEASURES OF SKEWNESS AND KURTOSIS, AND TESTS * 00052**9 C * ON MEANS AND VARIANCES. THE HITOGRAM OF STANDARDIZED OBSERVA- * 00053**9 C * TIONS IS ALSO PLOTTED. * 00054**9 C * * 00055**9 C * REQUIRED INPUT DATA: * 00056**9 C * -------------------- * 00057**9 C * * 00058**9 C * 1. NAME DESCRIBING TERRAIN TYPE. IT IS READ IN AS A CHARACTER * 00059**9 C * STRING. * 00060**9 C * FORMAT(25A1) * 00061**9 C * * 00062**9 C * 2. PARAMETERS DEFINING: * 00063**9 C * (I) GEGREE OF POLYNOMIAL, * 00064**9 C * (II) EXTENT OF BLOCK FOR PREDICTION (1 DEGREE, SAY), * 00065**9 C * (III) RADIAL INCREMENT FOR SELECTION OF DATA CLOSEST TO * 00066**9 C * PREDICTION POINT, * 00067**9 C * (IV) SCALE FOR X,Y COORDINATES, * 00068**9 C * (V) MAXIMUM NUMBER OF DATA TO BE USED FOR PREDICTION, AND * 00069**9 C * (VI) NUMBER FO POINTS IN POLYGON (INPUT INFORMATION FOR * 00070**9 C * SUBROUTINE 'INPOLY' USED FOR SELECTION OF DATA SURROUNDING * 00071**9 C * PREDICTION POINT) * 00072**9 C * FORMAT(I1,3D10.2,2I3) * 00073**9 C * * 00074**9 C * 3. PARAMETERS SPECIFYING: * 00075**9 C * (I) WHETHER (OR NOT) TO PRINT INTERMEDIATE RESULTS, * 00076**9 C * (II) WHETHER (OR NOT) TO USE WEIGHTED MEAN TECHNIQUE, * 00077**9 C * (III) THE VALUE OF POWER FOR WEIGHTING OBSERVATIONS IN THE * 00078**9 C * WEIGHTED MEAN TECHNIQUE, * 00079**9 C * (IV) WHETHER (OR NOT) TO SCREEN DATA FOR CLUSTERS OF POINTS * 00080**9 C * (V) WHETHER (OR NOT) TO USE REGRESSED DDTA. * 00081**9 C * PARAMETERS (I),(II),(IV) AND (V) CAN ONLY HAVE VALUES OF * 00082**9 C * 1 OR 0 CORRESPONDING TO THE ANSWERS YES OR NO. * 00083**9 C * FORMAT(5I2) * 00084**9 C * * 00085**9 C * 4. PARAMETERS SPECIFYING: * 00086**9 C * (I) LIMIT OF SEPARATION BETWEEN PAIRS OF DATA FOR SCREENING * 00087**9 C * (IF SEPARATION IS LESS THAN LIMIT, THE DATA ARE CONSIDERED * 00088**9 C * TOO CLOSE AND ARE REPLACED BY A WEIGHTED MEAN), * 00089**9 C * (II) WHETHER (OR NOT) TO COMBINE SURFACE FIT AND COLLOCA- * 00090**9 C * TION METHODS, * 00091**9 C * (III) WHETHER (OR NOT) TO COMBINE SURFACE FIT AND WEIGHTED * 00092**9 C * MEAN METHODS * 00093**9 C * PARAMETERS (II) AND (III) CAN ONLY HAVE VALUES OF 1 OR 0. * 00094**9 C * FORMAT(D10.2,2I5) * 00095**9 C * * 00096**9 C * 5. REGRESSION PARAMETERS: * 00097**9 C * (I) SLOPE OF REGRESSION LINE, AND * 00098**9 C * (II) INTERCEPT OF REGRESSION LINE. * 00099**9 C * PARAMETERS APPLY ONLY FOR MOUNTAINOUS TERRAIN DATA * 00100**9 C * BUT, THEY MUST BE GIVEN SOME VALUES FOR THE PROGRAM TO WORK * 00101**9 C * FORMAT(2D14.7) * 00102**9 C * * 00103**9 C * 6. PARAMETERS SPECIFYING: * 00104**9 C * (I) WHETHER (OR NOT) DATA FOR COMPUTATION OF COVARIANCE * 00105**9 C * FUNCTION (FOR MOUNTAINOUS TERRAIN) ARE REQUIRED, * 00106**9 C * (II) WHETHER (OR NOT) THE DATA ARE TO BE STORED TEMPORARILY * 00107**9 C * ON A FILE, AND * 00108**9 C * (III) FILE NUMBER IF THE DATA ARE ATO BE STORED. * 00109**9 C * FORMAT(2I2,I4) * 00110**9 C * * 00111**9 C * 7. VARIABLES FOR CHI-SQUARE TEST: * 00112**9 C * (I) NUMBER OF CLASSES/SEGMENTS, * 00113**9 C * (II) CHI-SQUARE VALUE FROM TABLES. * 00114**9 C * FORNAT(I2,D10.3) * 00115**9 C * * 00116**9 C * 8. FACTOR FOR WITHIN COTEXT TEST OF OUTLIERS. * 00117**9 C * FORMAT(D12.4) * 00118**9 C * * 00119**9 C * 9. SCALE FACTOR FOR VARIANCE OF PPEDICTED ANOMALY TO CHECK * 00120**9 C * OPTIMISM OR PESSIMISM IN COMPUTING THE VARIANCES. * 00121**9 C * FORMAT(D12.5) * 00122**9 C * * 00123**9 C * 10. (I) IDENTIFICATION NUMBER OF PREDICTION POINT * 00124**9 C * (II) POSITION COORDINATES OF PREDICTION POINT: LATITUDE(+VE * 00125**9 C * N), LONGITUDE(+VE W) IN DECIMAL DEGREE, * 00126**9 C * (III) OBSERVED FREE-AIR ANOMALY (MGAL), * 00127**9 C * (IV)OBSERVED BOUGUER ANOMALY (MGAL), * 00128**9 C * (V) HEIGHT (METRES), * 00129**9 C * (VI) VARIANCE OF FREE-AIR ANOMALY (MGAL**2), * 00130**9 C * (VII) VARIANCE OF BOUGUER ANOMALY (MGAL**2), AND * 00131**9 C * (VIII) VARIANCE OF HEIGHT (METRE**2). * 00132**9 C * FORMAT(I5,5F10.5,2F6.2,F8.2) * 00133**9 C * * 00134**9 C * 11. GRAVITY DATA FILES SPECIFIED BY THEIR APPROPRIATE JCL CARDS* 00135**9 C * * 00136**9 C * THE DATA USED IN THIS PROGRAM ARE RESIDING ON SEVEN DISCS, THE * 00137**9 C * DETAILS OF WHICH WOULD BE FOUND IN NASSAR AND VANICEK (1975) * 00138**9 C * ''LEVELLING AND GRAVITY'',LECTURE NOTE NO. 33, DEPT OF S.E., * 00139**9 C * U.N.B., FREDERICTON. * 00140**9 C * * 00141**9 C * * 00142**9 C * OUTPUT: * 00143**9 C * ------ * 00144**9 C * * 00145**9 C * 1. (I) IDENTIFICATION OF PREDICTION POINT, * 00146**9 C * (II) POSITION COORDINATES - LATITUDE(+VE N), LONGITUDE(+VE * 00147**9 C * W) IN DECIMAL DEGREE, * 00148**9 C * (III) OBSERVED ANOMALY (MGAL) * 00149**9 C * (IV) PREDICTED ANOMALY (MGAL), * 00150**9 C * (V) STANDARD DEVIATION OF OBSERVED ANOMALY (MGAL) * 00151**9 C * (VI) STANDARD DEVIATION OF PREDICTED ANOMALY (MGAL), * 00152**9 C * (VII) DIFFERENCE BETWEEN OBSERVED AND PREDICTED ANOMALY * 00153**9 C * (VIII) RATIO OF DIFFERENCE TO QUADRATIC SUM OF STD. DEVIA- * 00154**9 C * TIONS OF OBSERVED AND PREDICTED ANOMALIES, AND * 00155**9 C * (IX) NUMBER OF DATA USED FOR THE PREDICTION. * 00156**9 C * * 00157**9 C * 2. AT THE END OF A SERIES OF PREDICTIONS: * 00158**9 C * (I) MEAN OF DIFFERENCES .... MD * 00159**9 C * (II) STANDARD DEVIATION OF THE MEAN .... S * 00160**9 C * * 00161**9 C * 3. RESULTS OF: * 00162**9 C * (I) TEST OF OUTLIERS (PERHAPS WITH NEW MD & S ) * 00163**9 C * (II) TEST ON MEAN DIFFERENCE * 00164**9 C * (III) CHI-SQUARE GOODNESS OF FIT TEST (PRECEDED BY HISTO- * 00165**9 C * GRAM) * 00166**9 C * (IV) SKEWNESS AND KURTOSIS * 00167**9 C * (V) TEST FOR OPTIMISM AND PESSIMISM * 00168**9 C * * 00169**9 C * 4. PLOT OF DATA TO SEE THEIR DISTRIBUTION (IF REQUIRED) * 00170**9 C * * 00171**9 C * 5. INTERMEDIATE RESULTS FROM SUBROUTINES 'APPROX' AND 'GRAV' * 00172**9 C * (IF REQUIRED). * 00173**9 C * * 00174**9 C * * 00175**9 C * A COMMON BLOCK CONTAINING THE LOCAL COORDINATES X,Y FOR * 00176**9 C * ALL DATA POINTS, AS WELL AS THE ORDER OF THE POLYNOMIAL, MUST * 00177**9 C * BE DECLARED BETWEEN THE MAIN PROGRAM AND THE FUNCTION- * 00178**9 C * SUBPROGRAM 'A'. * 00179**9 C * * 00180**9 C * ANOTHER COMMON BLOCK CONTAINING THE PARAMETERS OF THE * 00181**9 C * SUBPROGRAM 'COVAX' IS REQUIRED BETWEEN THE MAIN PROGRAM AND * 00182**9 C * THE SUBPROGRAM 'GRAV' IF 'COVAX' IS TO BE USED. * 00183**9 C * * 00184**9 C ******************************************************************* 00185**9 C 00186 IMPLICIT REAL*8(A-H,O-Z),LOGICAL (L) 00187**5 C 00188 C 00189**2 C DECLARE EXPLICITLY THE LABELLED COMMON BLOCK BETWEEN THIS MAIN-PROGRAM00190**2 C AND THE FUNCTION-SUBPROGRAM 'A' (WHICH WILL EVALUATE THE INDIVIDUAL 00191**2 C ELEMENTS OF THE DESIGN MATRIX-A OF THE BASE-FUNCTIONS OF THE POLYNOMIA00192**2 C 00193**2 COMMON /XYMAIN/ XFY,YAL,POLDEG 00194**2 COMMON /CMCOV/CI(12),CR(51),SIGMA0(300),SIGMA(300), 00195**5 *KI(25),N1,LOCAL 00196**5 C 00197**2 C SPECIFY THE DIMENSIONS OF THE INVOLVED ARRAYS 00198 C 00199 C XFY : THE X-PLANE-COORDINATE (+VE NORTH), IN KILOMETERS. 00200 C YAL : THE Y-PLANE-COORDINATE (+VE EAST) , IN KILOMETERS. 00201 C DGFA : THE FREE-AIR ANOMALY, IN MGALS. 00202 C WFA : THE WEIGHT FUNCTION OF THE FREE-AIR ANOMALY 00203 C HT : THE ELEVATION OF THE POINT ABOVE M.S.L., IN METERS. 00204 C WHT : THE WEIGHT FUNCTION OF THE ELEVATIONS (HEIGHTS) 00205 C CIJ : THE COEFFICIENTS OF THE APPROXIMATING POLYNOMIAL FOR THE 00206 C FREE-AIR ANOMALY SURFACE. 00207 C SIGCIJ: THE VARIANCE-COVARIANCE MATRIX OF 'CIJ' 00208 C TIJ : THE COEFFICIENTS OF THE APPROXIMATING POLYNOMIAL FOR THE 00209 C TOPOGRAPHY (HEIGHT) SURFACE. 00210 C SIGTIJ: THE VARIANCE-COVARIANCE MATRIX OF 'TIJ' . 00211 C IQUADR: THE QUADRANT INDICATOR. 00212 C U : THE ABSOLUTE TERMS OF THE NORMAL EQUATIONS. 00213 C 00214 DIMENSION RES(500),VRES(500),CSX(50),CXX(1500),CX(50,50), 00215**6 *CX2(50,50),CSXA(1,50),CSXB(50,1),CX1(50,50),CSX1(50,1) 00216**6 DIMENSION XFY(500),YAL(500),XI(9),YI(9),WTT(500),MINE(25) 00217**6 DIMENSION DGFA(500),DGBG(500),WFA(500),WFB(500) 00218**6 DIMENSION ALATT(500),ALOGG(500),DGA(500),HTA(500),HAC(500), 00219**6 * DIS(500),XLT(500),XLG(500),DG(500),DGC(500),TDS(500), 00220**6 *ALTT(500),ALGG(500),DGG(500),DGGC(500),TDSS(500),ALT(60), 00221**6 *ALG(60),DG3(200),DIT(100),HT(500),TX(5),TY(5),PWT(50),GD(50) 00222**8 DIMENSION CIJ(25),SIGCIJ(25,25),U(25),CVR(12),SIGCVR(12), 00223**6 *AVGPSI(12),CVRSUM(12),VCVRSM(12),DG4(200),SEG(10),SUPDIS(50), 00224**8 *SNORM(200),Z(200) 00225**8 DATA GM,RE/3.98D14,6371.0D3/, 00226**5 *D0,D1,D2/0.0D0,1.0D0,2.0D0/,PI/3.1415926535D0/ 00227**5 DOUBLE PRECISION DSIN,DCOS,DSQRT,DABS,DSIGN 00228 REAL*4 AA(8) 00229**5 INTEGER HVEC(20),IQUADR(4),POLDEG,NB(11),ICVR(12),ICHI(10),JNC(50)00230**7 C 00231 C 00232 C SET-UP VALUES FOR THE DEGREE OF THE POLYNOMIAL 'POLDEG' (FOR THE FREE-00233 C ANOMALY AND THE HEIGHT), AND THE MAXIMUM DIMENSIONAL SPACE ALLOCATED 00234 C FOR EACH ARRAY 'MAXDIM', AND THE DESIRED INTEGRATION DISTANCE 'DIST' 00235 C (IN KILOMETERS) ALONG THE LEVELLING LINE, AND THE REASONABLE SCALE-FAC00236 C 'SCALXY' TO BE MULTIPLIED BY THE LINEAR UNITS; SPECIFICALLY: THE X & Y00237 C PLANE-COORDINATES, THE INTEGRATION LENGTH 'DIST' OF THE LEVELLING LINE00238 C AND THE MEAN RADIUS OF CURVATURE 'R0M' OF THE ELLIPSOID- IN ORDER TO G00239 C A WELL-CONDITIONED MATRIX-Q OF THE NORMAL EQUATIONS FOR THE INVERSION 00240 C PROCESS. 00241 C 00242 READ(5,37) MINE 00243**4 37 FORMAT(25A1) 00244**4 WRITE(6,38) (MINE(I),I=1,25) 00245**4 38 FORMAT(/,50X,25A1//) 00246**4 READ(5,30) POLDEG,EXTENT,EXTEND,SCALXY,MAXDIM,N11 00247**5 30 FORMAT(I1,3D10.2,2I3) 00248**3 WRITE(6,34) POLDEG,EXTENT,EXTEND,MAXDIM,SCALXY 00249**3 34 FORMAT(25X,'DEGREE OF POLYNOMIAL=',I5,/,25X,'INITIAL EXTENT=', 00250**3 *F10.2/,25X,'INCREMENT=',F10.2/,25X,'MAX. # OF POINTS=',I5,//, 00251**4 *25X,'SCALE FACTOR=',F10.2//) 00252**3 READ(5,35) IPRNT,IWM,IPW,ISCRN,IREG 00253**5 35 FORMAT(5I2) 00254**5 READ(5,41) DLIMIT,ICOMB,JCOMB 00255**7 41 FORMAT(D10.2,2I5) 00256**7 WRITE(6,36) IPRNT,IWM,IPW,ISCRN,DLIMIT,IREG,ICOMB,JCOMB 00257**7 36 FORMAT(25X,'PRINT INTERMEDIATE RESULTS=',I4/,25X,'WEIGHTED MEAN PR00258**4 *EDICTION=',I4/,25X,'POWER FOR WTD MEAN PREDICTION=',I4/,25X, 00259**5 *'SCREEN DATA=',I4,10X,'LIMIT FOR DATA SCREENING=',F5.1//, 00260**5 *25X,'USE HEIGHT INDEPENDENT DATA(REGRESSION)=',I4/,30X,18('-')//, 00261**5 *23X,'SUPPLEMENTARY PREDICTION WITH COLLOCATION=',I4//,23X,'SUPPLEM00262**7 *ENTARY PREDICTION WITH WTD. MEANS=',I4//) 00263**7 C 00264**7 C... READ REGRESSION PARAMETERS: 00265**5 C 00266**5 READ(5,42) B0,AH0 00267**5 42 FORMAT(2D14.7) 00268**5 WRITE(6,43) B0,AH0 00269**5 43 FORMAT(9X,'REGRESSION PARAMETERS:'/10X,'B0,AH0=',2D14.7//) 00270**5 C... 00271**5 EXTEND=EXTEND*110.4656 00272**4 MAXD=MAXDIM 00273**4 EXTT=EXTENT 00274**4 C SET-UP THE NECESSARY VALUE OF THE SCALE-FACTOR 'SCALE', WHICH WILL BE 00275 C MULTIPLIED (IF REQUIRED) BY THE NORMAL-EQUATIONS MATRIX-Q, TO OVERCOME00276 C THE OVERFLOW (OR UNDERFLOW) PROBLEM DURING THE INVERSION PROCESS. 00277 C 00278 C 00279 C COMPUTE THE TOTAL NUMBER OF UNKNOWN COEFFICIENTS (WHICH EQUALS THE 00280 C NUMBER OF BASE-FUNCTIONS) FOR THE APPROXIMATING POLYNOMIAL . 00281 C 00282 MPOLY=(POLDEG+1)*(POLDEG+1) 00283 RLIMIT=60.D0 00284**5 MHALF=MPOLY/2 00285**5 MPY=MHALF*2 00286**5 IF(MPY.EQ.MPOLY) MHALF=MHALF-1 00287**5 C... 00288**5 C READ THE PARAMETERS FOR SUBROUTINES COVV; GRAV: 00289**5 C 00290**5 103 READ(5,46) S,ACV,KI(3),KI(4),LOCAL,N,KT,LSUM,N2,HMAX,LAST1 00291**5 46 FORMAT(2D14.7,2I5,L2,2I5,L2,I5,D14.7,L2) 00292**5 READ(5,33) ICOVAR,IFILE,IUNIT2 00293**6 33 FORMAT(2I2,I4) 00294**6 C... READ VARIABLES FOR CHI-SQUARE TEST;TEST ON VARIANCE: 00295**8 READ(5,1000) NKI,CHI,CHI1,CHI2 00296**8 1000 FORMAT(I2,3D12.5) 00297**8 C... READ FACTOR FOR REJECTION OF OUTLIERS: 00298**8 READ(5,8) REJ 00299**8 8 FORMAT(D12.4) 00300**8 READ(5,1500) FACTOR 00301**7 1500 FORMAT(D12.5) 00302**7 IF(NKI.EQ.2) GO TO 1001 00303**7 SEG(1)=0.431D0 00304**7 SEG(2)=-0.431D0 00305**7 GO TO 1002 00306**7 1001 SEG(1)=0.0D0 00307**7 1002 CONTINUE 00308**7 IF(ICOVAR.EQ.1) GO TO 62 00309**6 IF(N.LT.2) N2=2 00310**5 IF(N2.GT.2001) N2=2000 00311**5 N2=N2+1 00312**5 KI(5)=KT 00313**5 WRITE(6,47) S,ACV,KI(3),KI(4),N,KT 00314**5 47 FORMAT('PARAMETERS SPECIFYING THE MODEL DEGREE VARIANCES:',/, 00315**5 *' S,A= ',2D14.7,/,' K0-K3,N,KT= -2 -1 ',4I5/) 00316**5 DO 48 I=1,300 00317**5 48 SIGMA0(I)=0.0D0 00318**5 RB2=RE*RE*S 00319**5 LMODEL=N.EQ.0 00320**5 CI(8)=ACV*RB2*1.0D-10 00321**5 CI(10)=S 00322**5 IF(N.NE.0) GO TO 49 00323**5 LOCAL=.TRUE. 00324**5 N=2 00325**5 49 N1=N+1 00326**5 IF(.NOT.LOCAL) READ(5,65) (SIGMA0(I),I=1,N1) 00327**5 65 FORMAT(9F8.5) 00328**5 COKO=1.0D0 00329**5 DO 66 I=1,N1 00330**5 66 SIGMA0(I)=SIGMA0(I)*COKO 00331**5 IF(.NOT.LOCAL) WRITE(6,67) (SIGMA0(I),I=1,N1) 00332**5 67 FORMAT(' EMPIRICAL ANOMALY DEGREE VARIANCES IN UNITS OF MGAL**2:',00333**5 */,25(12F6.2/)) 00334**5 N2=2001 00335**5 KI(6)=3 00336**5 KI(7)=3 00337**5 CR(2)=0.0D0 00338**5 CR(3)=0.0D0 00339**5 CR(10)=GM/(RE+CR(2))**2 00340**5 CR(11)=GM/(RE+CR(3))**2 00341**5 62 CONTINUE 00342**6 C 00343 C SET-UP VALUES FOR SOME CONVERSION CONSTANTS. 00344 C 00345 RHODEG=180.D 0/PI 00346 RHOSEC=RHODEG*3600.D 0 00347 FEETMS=0.3048D 0 00348 CONVRT=110.4656D 00 00349**4 IV=0 00350**3 IJ=0 00351**3 SUMM=0.0D 00 00352**4 SUM=0.0D 00 00353**4 C 00354**4 C... DEFINE THE PARAMETERS OF THE REFERENCE ELLIPSOID IN CANADA(SEMI- 00355**4 C MAJOR&MINOR AXES)-THE 1927-NAD OF CLARK 1866 ELLIPSOID IN KM 00356**4 C 00357**4 ANAD=6378.2064D 00 00358**4 BNAD=6356.5838D 00 00359**4 AN=ANAD+BNAD 00360**4 BN=ANAD-BNAD 00361**4 ANA=ANAD*ANAD 00362**4 ECC1SQ=(AN*BN)/ANA 00363**4 TESTV=1.0D-16 00364**4 IF(ICOVAR.EQ.1) GO TO 997 00365**6 C 00366**4 C... PRINT SOME HEADINGS FOR THE JOB: 00367**4 C 00368**4 IF((ICOMB.EQ.1).OR.(JCOMB.EQ.1)) GO TO 69 00369**7 WRITE(6,32) 00370**4 32 FORMAT(//,35X,'GRAVITY ANOMALY PREDICTION AT BENCH MARKS'/,35X, 00371**4 *41('=')//,10X,'(1)',5X,'(2)',7X,'(3)',8X,'(4)',7X,'(5)',8X,'(6)', 00372**4 *7X,'(7)',8X,'(8)',7X,'(9)',7X,'(10)'/,10X,'BM-#',2X,'LATITUDE',1X,00373**4 *'LONGITUDE',2X,'OBSERVED',2X,'PREDICTED',2X,'SIG(OBS)',1X,'SIG(PRO00374**4 *P)',2X,'SIG(DIFF)',2X,'(4)-(5)',3X,'(9)NORMD') 00375**4 GO TO 73 00376**6 69 CONTINUE 00377**6 WRITE(6,70) 00378**6 70 FORMAT(//,35X,'GRAVITY ANOMALY PREDICTION AT BENCH MARKS'/,35X, 00379**6 *41('=')//,1X,'(1)',4X,'(2)',7X,'(3)',8X,'(4)',7X,'(5)',6X,'(6)', 00380**6 *7X,'(7)',7X,'(8)',7X,'(9)',6X,'(10)',5X,'(11)',7X,'(12)',/,'BM-#',00381**6 *2X,'LATITUDE',1X,'LONGITUDE',2X,'OBSERVED',2X,'SIG(OBS)',1X, 00382**6 *'PREDCTD-1',2X,'SIGMA-1',2X,'PREDCTD-2',2X,'SIGMA-2',2X,'(4)-(6)',00383**6 *3X,'(4)-(8)',2X,'(11)NORMD') 00384**6 73 CONTINUE 00385**6 C 00386 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*00387 C 00388 997 CONTINUE 00389**6 IF(ICOVAR.EQ.0) GO TO 999 00390**6 IF(IFILE.EQ.1) REWIND IUNIT2 00391**6 IF(IFILE.EQ.1) IFL=0 00392**6 DO 998 I=1,12 00393**6 CVRSUM(I)=0.D0 00394**6 VCVRSM(I)=0.D0 00395**6 998 CONTINUE 00396**6 MCV=0 00397**6 ITWM=0 00398**6 ITSUM=0 00399**6 999 CONTINUE 00400 SUMWT=0.0D 00 00401**4 WTSUM=0.0D 00 00402**4 SUMSIG=0.0D 00 00403**4 WTSQSM=0.0D 00 00404**4 ITSUM=ITSUM+ITWM 00405**6 C 00406 C READ-IN THE LATITUDE (+VE N) AND THE LONGITUDE (+VE W ) IN DECIMAL DEG00407 C OF THE CENTRE OF THE DESIRED DEGREE-SQUARE -UNDER INVESTIGATION 00408 C NOTE THAT THIS CENTRE WILL BE REFERRED TO AS THE POINT OF INTEREST 00409 C 00410 READ(5,31,END=6000) ID,(AA(I),I=1,8) 00411**5 IF(ICOVAR.EQ.1) WRITE(6,68) ID,(AA(I),I=1,2) 00412**6 68 FORMAT(//,5X,I5,2F10.2/,5X,25('*')//) 00413**6 31 FORMAT(I5,5F10.5,2F6.2,F8.2) 00414**4 PHI=AA(1) 00415**5 ALONG=AA(2) 00416**5 DGBK=AA(4) 00417**5 HT1=AA(5) 00418**5 WFBINV=0.05D 00**2+0.1967D 00**2*AA(8) 00419**5 C 00420**5 IF(IREG.EQ.1) GO TO 44 00421**5 GO TO 45 00422**5 44 CONTINUE 00423**5 DGBK=AA(3) 00424**5 WFBINV=AA(6) 00425**5 HT2=AA(8) 00426**5 45 CONTINUE 00427**5 C CONVERT THE LATITUDE OF THE GIVEN POINT OF INTEREST TO RADIANS 00428 C 00429 PHI0=PHI/RHODEG 00430 C 00431 C COMPUTATION OF THE MEAN RADIUS OF CURVATURE (IN KILOMETERS) OF THE NAD00432 C ELLIPSOID, AT THE LATITUDE (PHI-DEGREES) OF THE GIVEN POINT OF INTERES00433 C 00434 DENO=1.0D 0-ECC1SQ*DSIN(PHI0)**2 00435 RCN=ANAD/DSQRT(DENO) 00436 RCM=RCN*(1.D 0-ECC1SQ)/DENO 00437 R0M=DSQRT(RCN*RCM) 00438 C 00439 C COMPUTE THE CONVERSION FACTORS FROM CURVILINEAR DISTANCES ON THE ELLIP00440 C IN THE LATITUDE AND LONGITUDE DIRECTIONS, TO THE CORRESPONDING LINEAR 00441 C DISTANCES ON THE PLANE IN X AND Y DIRECTIONS, RESPECTIVELY. 00442 C 00443 FYFACT=R0M/RHODEG 00444 ALFACT=DCOS(PHI0)*FYFACT 00445 C 00446 C SET-UP THE ALLOWABLE INTERPOLATION BOUNDARIES (IN BOTH LATITUDE AND 00447 C LONGITUDE) WITHIN EACH OF THE SEVEN PARTITIONED DATA SETS. 00448 C 00449 PHIL21=41.5D 0 00450 PHIU21=63.5D 0 00451 ALGL21=52.5D 0 00452 ALGU21=73.5D 0 00453 PHIL22=43.5D 0 00454 PHIU22=63.5D 0 00455 ALGL22=73.5D 0 00456 ALGU22=91.5D 0 00457 PHIL23=48.5D 0 00458 PHIU23=63.5D 0 00459 ALGL23=91.5D 0 00460 ALGU23=118.5D 0 00461 PHIL24=48.5D 0 00462 PHIU24=63.5D 0 00463 ALGL24=118.5D 0 00464 ALGU24=140.5D 0 00465 PHIL25=63.5D 0 00466 PHIU25=70.5D 0 00467 ALGL25=61.5D 0 00468 ALGU25=142.5D 0 00469 PHIL16=70.5D 0 00470 PHIU16=80.5D 0 00471 ALGL16=66.5D 0 00472 ALGU16=142.5D 0 00473 PHIL17=80.5D 0 00474 PHIU17=82.5D 0 00475 ALGL17=51.5D 0 00476 ALGU17=78.5D 0 00477 C 00478 C TESTING THE GIVEN LOCATION (LONGITUDE & LATITUDE) OF THE BENCH MARK 00479 C AGAINST THE ABOVE SPECIFIED INTERPOLATION BOUNDARIES , TO DETERMINE TH00480 C APPROPRIATE FILE (WHOSE DATA REFERENCE NUMBER IS DENOTED BY 'IUNIT'), 00481 C THAT CONTAINS THE NECESSARY DATA REQUIRED FOR INTERPOLATION. 00482 C 00483 IF(ALONG.GE.ALGL21.AND.ALONG.LE.ALGU21) GO TO 1 00484 IF(ALONG.GE.ALGL22.AND.ALONG.LE.ALGU22) GO TO 2 00485 IF(ALONG.GE.ALGL23.AND.ALONG.LE.ALGU23) GO TO 3 00486 IF(ALONG.GE.ALGL24.AND.ALONG.LE.ALGU24) GO TO 4 00487 C 00488 C PRINT-OUT A MESSAGE TO THE USER IF THE INPUTTED BENCH MARK IS LOCATED00489 C OUTSIDE THE RANGE OF OBSERVED GRAVITY DATA EXISTING ON THE AVAILABLE 00490 C DATA FILES. 00491 C 00492 WRITE(6,10) 00493 10 FORMAT('1'///,5X,'SORRY: THE GIVEN POINT IS LOCATED IN AN AREA WHE00494 1RE NO OBSERVED GRAVITY DATA EXIST ON THE AVAILABLE FILES'//) 00495 GO TO 9999 00496 1 IF(PHI.GE.PHIL21.AND.PHI.LE.PHIU21) GO TO 11 00497 IF(PHI.GE.PHIL25.AND.PHI.LE.PHIU25) GO TO 15 00498 IF(PHI.GE.PHIL16.AND.PHI.LE.PHIU16) GO TO 6 00499 IF(PHI.GE.PHIL17.AND.PHI.LE.PHIU17) GO TO 7 00500 WRITE(6,10) 00501 GO TO 9999 00502 2 IF(PHI.GE.PHIL22.AND.PHI.LE.PHIU22) GO TO 12 00503 IF(PHI.GE.PHIL25.AND.PHI.LE.PHIU25) GO TO 15 00504 IF(PHI.GE.PHIL16.AND.PHI.LE.PHIU16) GO TO 6 00505 IF(PHI.GE.PHIL17.AND.PHI.LE.PHIU17) GO TO 7 00506 WRITE(6,10) 00507 GO TO 9999 00508 3 IF(PHI.GE.PHIL23.AND.PHI.LE.PHIU23) GO TO 13 00509 IF(PHI.GE.PHIL25.AND.PHI.LE.PHIU25) GO TO 15 00510 IF(PHI.GE.PHIL16.AND.PHI.LE.PHIU16) GO TO 6 00511 WRITE(6,10) 00512 GO TO 9999 00513 4 IF(PHI.GE.PHIL24.AND.PHI.LE.PHIU24) GO TO 14 00514 IF(PHI.GE.PHIL25.AND.PHI.LE.PHIU25) GO TO 15 00515 IF(PHI.GE.PHIL16.AND.PHI.LE.PHIU16) GO TO 6 00516 WRITE(6,10) 00517 GO TO 9999 00518 C 00519 C SPECIFY THE APPROPRIATE DATA-SET REFERENCE-NUMBER ,WHICH CONTAINS THE 00520 C NECESSARY GRAVITY DATA AROUND THE POINT OF INTEREST. 00521 C 00522 11 IUNIT=21 00523 GO TO 5 00524 12 IUNIT=22 00525 GO TO 5 00526 13 IUNIT=23 00527 GO TO 5 00528 14 IUNIT=24 00529 GO TO 5 00530 15 IUNIT=25 00531 GO TO 5 00532 6 IUNIT=16 00533 GO TO 5 00534 7 IUNIT=17 00535 5 CONTINUE 00536 C 00537 39 CONTINUE 00538**3 RADIUS=EXTENT/2.D0 00539**3 RR=RADIUS/DSQRT(2.D0) 00540**3 XI(1)=PHI 00541**3 XI(2)=PHI+RR 00542**3 XI(3)=PHI+RADIUS 00543**3 XI(4)=XI(2) 00544**3 XI(5)=XI(1) 00545**3 XI(6)=PHI-RR 00546**3 XI(7)=PHI-RADIUS 00547**3 XI(8)=XI(6) 00548**3 XI(9)=XI(1) 00549**3 C 00550**3 YI(1)=ALONG-RADIUS 00551**3 YI(2)=ALONG-RR 00552**3 YI(3)=ALONG 00553**3 YI(4)=ALONG+RR 00554**3 YI(5)=ALONG+RADIUS 00555**3 YI(6)=YI(4) 00556**3 YI(7)=YI(3) 00557**3 YI(8)=YI(2) 00558**3 YI(9)=YI(1) 00559**3 C 00560**3 C 00561 C INITIALIZE ZERO VALUES FOR THE 4-QUADRANTS INDICATOR 'IQUADR' 00562 C 00563 DO 40 I=1,4 00564 40 IQUADR(I)=0 00565 C 00566 C READ-IN THE DATA RESIDING ON DISC AS DATA SET NO. 'IUNIT' CONTAINING 00567 C FOR EACH POINT:- LAT.+VE N AND LONG.+VE W ,BOTH IN DECIMAL DEGREES ; 00568 C OBSERVED GRAVITY IN GALS ; FREE-AIR GRAVITY ANOMAL IN MGALS ; ELEVATIO00569 C AND ITS ACCURACY IN FEET. 00570 C 00571 150 FORMAT(9X,I5,8F10.5) 00572**5 152 FORMAT(106X,2I10) 00573**5 153 FORMAT(106X,F10.2,I10) 00574**5 154 FORMAT(106X,I10) 00575**5 367 FORMAT(//) 00576**5 REWIND IUNIT 00577 J=0 00578**3 CALL CPUTIM(ITI,JJR) 00579**3 50 READ(IUNIT,52,END=60) (AA(I),I=1,8) 00580**5 IF((AA(6).EQ.0.).OR.(AA(7).EQ.0.).OR.(AA(8).EQ.0.)) GO TO 50 00581**5 52 FORMAT(5F12.5,2F6.2,F8.2) 00582**4 ALAT=AA(1) 00583**5 ALOG=AA(2) 00584**5 FA=AA(3) 00585**5 BA=AA(4) 00586**5 AHT=AA(5) 00587**5 SFA=AA(6) 00588**5 SHTA=AA(8) 00589**5 SBA=0.05D 00**2+0.1967D 00**2*SHTA 00590**4 C 00591 C CHECK, SELECT, AND STORE THE INFORMATION FOR ONLY THOSE POINTS LOCATED00592 C WITHIN THE SPECIFIED FOUR DEGREE-SQUARES AREA AROUND THE POINT OF INTE00593 C 00594 IF((ALAT.EQ.PHI).AND.(ALOG.EQ.ALONG)) GO TO 50 00595**4 CALL INPOLY(XI,YI,N11,ALAT,ALOG,IN) 00596**5 IF(IN.EQ.1) GO TO 51 00597**3 GO TO 50 00598**3 51 CONTINUE 00599**3 J=J+1 00600**3 ALATT(J)=ALAT 00601**3 ALOGG(J)=ALOG 00602**3 DGA(J)=BA 00603**3 HTA(J)=SBA 00604**3 IF(IREG.EQ.1) GO TO 155 00605**5 GO TO 156 00606**5 155 CONTINUE 00607**5 DGA(J)=FA-B0*(AHT-AH0) 00608**5 HTA(J)=SFA-SHTA*B0**2 00609**5 156 CONTINUE 00610**5 DIFF1=ALAT-PHI 00611**4 DIFF2=ALONG-ALOG 00612**4 DS1=DIFF1*CONVRT 00613**4 DS2=DIFF2*CONVRT*DCOS(PHI0) 00614**4 DS1SQ=DS1*DS1 00615**4 DS2SQ=DS2*DS2 00616**4 DIS(J)=DSQRT(DS1SQ+DS2SQ) 00617**4 IF(ICOVAR.EQ.1) GO TO 63 00618**6 GO TO 64 00619**6 63 CONTINUE 00620**6 XFY(J)=DS1 00621**6 YAL(J)=DS2 00622**6 IF(J.EQ.500) GO TO 60 00623**6 64 CONTINUE 00624**6 GO TO 50 00625**3 60 CONTINUE 00626**3 CALL CPUTIM(ITI,JJR) 00627**3 WRITE(6,152) ITI,JJR 00628**3 WRITE(6,154) J 00629**5 IF(ICOVAR.EQ.1) GO TO 5200 00630**6 IF(ISCRN.EQ.0) GO TO 61 00631**5 C 00632**5 C... SCREEN DATA COLLECTED; WHERE CLUSTER OF POINTS EXISTS, REPLACE 00633**5 C WITH THE STRAIGHT MEAN: 00634**5 C 00635**5 IQRR=1 00636**5 KL=0 00637**5 NII=5 00638**5 SMALL=0.0025D 00 00639**5 221 CONTINUE 00640**5 C... 00641**5 TX(1)=PHI 00642**5 TX(2)=TX(1) 00643**5 TX(3)=TX(1)+RADIUS 00644**5 IF((IQRR.EQ.3).OR.(IQRR.EQ.4)) TX(3)=TX(1)-RADIUS 00645**5 TX(4)=TX(3) 00646**5 TX(5)=TX(1) 00647**5 C... 00648**5 TY(1)=ALONG 00649**5 TY(2)=TY(1)+RADIUS 00650**5 IF((IQRR.EQ.2).OR.(IQRR.EQ.3)) TY(2)=TY(1)-RADIUS 00651**5 TY(3)=TY(2) 00652**5 TY(4)=TY(1) 00653**5 TY(5)=TY(1) 00654**5 C... 00655**5 IT=0 00656**5 DO 222 I=1,J 00657**5 IF(ALATT(I).EQ.0.) GO TO 222 00658**5 AAL=ALATT(I) 00659**5 AAG=ALOGG(I) 00660**5 CALL INPOLY(TX,TY,NII,AAL,AAG,IN) 00661**5 IF(IN.EQ.1) GO TO 223 00662**5 GO TO 222 00663**5 223 CONTINUE 00664**5 IF(HTA(I).LE.SMALL) GO TO 222 00665**5 IT=IT+1 00666**5 XLT(IT)=ALATT(I) 00667**5 XLG(IT)=ALOGG(I) 00668**5 DG(IT)=DGA(I) 00669**5 DGC(IT)=HTA(I) 00670**5 TDS(IT)=DIS(I) 00671**5 ALATT(I)=0. 00672**5 222 CONTINUE 00673**5 IF(IT.EQ.0) GO TO 229 00674**5 ITT=1 00675**5 ICT=1 00676**5 K=1 00677**5 IF(IT.EQ.1) GO TO 224 00678**5 KP1=K+1 00679**5 GD(ICT)=DG(K) 00680**5 IF(DGC(K).LE.TESTV) GO TO 231 00681**5 GO TO 232 00682**5 231 PWT(ICT)=1.D 16 00683**5 232 PWT(ICT)=1.D 0/DGC(K) 00684**5 DSM1=XLT(K)*PWT(ICT) 00685**5 DSM2=XLG(K)*PWT(ICT) 00686**5 DSM3=DG(K)*PWT(ICT) 00687**5 DSM4=PWT(ICT) 00688**5 DSM5=DGC(K) 00689**5 228 CONTINUE 00690**5 DO 224 I=KP1,IT 00691**5 IF(XLT(K).EQ.0.) GO TO 227 00692**5 IF(XLT(I).EQ.0.) GO TO 224 00693**5 DT1=(XLT(K)-XLT(I))*CONVRT 00694**5 DT2=(XLG(I)-XLG(K))*CONVRT*DCOS(PHI0) 00695**5 DTEST=DSQRT(DT1**2+DT2**2) 00696**5 IF(DTEST.LE.DLIMIT) GO TO 225 00697**5 IF((ITT.GT.1).AND.(I.EQ.IT)) GO TO 226 00698**5 GO TO 224 00699**5 225 CONTINUE 00700**5 ICT=ICT+1 00701**5 IF(DGC(I).LE.TESTV) GO TO 233 00702**5 GO TO 234 00703**5 233 PWT(ICT)=1.D 16 00704**5 234 PWT(ICT)=1.D 0/DGC(I) 00705**5 GD(ICT)=DG(I) 00706**5 C 00707**5 DSM1=DSM1+XLT(I)*PWT(ICT) 00708**5 DSM2=DSM2+XLG(I)*PWT(ICT) 00709**5 DSM3=DSM3+DG(I)*PWT(ICT) 00710**5 DSM4=DSM4+PWT(ICT) 00711**5 DSM5=DSM5+DGC(I) 00712**5 C 00713**5 XLT(I)=0. 00714**5 ITT=ITT+1 00715**5 IF(I.EQ.IT) GO TO 226 00716**5 GO TO 224 00717**5 C... WEIGHTED MEAN: 00718**5 226 CONTINUE 00719**5 KL=KL+1 00720**5 ALTT(KL)=DSM1/DSM4 00721**5 ALGG(KL)=DSM2/DSM4 00722**5 DGG(KL)=DSM3/DSM4 00723**5 GDSM=0. 00724**5 DO 230 IM=1,ICT 00725**5 230 GDSM=GDSM+PWT(IM)*(GD(IM)-DGG(KL))**2 00726**5 COUNT=DFLOAT(ICT-1) 00727**5 DGGC(KL)=GDSM/(DSM4*COUNT) 00728**5 DT3=(ALTT(KL)-PHI)*CONVRT 00729**5 DT4=(ALONG-ALGG(KL))*CONVRT*DCOS(PHI0) 00730**5 TDSS(KL)=DSQRT(DT3**2+DT4**2) 00731**5 GO TO 227 00732**5 224 CONTINUE 00733**5 KL=KL+1 00734**5 ALTT(KL)=XLT(K) 00735**5 ALGG(KL)=XLG(K) 00736**5 DGG(KL)=DG(K) 00737**5 DGGC(KL)=DGC(K) 00738**5 TDSS(KL)=TDS(K) 00739**5 IF(KP1.EQ.IT) GO TO 235 00740**5 227 CONTINUE 00741**5 K=K+1 00742**5 KP1=K+1 00743**5 ITT=1 00744**5 ICT=1 00745**5 GD(ICT)=DG(K) 00746**5 PWT(ICT)=1.D0/DGC(K) 00747**5 DSM1=XLT(K)*PWT(ICT) 00748**5 DSM2=XLG(K)*PWT(ICT) 00749**5 DSM3=DG(K)*PWT(ICT) 00750**5 DSM4=PWT(ICT) 00751**5 DSM5=DGC(K) 00752**5 IF(K.LT.IT) GO TO 228 00753**5 235 CONTINUE 00754**5 IF(XLT(IT).EQ.0.) GO TO 229 00755**5 KL=KL+1 00756**5 ALTT(KL)=XLT(IT) 00757**5 ALGG(KL)=XLG(IT) 00758**5 DGG(KL)=DG(IT) 00759**5 DGGC(KL)=DGC(IT) 00760**5 TDSS(KL)=TDS(IT) 00761**5 229 CONTINUE 00762**5 IQRR=IQRR+1 00763**5 IF(IQRR.LE.4) GO TO 221 00764**5 J=KL 00765**5 C 00766**5 C... NOW SELECT DATA FOR PREDICTION: 00767**5 C 00768**5 WRITE(6,154) J 00769**5 DO 58 I=1,J 00770**5 ALATT(I)=ALTT(I) 00771**5 ALOGG(I)=ALGG(I) 00772**5 DGA(I)=DGG(I) 00773**5 HTA(I)=DGGC(I) 00774**5 DIS(I)=TDSS(I) 00775**5 58 CONTINUE 00776**5 CALL CPUTIM(ITI,JJR) 00777**5 WRITE(6,152) ITI,JJR 00778**5 61 CONTINUE 00779**5 TESTER=(5.D 00/60.D 00)*110.4656 00780**4 NP=0 00781**3 ITWM=0 00782**6 53 CONTINUE 00783**3 DO 54 I=1,J 00784**3 IF(NP.GE.MAXDIM) GO TO 78 00785**5 IF(ALATT(I).EQ.0.0D0) GO TO 54 00786**3 IF(DIS(I).GT.TESTER) GO TO 54 00787**3 IF(HTA(I).LE.SMALL) GO TO 54 00788**5 NP=NP+1 00789**3 ALT(NP)=ALATT(I) 00790**3 ALG(NP)=ALOGG(I) 00791**3 DGBG(NP)=DGA(I) 00792**3 WFB2=HTA(I) 00793**4 IF(JCOMB.EQ.1) SUPDIS(NP)=DIS(I) 00794**7 IF(IWM.EQ.1) GO TO 56 00795**4 WFB(NP)=1.D0/WFB2 00796**3 55 FORMAT(I5,5F10.2) 00797**3 GO TO 57 00798**4 56 CONTINUE 00799**4 CALL CPUTIM(ITI,JJR) 00800**6 IF(DIS(I).EQ.0.D0) GO TO 54 00801**4 EDIST=DIS(I)**IPW 00802**4 IF(EDIST.LE.TESTV) GO TO 540 00803**4 WTT(NP)=1.0D 00/EDIST 00804**4 GO TO 541 00805**4 540 WTT(NP)=1.0D 16 00806**4 541 CONTINUE 00807**4 WT=WTT(NP) 00808**4 WTSQ=WT*WT 00809**4 SUMWT=SUMWT+WT 00810**4 WTSUM=WTSUM+DGBG(NP)*WT 00811**4 SUMSIG=SUMSIG+WTSQ*WFB2 00812**4 WTSQSM=WTSQSM+WTSQ 00813**4 CALL CPUTIM(ITI,JJR) 00814**6 ITWM=ITWM+ITI 00815**6 IF((IWM.EQ.1).AND.(IPRNT.EQ.0)) GO TO 542 00816**5 57 CONTINUE 00817**4 C 00818**4 C 00819**3 XFY(NP)=(ALT(NP)-PHI)*FYFACT*SCALXY 00820**3 YAL(NP)=(ALONG-ALG(NP))*ALFACT*SCALXY 00821**3 IF(IPRNT.EQ.1) WRITE(6,198) NP,ALT(NP),ALG(NP),DGBG(NP),WFB2, 00822**5 *WFB(NP),XFY(NP),YAL(NP),DIS(I) 00823**4 198 FORMAT(9X,I5,8F10.5) 00824**4 C 00825**3 C 00826 C STORE ANY INTEGER NUMBER, SAY 1, IN THE POSITION 'JACK' OF THE QUADRAN00827 C INDICATOR ARRAY 'IQUADR', TO HELP TESTING FOR THE EXISTANCE OF AT LEAS00828 C ONE DATA POINT IN EACH QUADRANT AROUND THE POINT OF INTEREST 00829 C 00830 C 'JACK' TAKES THE VALUES: 1,2,3 OR 4 ,ACCORDING TO THE QUADRANTS WITH00831 C X,Y=(-,-),(-,+),(+,-) AND (+,+) RESPECTIVELY. 00832 C 00833 542 CONTINUE 00834**5 IF((ALT(NP).GE.PHI).AND.(ALG(NP).GE.ALONG)) IQUADR(1)=IQUADR(1)+1 00835**4 IF((ALT(NP).GE.PHI).AND.(ALG(NP).LT.ALONG)) IQUADR(2)=IQUADR(2)+1 00836**4 IF((ALT(NP).LT.PHI).AND.(ALG(NP).LT.ALONG)) IQUADR(3)=IQUADR(3)+1 00837**4 IF((ALT(NP).LT.PHI).AND.(ALG(NP).GE.ALONG)) IQUADR(4)=IQUADR(4)+1 00838**4 C 00839 ALATT(I)=0.0D0 00840**3 IF(NP.GE.MAXDIM) GO TO 72 00841**3 54 CONTINUE 00842**3 78 CONTINUE 00843**5 NUMDF=NP-MPOLY 00844**3 IF(NUMDF.LE.MHALF) GO TO 74 00845**5 C 00846 C COMPUTE AND PRINT-OUT THE TOTAL NUMBER OF AVAILABLE DATA POINTS 00847 C 00848 IF(IPRNT.EQ.0) GO TO 72 00849**5 WRITE(6,71) NP 00850 71 FORMAT(10X,'THE TOTAL NUMBER OF AVAILABLE DATA POINTS IS =',I5,2X,00851 1'POINTS'/) 00852 GO TO 72 00853**3 74 CONTINUE 00854**3 IF(NP.EQ.J) GO TO 9999 00855**3 IF(TESTER.GE.RLIMIT) GO TO 9999 00856**5 TESTER=TESTER+EXTEND 00857**3 GO TO 53 00858**3 72 CONTINUE 00859**3 WRITE(6,153) TESTER,NP 00860**3 C 00861 C TEST FOR THE EXISTANCE OF AT LEAST ONE DATA POINT IN EACH OF THE FOUR 00862 C QUADRANTS AROUND THE POINT OF INTEREST, OTHERWISE PRINT-OUT A MESSAGE 00863 C THE USER, AND THEN STOP EXECUTION. 00864 C 00865 IQVAL=IQUADR(1)*IQUADR(2)*IQUADR(3)*IQUADR(4) 00866 IF(IQVAL.EQ.0) GO TO 75 00867**4 GO TO 77 00868**4 75 CONTINUE 00869**4 IQR=0 00870**4 DO 76 I=1,4 00871**4 IF(IQUADR(I).EQ.0) IQR=IQR+1 00872**4 76 CONTINUE 00873**4 IF(IQR.GT.1) GO TO 990 00874**4 77 CONTINUE 00875**4 IF(IPRNT.EQ.0) GO TO 5000 00876**5 C***********************************************************************00877**3 C 00878**3 C CALL THE SUBROUTINE 'GCPLOT' TO PLOT A SKETCH SHOWING THE APPROXIMATE 00879**3 C DISTRIBUTION OF THE USED DATA POINTS WITHIN THE ONE-AREA-SQUARED 00880**3 C SPECIFIED AROUND THE POINT OF INTEREST, FOR THE INTERPOLATION PURPOSES00881**3 C 00882**3 AREA=EXTENT/2.D0 00883**3 PLOTFY= AREA *FYFACT*SCALXY 00884**3 PLOTAL= AREA *ALFACT*SCALXY 00885**3 XFYMIN=-PLOTFY 00886**3 YALMIN=-PLOTAL 00887**3 XFYMAX=+PLOTFY 00888**3 YALMAX=+PLOTAL 00889**3 C 00890**3 CALL GCPLOT(YAL,XFY,NP,YALMIN,XFYMIN,YALMAX,XFYMAX) 00891**3 C 00892**3 CMFLSC=( EXTD*60.)/15. 00893**3 WRITE(6,554) CMFLSC 00894**3 554 FORMAT(5X,'SCALE OF THE ABOVE PLOT IS: 1.0 CENTIMETER REPRESENTS',00895**3 1F5.2,2X,'MINUTES (IN BOTH LAT.& LONG.)'/) 00896**3 C***********************************************************************00897 5000 CONTINUE 00898**4 C 00899 C CALL THE SUBROUTINE 'APPROX' TO COMPUTE AND PRINt-OUT THE ESTIMATED 00900 C COEFFICIENTS OF THE BEST-FITTING POLYNOMIAL, APPROXIMATING THE FREE-AI00901 C GRAVITY SURFACE, ALONG WITH THEIR ASSOCIATED VARIANCE-COVARIANCE MATRI00902 C 00903 CALL CPUTIM(ITI,JJR) 00904**4 IV=IV+1 00905**3 199 FORMAT(9X,I5,8F10.5,I5,3F10.5) 00906**3 IF(IWM.EQ.1) GO TO 5001 00907**4 C 00908**4 C... 00909**4 CALL APPROX(NP,MPOLY,DGBG,WFB,SCALXY,CIJ,SIGCIJ,EVFCIJ,U,DETQ, 00910**4 *IDEXP,IPRNT,RES,VRES,ICOVAR,ICOMB,INDEX) 00911**6 IF(INDEX.EQ.1) GO TO 9999 00912**6 IF(IREG.EQ.1) GO TO 200 00913**5 GO TO 201 00914**5 200 CONTINUE 00915**5 CIJ(1)=CIJ(1)+B0*(HT1-AH0) 00916**5 SIGCIJ(1,1)=SIGCIJ(1,1)+HT2*B0**2 00917**5 201 CONTINUE 00918**5 C 00919 C EXTRACT AND PRINT-OUT THE PREDICTED VALUE OF THE FREE-AIR ANOMALY, ALO00920 C WITH ITS ESTIMATED STANDARD DEVIATION, AT THE CENTRE OF THE ONE-DEGREE00921 C SQUARE UNDER INVESTIGATION. 00922 C 00923 GO TO 5002 00924**4 5001 CONTINUE 00925**6 CALL CPUTIM(ITI,JJR) 00926**6 CIJ(1)=WTSUM/SUMWT 00927**6 WTSQSM=SUMWT*SUMWT 00928**4 SUMV=0.D0 00929**4 DO 5003 I=1,NP 00930**4 VV=DGBG(I)-CIJ(1) 00931**4 VVSQP=VV*VV*WTT(I) 00932**4 SUMV=SUMV+VVSQP 00933**4 5003 CONTINUE 00934**4 SJ=SUMSIG/WTSQSM 00935**4 SUMVD=(NP-1)*SUMWT 00936**4 SJD=SUMV/SUMVD 00937**4 SIGCIJ(1,1)=SJ 00938**4 SJ=DSQRT(SJ) 00939**4 SJD=DSQRT(SJD) 00940**4 5002 CONTINUE 00941**4 DG3(IV)=DGBK-CIJ(1) 00942**3 IF(ICOMB.EQ.0) SUMM=SUMM+DG3(IV) 00943**5 IF(ICOMB.EQ.0) SUM=SUM+DG3(IV)**2 00944**5 VK=DSQRT(WFBINV) 00945**4 FACTSQ=FACTOR*FACTOR 00946**8 SS=DSQRT(FACTSQ*SIGCIJ(1,1)+WFBINV) 00947**8 SNORM(IV)=DG3(IV)/SS 00948**7 TEST1=DABS(SNORM(IV)) 00949**7 IF((TEST1.LT.1.).AND.(ICOMB.EQ.0)) IJ=IJ+1 00950**5 SJ=DSQRT(SIGCIJ(1,1)) 00951**4 CALL CPUTIM(ITI,JJR) 00952**4 DG5=DG3(IV) 00953**7 SJ1=SJ 00954**6 IF((ICOMB.EQ.1).OR.(JCOMB.EQ.1)) GO TO 16 00955**7 IF(IWM.EQ.1) GO TO 19 00956**6 ITWM=ITI 00957**6 GO TO 18 00958**6 19 ITWM=ITWM+ITI 00959**6 18 CONTINUE 00960**6 WRITE(6,152) ITWM,JJR 00961**6 WRITE(6,9) ID,PHI,ALONG,DGBK,CIJ(1),VK,SJ,SS,DG3(IV),TEST1,NP 00962**8 9 FORMAT(9X,I5,9F10.5,I5) 00963**4 GO TO 9999 00964**6 16 CONTINUE 00965**6 IF(IPRNT.EQ.1) WRITE(6,367) 00966**5 IF(ICOMB.EQ.1) GO TO 5110 00967**7 C 00968**7 C<<< USE WEIGHTED MEAN METHOD ON RESIDUALS FROM POLYNOMIAL APPROX: 00969**7 C 00970**7 CALL WTDAVG(NP,RES,VRES,SUPDIS,IPW,TESTV,DGP,VAR) 00971**7 GO TO 5120 00972**7 5110 CONTINUE 00973**7 C... 00974**5 C USE COLLOCATION ON RESIDUALS FROM POLYNOMIAL APPROXIMATION: 00975**5 C 00976**5 PHI=PHI/RHODEG 00977**5 ALONG=ALONG/RHODEG 00978**5 VRFACT=EVFCIJ*NUMDF/NP 00979**5 5105 FORMAT(9X,I5,4F10.5) 00980**5 DO 5100 IJK=1,NP 00981**5 IF(IPRNT.EQ.1) WRITE(6,5105) IJK,ALT(IJK),ALG(IJK),RES(IJK), 00982**5 *VRES(IJK) 00983**5 ALT(IJK)=ALT(IJK)/RHODEG 00984**5 ALG(IJK)=ALG(IJK)/RHODEG 00985**5 5100 CONTINUE 00986**5 IF(IPRNT.EQ.1) WRITE(6,367) 00987**5 NBRR=(NP*(NP+1))/2 00988**5 C... 00989**5 C CALL SUBROUTINE GRAV FOR SUPPLEMENTARY PREDICTION: 00990**5 C 00991**5 CALL GRAV(NP,VRES,ALT,ALG,RES,CXX,CSX,CSXA,CSXB,PHI,ALONG,DGP, 00992**5 *VAR,NBRR,CX,CX2,GG,CX1,CSX1,IPRNT,INDEX) 00993**6 IF(INDEX.EQ.1) GO TO 9999 00994**6 PHI=PHI*RHODEG 00995**6 ALONG=ALONG*RHODEG 00996**5 5120 CONTINUE 00997**7 DGP1=CIJ(1)+DGP 00998**5 SJSQ=VAR+SIGCIJ(1,1) 00999**5 SJ=DSQRT(SJSQ) 01000**5 DG3(IV)=DGBK-DGP1 01001**5 SUMM=SUMM+DG3(IV) 01002**5 SUM=SUM+DG3(IV)**2 01003**5 SS1=DSQRT(WFBINV+SJSQ) 01004**5 TEST1=DABS(DG3(IV)/SS1) 01005**5 IF(TEST1.LT.1.) IJ=IJ+1 01006**5 CALL CPUTIM(ITI,JJR) 01007**6 ITWM=ITI 01008**6 WRITE(6,152) ITWM,JJR 01009**6 WRITE(6,17) ID,PHI,ALONG,DGBK,VK,CIJ(1),SJ1,DGP1,SJ,DG5,DG3(IV), 01010**7 *TEST1 01011**6 17 FORMAT(I4,11F10.5) 01012**6 GO TO 9999 01013**3 5200 CONTINUE 01014**6 C 01015**9 C... OBTAIN RESIDUALS FROM SURFACE FIT TO DATA FOR COMPUTATION OF 01016**9 C COVARIANCE FUNCTION: 01017**9 C 01018**9 NUMDF=J-MPOLY 01019**6 IF(NUMDF.LE.MHALF) GO TO 9999 01020**6 ND=J 01021**6 CALL APPROX(ND,MPOLY,DGA,HTA,SCALXY,CIJ,SIGCIJ,EVFCIJ,U,DETQ,IDEXP01022**6 *,IPRNT,RES,VRES,ICOVAR,ICOMB,INDEX) 01023**6 IF(INDEX.EQ.1) GO TO 9999 01024**6 DO 5201 I=1,ND 01025**6 WRITE(6,5202) ALATT(I),ALOGG(I),RES(I),VRES(I),XFY(I),YAL(I) 01026**6 IF(IFILE.EQ.1) WRITE(IUNIT2,5206) ALATT(I),ALOGG(I),RES(I),VRES(I)01027**6 IF(IFILE.EQ.1) IFL=IFL+1 01028**6 5201 CONTINUE 01029**6 5202 FORMAT(6F12.5) 01030**6 WRITE(6,5205) EVFCIJ 01031**6 WRITE(6,5203) 01032**6 5203 FORMAT(//) 01033**6 5205 FORMAT(/,75X,F13.5) 01034**6 5206 FORMAT(4F20.5) 01035**6 GO TO 9999 01036**6 5204 CONTINUE 01037**6 CALL PRECOV(ND,ALATT,ALOGG,RES,VRES,CVR,SIGCVR,CLASS,AVGPSI,ICVR) 01038**6 C... 01039**6 DO 5210 I=1,12 01040**6 CVRSUM(I)=CVRSUM(I)+CVR(I) 01041**6 VCVRSM(I)=VCVRSM(I)+SIGCVR(I) 01042**6 5210 CONTINUE 01043**6 MCV=MCV+1 01044**6 GO TO 9999 01045**6 6000 CONTINUE 01046**4 IF(ICOVAR.EQ.1) GO TO 950 01047**6 TIME=ITSUM*1.D-4 01048**6 DIV=DFLOAT(IV) 01049**4 SUMM1=SUMM/DIV 01050**4 SUM1=(SUM*DIV-SUMM**2)/(DIV*(DIV-1)) 01051**4 SUM2=DSQRT(SUM1) 01052**4 IPCT=0 01053**7 VSM=SUM/DIV 01054**8 RMS2=DSQRT(VSM) 01055**7 DO 804 I=1,IV 01056**7 TEST4=DG3(I)/RMS2 01057**7 DG4(I)=TEST4 01058**7 IF(TEST4.LT.1.) IPCT=IPCT+1 01059**7 804 CONTINUE 01060**7 WRITE(6,800) SUMM1,SUM1,SUM2,RMS2 01061**7 800 FORMAT(//,50X,'MEAN DIFFERENCE =',F12.5//,50X,'MEAN SQ. ERROR 01062**7 * =',F12.5/,50X,'ROOT MEAN SQ. ERROR=',F12.5//,50X,'RMS WRPT POP-01063**7 *MEAN =',F12.5) 01064**7 TEST2=(DFLOAT(IJ)*100.)/DIV 01065**4 TEST5=(IPCT*100.)/DIV 01066**7 WRITE(6,801) IV,IJ,TEST2,TEST5 01067**7 801 FORMAT(//,50X,'TOT. # OF PREDICTIONS=',I10/,50X,'TOT. PASS PREDICT01068**4 *IONS=',I10/,50X,'PERCENTAGE',11X,'=',F5.1/,50X,'PERCENT WRPT POP-M01069**7 *EAN=',F5.1) 01070**7 WRITE(6,802) TIME 01071**6 802 FORMAT(//,50X,'TOTAL CPU TIME FOR PREDICTIONS=',F12.2,2X,'SECONDS'01072**6 *//,50X,30('*')//) 01073**6 CALL REJECT(DG3,IV,SUM2,REJ,JNC,NLCN,50) 01074**7 WRITE(6,803) REJ,NLCN 01075**8 803 FORMAT(//,30X,F10.2,I5//) 01076**8 IF(NLCN.EQ.0) GO TO 808 01077**7 WRITE(6,807) 01078**7 DO 805 I=1,NLCN 01079**7 IRJ=JNC(I) 01080**7 WRITE(6,806) IRJ,DG3(IRJ) 01081**7 SUMM=SUMM-DG3(IRJ) 01082**8 SUM=SUM-DG3(IRJ)**2 01083**8 805 CONTINUE 01084**7 IV=IV-NLCN 01085**8 DIV=DFLOAT(IV) 01086**8 SUMM1=SUMM/DIV 01087**8 SUM1=(SUM*DIV-SUMM**2)/(DIV*(DIV-1.)) 01088**8 SUM2=DSQRT(SUM1) 01089**8 WRITE(6,752) 01090**8 WRITE(6,810) SUMM1,SUM2 01091**8 752 FORMAT(//,40X,'NEW MEAN AND RMS VALUES:'/,40X,23('=')//) 01092**8 810 FORMAT(40X,'MEAN DIFFERENCE=',F10.5/,40X,'STANDARD DEVIATION=', 01093**8 *F10.5//) 01094**8 806 FORMAT(30X,'OBSERVATION #',I4,2X,'(VALUE:',F7.2,') IS REJECTED') 01095**7 807 FORMAT(//,30X,'REJECTED OBSERVATIONS'/,30X,21('=')//) 01096**7 808 CONTINUE 01097**7 WRITE(6,850) 01098**7 850 FORMAT('1') 01099**7 C 01100**8 C... STUDENTS-T TEST ON MEAN DIFFERENCE: 01101**8 CALL TONMN(SUMM1,SUM2,IV,1) 01102**8 DO 601 I=1,IV 01103**8 Z(I)=DG3(I)/SUM2 01104**8 601 CONTINUE 01105**8 CALL HSTPLT(Z,IV,NCLASS,HVEC,NB) 01106*10 WRITE(6,129) 01107**8 C 01108**8 C... CHI-SQUARE GOODNESS OF FIT TEST: 01109**8 C 01110**8 CALL CHISQ(Z,IV,SEG,ICHI,NKI,CHIVAL,CHI,NKI1) 01111**8 C 01112**8 IAN=IV 01113**8 AM1=0.D0 01114**8 AM2=0.D0 01115**8 AM3=0.D0 01116**8 AM4=0.D0 01117**8 DO 1007 I=1,IAN 01118**8 DELMN=DG3(I)-SUMM1 01119**8 DMN2=DELMN*DELMN 01120**8 DMN3=DMN2*DELMN 01121**8 DMN4=DMN3*DELMN 01122**8 AM1=AM1+DABS(DELMN) 01123**8 AM2=AM2+DMN2 01124**8 AM3=AM3+DMN3 01125**8 AM4=AM4+DMN4 01126**8 1007 CONTINUE 01127**8 RIAM=DFLOAT(IAN) 01128**8 AM0=DSQRT(AM2) 01129**8 AM2=AM2/RIAM 01130**8 AM3=AM3/RIAM 01131**8 AM4=AM4/RIAM 01132**8 ALPHA3=AM3/((DSQRT(AM2))**3) 01133**8 ALPHA4=AM4/(AM2*AM2) 01134**8 A4=AM1/(RIAM*AM0) 01135**8 WRITE(6,1008) ALPHA3,ALPHA4,A4 01136**8 1008 FORMAT(///,40X,'MEASURE OF SKEWNESS=',F12.3/,40X,'MEASURE OF KURTO01137**8 *SIS=',F12.3,5X,F12.3//) 01138**8 C 01139**7 J0=IV/10.+1 01140**7 J1=1 01141**7 J2=J1+9 01142**7 DO 1003 I=1,J0 01143**7 WRITE(6,1004) (SNORM(J),J=J1,J2) 01144**7 J1=J2+1 01145**7 J2=J1+9 01146**7 1003 CONTINUE 01147**7 1004 FORMAT(10X,10F10.3) 01148**7 CALL HSTPLT(SNORM,IV,NCLASS,HVEC,NB) 01149*10 WRITE(6,137) 01150**7 CALL AMNSDR(SNORM,IV,AMEAN,SDEV) 01151**7 WRITE(6,138) AMEAN,SDEV 01152**7 CALL CHISQ(SNORM,IV,SEG,ICHI,NKI,CHIVAL,CHI) 01153**7 C 01154**7 C... TEST ON MEAN: 01155**7 CALL TONMN(AMEAN,SDEV,IV,0) 01156**7 C 01157**7 C... TEST ON VARIANCE: 01158**7 C 01159**7 CHI0=1.D0 01160**7 SDEVSQ=SDEV*SDEV 01161**7 ANVAR=IV*SDEVSQ 01162**7 ANVAR1=ANVAR/CHI1 01163**7 ANVAR2=ANVAR/CHI2 01164**7 WRITE(6,124) 01165**7 WRITE(6,125) ANVAR1,CHI0,ANVAR2 01166**7 IF((CHI0.GT.ANVAR1).AND.(CHI0.LT.ANVAR2)) GO TO 1005 01167**7 WRITE(6,126) 01168**7 GO TO 1006 01169**7 1005 WRITE(6,127) 01170**7 1006 CONTINUE 01171**7 124 FORMAT(//,40X,'TEST ON VARIANCE'/,40X,14('=')//) 01172**7 125 FORMAT(40X,F12.5,2X,'<',F12.5,2X,'<',F12.5) 01173**7 126 FORMAT(//,40X,'** TEST FAILS **'//) 01174**7 127 FORMAT(//,40X,'** TEST PASSES **'//) 01175**7 129 FORMAT(//,40X,'HISTOGRAM PLOT OF DIFFERENCES'/,40X,29('=')//) 01176**7 130 FORMAT(//,30X,F5.2,4X,F5.2,6X,I2,12X,I2,12X,F5.1,8X,F5.2) 01177**7 131 FORMAT(//,39X,F5.2,6X,I2,12X,I2,12X,F5.1,8X,F5.2) 01178**7 135 FORMAT(//,91X,F5.2) 01179**7 136 FORMAT(//,37X,'SEGMENT',3X,'OBSERVED',6X,'EXPECTED',8X,'(O-E)',6X,01180**7 *'(O-E)**2/E') 01181**7 137 FORMAT(//,40X,'HISTOGRAM PLOT OF STANDARDIZED DIFFERENCES'/,40X,4201182**7 *('=')//) 01183**7 138 FORMAT(40X,'MEAN VALUE=',F10.2/,40X,'STD. DEV=',F10.2//) 01184**7 GO TO 99 01185**4 950 CONTINUE 01186**6 IF(IFILE.EQ.1) GO TO 954 01187**6 IF(MCV.EQ.0) GO TO 99 01188**6 WRITE(6,951) 01189**6 951 FORMAT(//,25X,'MEAN COVARIANCE',4X,'STD. DEV.'/) 01190**6 AMCV=DFLOAT(MCV) 01191**6 DO 952 I=1,12 01192**6 CVRSUM(I)=CVRSUM(I)/AMCV 01193**6 VCVRSM(I)=VCVRSM(I)/AMCV 01194**6 VCVRSM(I)=DSQRT(VCVRSM(I)) 01195**6 WRITE(6,953) I,CVRSUM(I),VCVRSM(I) 01196**6 952 CONTINUE 01197**6 953 FORMAT(20X,I5,2D15.5) 01198**6 GO TO 99 01199**6 954 WRITE(6,955) IFL 01200**6 955 FORMAT(//,20X,'TOTAL # OF RECORDS IN FILE=',I10//) 01201**6 GO TO 99 01202**6 C**************R********************************************************01203 C 01204 960 WRITE(6,961) 01205 961 FORMAT(///,10X,'ERROR MESSAGE : THE GRAMS MATRIX OF NORMAL EQUATIO01206 1NS (Q) FOR ESTIMATING THE FREE-AIR GRAVITY'/,15X,'ANOMALY SURFACE 01207 2TURNED-OUT TO BE SINGULAR, SOLUTION IS IMPOSSIBLE'/,10X,'PLEASE CH01208 3ECK YOUR SCALING OF THE (X,Y) LOCAL COORDINATES OF DATA POINTS AND01209 4/OR'/15X,'THE SCALING OF THE MATRIX-Q ITSELF') 01210 GO TO 9999 01211 970 WRITE(6,971) 01212 971 FORMAT(///,10X,'ERROR MESSAGE : THE GRAMS MATRIX OF NORMAL EQUATIO01213 1NS (Q) FOR ESTIMATING THE HEIGHT (TOPOGRAPHY) '/,15X,'SURFACE 01214 2TURNED-OUT TO BE SINGULAR, SOLUTION IS IMPOSSIBLE'/,10X,'PLEASE CH01215 3ECK YOUR SCALING OF THE (X,Y) LOCAL COORDINATES OF DATA POINTS AND01216 4/OR'/15X,'THE SCALING OF THE MATRIX-Q ITSELF') 01217 GO TO 9999 01218 980 WRITE(6,981) 01219 981 FORMAT(///,10X,'ERROR MESSAGE: THE LEAST SQUARES ESTIMATION FOR TH01220 1E POLYNOMIAL COEFFICIENTS IS IMPOSSIBLE'/,25X,'BECAUSE THE NUMBER 01221 2OF AVAILABLE DATA POINTS IS LESS THAN THE NUMBER OF UNKNOWN COEFFI01222 2CIENTS (-VE DEGREES OF FREEDOM)') 01223 GO TO 9999 01224 990 CONTINUE 01225**4 992 WRITE(6,991) 01226**4 991 FORMAT(///,10X,'ERROR MESSAGE: THE CHOSEN SIZE OF THE AREA AROUND 01227 1THE ABOVE POINT OF INTEREST'/,25X,'DOES NOT HAVE A RELIABLE DISTRI01228 2BUTION OF AT LEAST ONE DATA POINT IN EACH QUADRANT') 01229 777 FORMAT('1') 01230 C 01231 9999 GO TO 999 01232 C 01233 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*01234 C 01235 99 STOP 01236 END 01237 C 01238 C 01239 C***********************************************************************01240**3 SUBROUTINE INPOLY(XPOLY,YPOLY,NP1,XI,YI,IN) 01241**3 C 01242**3 C THIS SUBROUTINE RETURNS A VALUE OF IN=1 IF POINT (XI,YI) IS 01243**3 C WITHIN THE POLYGON SPECIFIED BY XPOLY AND YPOLY. A VALUE OF 01244**3 C IN=0 IS RETURNED IF THE POINT IS NOT WITHIN THE POLYGON. 01245**3 C 01246**3 C NP1 = NUMBER OF POLYGON POINTS + 1 (THE POLYGON MUST BE CLOSED , 01247**3 C I.E. THE FIRST AND LAST POINT MUST BE THE SAME IN THE DATA 01248**3 C VECTORS XPOLY AND YPOLY) 01249**3 C 01250**3 C********************************************************************* 01251**9 C 01252**9 IMPLICIT REAL*8(A-H,O-Z) 01253**3 DIMENSION XPOLY(NP1),YPOLY(NP1) 01254**3 IN=0 01255**3 N12 = NP1 - 1 01256**5 DO 625 I=1,N12 01257**5 IF(YI.LE.YPOLY(I))GO TO 623 01258**3 622 IF(YI.GT.YPOLY(I+1)) GO TO 625 01259**3 624 IF((XI-XPOLY(I) -(YI-YPOLY(I))*(XPOLY(I+1)-XPOLY(I))/(YPOLY(I+1) 01260**3 & -YPOLY(I))).LT.0.0) IN=1-IN 01261**3 GO TO 625 01262**3 623 IF(YI.GT.YPOLY(I+1)) GO TO 624 01263**3 625 CONTINUE 01264**3 RETURN 01265**3 END 01266**3 C***********************************************************************01267 C* 01268 C* SUBROUTINE 'A P P R O X' :- (BY: M.M. NASSAR, JANUARY 1975.)01269 C* ------------------------ 01270 C* 01271 C* THIS SUBROUTINE HAS BEEN DESIGNED AS A 01272 C* GENERAL SUBROUTINE FOR PERFORMING A LEAST-SQUARES APPROXIMATION TO AN01273 C* DATA SET, USING A MIXED-ALGEBRAIC POLYNOMIAL OF ANY ARBITRARY ORDER 01274 C* IN TWO DIMENSIONS, WHICH IMPLIES THAT THIS SUBROUTINE HAS THE ADVANTA01275 C* OF EXECUTION-TIME-DIMENSIONING. 01276 C* 01277 C* FOR DETAILS CONCERNING THE MATHEMATICAL MODELS AND RELATED DEVELOPM- 01278 C* ENTS, SEE : NASSAR AND VANICEK (1975).'LEVELLING AND GRAVITY'-TECHN- 01279 C* ICAL REPORT # 33, SURVEYING ENGINEERING DEPT., U.N.B., FREDERICTON. 01280 C* 01281 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*01282 C* 01283 C* THE INDIVIDUAL ELEMENTS OF THE DESIGN MATRIX-A ARE GENERATED (VIA THE01284 C* FUNCTION-SUBPROGRAM 'A') ONE AT A TIME WHENEVER NEEDED IN THE COMPUTA01285 C* IONS WITHOUT THE NECESSITY OF STORING THEM, TO SAVE COMPUTER STORAGE.01286 C* (SEE THE LISTINGS OF THE FUNCTION-SUBPROGRAM 'A' CONCERNING THE REQUI01287 C* COMMON BLOCK WITH THE MAIN PROGRAM). 01288 C* 01289 C* THE INVERSION OF THE NORMAL EQUATIONS MATRIX IS DONE VIA THE SUBROUTI01290 C* 'CHOLD', WHICH IS BASED ON THE CHOLESKY DECOMPOSITION METHOD. IN THI01291 C* RESPECT, APPROX CONTAINS THE POSSIBILITY OF SCALING THE NORMAL EQUATI01292 C* MATRIX (AS A WHOLE) BY THE DESIRED SCALE FACTOR, IF REQUIRED BEFORE 01293 C* AND AFTER THE INVERSION PROCESS. 01294 C* 01295 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*01296 C* 01297 C* USAGE:- CALL APPROX(N,M,F,W,SCALE,C,Q,EVF,U,DETQ) 01298 C* 01299 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*01300 C* 01301 C* THE REQUIRED INPUT DATA :- (VIA THE CALLING ROUTINE) 01302 C* ----------------------- 01303 C 01304 C* N : TOTAL NUMBER OF AVAILABLE DATA POINTS. 01305 C* M : NUMBER OF UNKNOWN COEFFICIENTS OF THE APPROXIMATING POLYNOMIA01306 C* F : VECTOR OF OBSERVED DATA AT THE N DISCRETE POINTS. 01307 C* W : VECTOR OF WEIGHTS FOR THE N OBSEVED QUANTITIES. 01308 C* SCALE : DESIRED SCALE FACTOR FOR THE NORMAL EQUATIONS MATRIX. 01309 C* 01310 C* IT SHOULD BE NOTED HERE THAT THE VECTORS OF X & Y LOCAL COORDINATES ,01311 C* AS WELL AS THE SPECIFIED ORDER OF THE POLYNOMIAL ARE PART OF THE INPU01312 C* DATA ,HOWEVER THEY MUST BE DECLARED AS A COMMON BLOCK BETWEEN THE MAI01313 C* PROGRAM AND THE FUNCTION 'A'. 01314 C* 01315 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*01316 C* 01317 C* THE MAIN OUTPUT ARE :- 01318 C* ------------------- 01319 C* 01320 C* C : VECTOR OF ESTIMATED COEFFICIENTS OF THE APPROXIMATING POLYNOMI01321 C* Q : ESTIMATED VARIANCE-COVARIANCE MATRIX OF THE COEFFICIENTS 'C'. 01322 C* EVF : THE A POSTERIORI VARIANCE FACTOR. 01323 C* U : VECTOR OF ABSOLUTE TERMS OF THE NORMAL EQUATIONS. 01324 C* DETQ : DETERMINANT OF THE NORMAL EQUATIONS MATRIX. 01325 C* 01326 C* NOTE HERE THAT THE REASON FOR INCLUDING THE VECTOR-U IN THE ARGUMENT 01327 C* LIST IS SOLELY TO ACHIEVE THE ADVANTAGE OF ITS EXECUTION TIME DIMENSI01328 C* 01329 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*01330 C* 01331 C* ON THE OTHER HAND THE PRINTED OUTPUT FROM THE LEAST-SQUARES APPROXIMA01332 C* ION TECHNIQUE ,BY THIS SUBROUTINE, INCLUDES :- 01333 C* 01334 C* 1-THE NORMAL EQUATIONS MATRIX. 01335 C* 2-THE VECTOR-U OF ABSOLUTE TERMS OF THE NORMAL EQUATIONS. 01336 C* 3-THE DETERMINANT OF THE NORMAL EQUATIONS MATRIX. 01337 C* 4-THE VECTOR OF ESTIMATED COEFFICIENTS OF THE APPROXIMATING POLYNOMIA01338 C* 5-THE VECTOR-V OF RESIDUALS. 01339 C* 6-THE NUMBER OF DEGREES OF FREEDOM. 01340 C* 7-THE ESTIMATED SQUARE-SUM OF THE WEIGHTED RESIDUALS. 01341 C* 8-THE A POSTERIORI (ESTIMATED) VARIANCE FACTOR. 01342 C* 01343 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*01344 C* 01345 C* FINALLY, IT SHOULD BE NOTED THAT THE ORDER OF THE APPROXIMATING POLYN01346 C* IAL IS LEFT ARBITRARY, HOWEVER, THE NUMBER OF DATA POINTS IS ALLOWED 01347 C* ONLY UP TO 2000, AND IF EXCEEDED THE CORRESPONDING DIMENSION STATEMEN01348 C* IN THE MAIN PROGRAM AND THE FUNCTION 'A' MUST BE MODIFIED. 01349 C* 01350 C***********************************************************************01351 C 01352 SUBROUTINE APPROX(N,M,F,W,SCALE,C,Q,EVF,U,DETQ,IDEXP,IPRT,VIV, 01353**6 *VVIV,NCVR,NCOMB,INDEX) 01354**6 IMPLICIT REAL*8(A-H,O-Z) 01355 DIMENSION F(500),W(500),VIV(50) 01356**6 DIMENSION Q(25,25),U(25),C(25) 01357**4 DIMENSION QA(500,25),AQ(25),AQAT(500),VVIV(500) 01358**6 C 01359 C INITIALIZE ZERO-VALUES FOR THE ELEMENTS OF SOME RESULTANT VECTORS & MA01360 C 01361 DO 15 I=1,M 01362 DO 15 J=1,M 01363 15 Q(I,J)=0.0D 0 01364 DO 25 I=1,M 01365 U(I)=0.0D 0 01366 25 C(I)=0.0D 0 01367 C 01368 C***********************************************************************01369 888 FORMAT('1') 01370 C***********************************************************************01371 C 01372 C FORMULATE THE NORMAL-EQUATIONS MATRIX (WHICH IS KNOWN AS 'GRAM'S-MATRI01373 C IN THE LEAST-SQUARES TERMINOLOGY, AND WHOSE ELEMENTS ARE THE COMBINATI01374 C OF THE SCALAR-PRODUCT OF THE BASE-FUNCTIONS). THIS MATRIX 'Q' IS GIVEN01375 C Q = AT * P * A , 01376 C WHERE 'A' IS THE DESIGN MATRIX, AND 'P' IS THE WEIGHT MATRIX. 'A' IS K01377 C AS VANDERMONDE'S-MATRIX IN THE LEAST-SQUARES APPROXIMATION TERMINOLOGY01378 C THE ELEMENTS IN EACH ROW OF 'A' PERTAIN TO ONLY ONE OBSERVED DATA POIN01379 C AND ARE SIMPLY THE FUNCTIONAL VALUES OF THE M CHOSEN BASE-FUNCTIONS I01380 C THE MIXED (ALGEBRAIC, NON-LINEAR, THIRD-ORDER, IN TWO-DIMENSIONS) 01381 C APPROXIMATING POLYNOMIAL. I.E. ONE ROW IN A-MATRIX FOR EACH DATA POIN01382 C 01383 C IT SHOULD BE NOTED HERE THAT THE INDIVIDUAL ELEMENTS 'A(ROW,COL.)' OF 01384 C A-MATRIX WHICH ARE NEEDED FOR COMPUTING THE Q-MATRIX, WILL BE EVALUATE01385 C ONE AT A TIME THROUGH THE FUNCTION-SUBPROGRAM ''A'' WHICH IS AUTOMATIC01386 C CALLED AS SOON AS THE EXECUTION HITS ANY EXPRESSION INVOLVING ANY ELEM01387 C OF THE DESIGN MATRIX-A. 01388 C 01389 C THE 'P' MATRIX IS A DIAGONAL WEIGHT MATRIX, CONTAINING THE WEIGHT FUNC01390 C W'S ON ITS MAIN DIAGONAL. 01391 C 01392 C... THE A-MATRIX: 01393**6 DO 16 I=1,N 01394**6 DO 16 J=1,M 01395**6 QA(I,J)=A(I,J) 01396**6 16 CONTINUE 01397**6 C... NORMAL-MATRIX: 01398**6 EPS=2.3D-16 01399**6 DO 17 I=1,M 01400**6 DO 17J=I,M 01401**6 DO 17 K=1,N 01402**6 QB=QA(K,I) 01403**6 QC=QA(K,J) 01404**6 Q(I,J)=Q(I,J)+QB*QC*W(K) 01405**6 17 CONTINUE 01406**6 C... FILL-IN REMAINING ELEMENTS OF NORMAL-MATRIX: 01407**6 DO 18 I=2,M 01408**6 IMN1=I-1 01409**6 DO 18 J=1,IMN1 01410**6 Q(I,J)=Q(J,I) 01411**6 18 CONTINUE 01412**6 C 01413 C PRINT-OUT THE NORMAL-EQUATIONS MATRIX 01414 C 01415 IF(IPRT.EQ.0) GO TO 121 01416**5 WRITE(6,81) M,M 01417 81 FORMAT(/,5X,'THE NORMAL EQUATIONS MATRIX-Q'/,10X,'DIMENSION : (',I01418 13,' ,',I3,' )') 01419 WRITE(6,82) 01420 82 FORMAT(/,10X,'1',12X,'2',12X,'3',12X,'4',12X,'5',12X,'6',12X,'7',101421 12X,'8'/) 01422 MHALF=M/2 01423 MHALF1=MHALF+1 01424 DO 110 I=1,M 01425 110 WRITE(6,83) I,(Q(I,J),J=1,MHALF) 01426 83 FORMAT(I4,8D13.4) 01427 WRITE(6,84) 01428 84 FORMAT(/,10X,'9',11X,'10',11X,'11',11X,'12',11X,'13',11X,'14',11X,01429 1'15',11X,'16'/) 01430 DO 120 I=1,M 01431 120 WRITE(6,83) I,(Q(I,J),J=MHALF1,M) 01432 121 CONTINUE 01433**4 C***********************************************************************01434 C 01435 C FORMULATE THE VECTOR-U OF ABSOLUTE TERMS OF THE NORMAL EQUATIONS (WHOS01436 C ELEMENTS ARE THE SCALAR PRODUCTS OF THE BASE-FUNCTIONS AND THE OBSERVE01437 C FUNCTION (F-VECTOR), I.E. U=AT*P*F WHERE A & P ARE AS DEF01438 C ABOVE, AND 'F' IS THE VECTOR OF OBSERVED QUANTITIES -TO BE USED IN COM01439 C THE BEST-FITTING POLYNOMIAL. 01440 C 01441 DO 30 I=1,M 01442 U(I)=0.0D 0 01443 DO 30 J=1,N 01444 U(I)=U(I)+A(J,I)*W(J)*F(J) 01445 30 CONTINUE 01446 IF(IPRT.EQ.0) GO TO 123 01447**5 C 01448 C PRINT-OUT THE VECTOR 'U' 01449 C 01450 WRITE(6,888) 01451 WRITE(6,85) M 01452 85 FORMAT(//,5X,'THE VECTOR-U OF ABSOLUTE TERMS OF THE NORMAL EQUATIO01453 1NS'/,10X,'DIMENSIONS : (',I3,' , 1 )') 01454 WRITE(6,86) 01455 86 FORMAT(/,12X,'1'/) 01456 DO 130 I=1,M 01457 130 WRITE(6,87) I,U(I) 01458 87 FORMAT(I4,D15.5) 01459 123 CONTINUE 01460**4 IF(SCALE.EQ.1.0D 0) GO TO 215 01461**4 C***********************************************************************01462 C 01463 C SCALE THE NORMAL-EQUATIONS MATRIX-Q, BEFORE INVERSION, TO OVERCOME THE01464 C OVERFLOW PROBLEM WHICH IS LIKELY TO OCCUR DURING THE INVERSION PROCESS01465 C 01466 DO 5 I=1,M 01467 DO 5 J=1,M 01468 5 Q(I,J)=Q(I,J)*SCALE 01469 215 CONTINUE 01470 C 01471 C CALL THE SUBROUTINE 'CHOLD' TO INVERT THE NORMAL-EQUATIONS MATRIX 01472 C 'Q' (I.E. GRAM'S-MATRIX), UNDER THE SAME NAME 'Q' AND IN THE SAME LOCA01473 C THIS INVERSE CAN BE REFERRED TO AS THE WEIGHT-COEFFICIENT MATRIX (IN T01474 C LEAST-SQUARES SENSE) OF THE COEFFICIENTS OF THE BEST-FITTING POLYNOMIA01475 C ***FOR DETAILS REGARDING THE DESCRIPTION OF THE INPUT/OUTPUT DATA FOR 01476 C 'CHOLD', SEE ITS DOCUMENTATION.*** 01477 C 01478 C 01479**3 CALL SPIN(Q,M,25,DETQ,IDEXP,INDEX) 01480**6 C 01481**6 IF(INDEX.EQ.1) GO TO 99 01482**6 IF(SCALE.EQ.1.0D 0) GO TO 225 01483 C 01484 C RE-SCALE THE INVERSE MATIX-Q (OF THE NORMAL EQUATIONS), AFTER THE INVE01485 C PROCESS, IN ORDER TO GET THE TRUE (CORRECT) WEIGHT-COEFFICIENT MATRIX 01486 C OF THE ESTIMATED POLYNOMIAL-COEFFICIENTS. 01487 C 01488 DO 210 I=1,M 01489 DO 210 J=1,M 01490 210 Q(I,J)=Q(I,J)*SCALE 01491 225 CONTINUE 01492 C***********************************************************************01493 C 01494 C COMPUTE THE ESTIMATED VALUES OF THE COEFFICIENTS (C'S) OF THE BEST-FIT01495 C POLYNOMIAL. C=Q(INVERSE)*U 01496 C 01497 DO 40 I=1,M 01498 C(I)=0.0D 0 01499 DO 40 J=1,M 01500 C(I)=C(I)+Q(I,J)*U(J) 01501 40 CONTINUE 01502 IF(IPRT.EQ.0) GO TO 141 01503**5 C 01504 C PRINT-OUT THE VECTOR-C OF ESTIMATED COEFFICIENTS 01505 C 01506 WRITE(6,88) M 01507 88 FORMAT(//,5X,'THE VECTOR-C OF ESTIMATED POLYNOMIAL-COEFFICIENTS'/,01508 110X,'DIMENSIONS : (',I3,' , 1 )') 01509 WRITE(6,86) 01510 DO 140 I=1,M 01511 140 WRITE(6,87) I,C(I) 01512 C***********************************************************************01513 C 01514 C COMPUTE THE SQUARE-SUM OF THE WEIGHTED RESIDUALS VTPV=VT*P*V, WHER01515 C IS AS DEFINED EARLIER, AND 'V' IS THE VECTOR OF RESIDUALS (OR DISCREP01516 C WHICH ARE GIVEN BY THE DIFFERENCE BETWEEN THE OBSERVED VALUES (VECTOR-01517 C AND THE CORRESPONDING PREDICTED VALUES 'A*C' (AS COMPUTED FROM THE BES01518 C FITTING POLYNOMIAL), I.E. V=A*C-F, IT SHOULD BE NOTED HERE T01519 C THESE RESIDUALS HAVE THE SAME PHYSICAL UNITS AS THE OBSERVED VALUES 'F01520 C 01521 C IT SHOULD BE REMEMBERED HERE AGAIN THAT THE INDIVIDUAL ELEMENTS OF THE01522 C DESIGN MATRIX-A ARE EVALUATED VIA THE FUNCTION-SUBPROGRAM ''A'' 01523 C 01524 WRITE(6,888) 01525 WRITE(6,122) N 01526 122 FORMAT(//,5X,'THE VECTOR-V OF RESIDUALS'/,5X,'DIMENSIONS : (',I4,'01527 1 , 1 )') 01528 WRITE(6,86) 01529 141 CONTINUE 01530**4 VTPV=0.D0 01531**5 VSUM=0.D0 01532**5 DO 60 I=1,N 01533 AC=0.0D 0 01534 DO 50 J=1,M 01535 50 AC=AC+A(I,J)*C(J) 01536 C 01537 C COMPUTE AND PRINT-OUT THE VECTOR OF RESIDUALS -DISCREPANCIES BETWEEN 01538 C THE OBSERVED QUANTITIES AND THEIR CORRESPONDING PREDICTED VALUES. 01539 C 01540 V=F(I)-AC 01541 VSUM=VSUM+V 01542**5 VIV(I)=V 01543**5 VTPV=VTPV+V*W(I)*V 01544 60 CONTINUE 01545 C 01546 C COMPUTE THE NUMBER OF DEGREES OF FREEDOM, WHICH IS THE DIFFERENCE BETW01547 C THE TOTAL NUMBER OF USED OBSERVATION POINTS AND THE NUMBER OF UNKNOWN 01548 C POLYNOMIAL-COEFFICIENTS 01549 C 01550 NUMDF=N-M 01551 C 01552 C COMPUTE THE ESTIMATED (IN THE LEAST-SQUARES SENCE) VARIANCE FACTOR 'EV01553 C -DENOTED BY SIGMA-ZERO-SQUARE-CAP. IT SHOULD BE NOTED HERE, THAT THE 01554 C IS UNITLESS, DUE TO OUR STIPULATION OF THE WEIGHTS (INVERSE OF VARIANC01555 C 01556 EVF=VTPV/FLOAT(NUMDF) 01557 C 01558 C COMPUTE THE ESTIMATED VARIANCE-COVARIANCE MATRIX 'EC' (SIGMA-C) OF THE01559 C ESTIMATED POLYNOMIAL-COEFFICIENTS (C'S), SIMPLY BY MULTIPLYING THE 01560 C ESTIMATED VARIANCE-FACTOR BY THE INVERSE OF THE NOMAL-EQUATIONS MATRIX01561 C I.E. EC=EVF*Q(INVERSE) 01562 C 01563 IF((NCVR.EQ.1).OR.(NCOMB.EQ.1)) GO TO 98 01564**6 IF(IPRT.EQ.0) GO TO 96 01565**5 IF(IPRT.EQ.1) GO TO 75 01566**6 98 CONTINUE 01567**6 DO 71 I=1,N 01568**6 AQAT(I)=0.D0 01569**6 DO 72 J=1,M 01570**6 AQ(J)=0.D0 01571**6 DO 72 K=1,M 01572**6 QD=QA(I,K) 01573**6 QE=Q(K,J) 01574**6 AQ(J)=AQ(J)+QD*QE 01575**6 72 CONTINUE 01576**6 DO 71 KK=1,M 01577**6 AQAT(I)=AQAT(I)+AQ(KK)*QA(I,KK) 01578**6 71 CONTINUE 01579**6 C... COVARIANCE MATRIX OF RESIDUALS: 01580**6 DO 73 I=1,N 01581**6 VVIV(I)=1.D0/W(I)-AQAT(I) 01582**6 IF(IPRT.EQ.1) WRITE(6,74) I,VIV(I),VVIV(I) 01583**6 73 CONTINUE 01584**6 74 FORMAT(20X,I5,2D15.5) 01585**6 IF(IPRT.EQ.0) GO TO 97 01586**6 75 CONTINUE 01587**6 C 01588**4 DO 70 I=1,M 01589 DO 70 J=1,M 01590 Q(I,J)=EVF*Q(I,J) 01591 70 CONTINUE 01592 C 01593 C PRINT-OUT THE 'NUMDF' , 'VTPV' AND 'EVF' SCALARS, AND THE 'EC' MATRIX 01594 C 01595 WRITE(6,14) VSUM 01596**6 14 FORMAT(//,5X,'SUM OF RESIDUALS=',D15.5//) 01597**6 WRITE(6,89) NUMDF,VTPV,EVF 01598 89 FORMAT(/,/,5X,'THE NUMBER OF DEGREES OF FREEDOM =',I5//,5X,'THE ES01599 1TIMATED SQUARE-SUM OF THE WEIGHTED RESIDUALS =',D15.5//,5X,'THE ES01600 2TIMATED VARIANCE FACTOR =',D15.5/) 01601 WRITE(6,888) 01602 WRITE(6,91) M,M 01603 91 FORMAT(/,5X,'THE ESTIMATED VARIANCE-COVARIANCE MATRIX-EC OF THE ES01604 1TIMATED POLYNOMIAL-COEFFICIENTS'/,15X,'DIMENSIONS : (',I3,' ,',I3,01605 2' )') 01606 WRITE(6,82) 01607 DO 150 I=1,M 01608 150 WRITE(6,83) I,(Q(I,J),J=1,MHALF) 01609 WRITE(6,84) 01610 DO 160 I=1,M 01611 160 WRITE(6,83) I,(Q(I,J),J=MHALF1,M) 01612 GO TO 97 01613**4 96 Q(1,1)=EVF*Q(1,1) 01614**4 97 RETURN 01615 99 CONTINUE 01616**6 WRITE(6,100) 01617**6 100 FORMAT(//,20X,'NEGATIVE OR ZERO ELEMENT ON THE DIAGONAL OF NORMAL-01618**6 *MATRIX: CHECK CAUSE:'//) 01619**6 RETURN 01620**6 END 01621 C***********************************************************************01622 C* 01623 C* FUNCTION-SUBPROGRAM ' A ' :- 01624 C* ------------------------- 01625 C* 01626 C* THIS FUNCTION-SUBPROGRAM IS DESIGNED AS 01627 C* GENERAL FUNCTION, TO COMPUTE (IN DOUBLE PRECISION) SINGLE ELEMENTS 01628 C* A(I,J) OF THE DESIGN MATRIX-A (WHICH IS ALSO KNOWN AS VANDERMONDE'S 01629 C* MATRIX), NEEDED IN THE FORMATION OF THE NORMAL EQUATIONS SYSTEM IN TH01630 C* LEAST-SQUARES APPROXIMATION TECHNIQUE- USING A TWO-DIMENSIONAL MIXED-01631 C* ALGEBRAIC POLYNOMIAL OF ANY DESIRED ORDER (IDEG), WHERE I IS THE ROW 01632 C* POSITION AND J IS THE COLUMN POSITION OF THE SOUGHT ELEMENT A(I,J). 01633 C* 01634 C* NOTE THAT THE MAIN IDEA BEHIND THIS FUNCTION IS TO SAVE A LARGE AMOUN01635 C* OF COMPUTER STORAGE, DUE TO THE FACT THAT MATRIX-A DOES NOT NEED TO B01636 C* GENERATED IN FULL, SAY WITH DIMENSIONS (N BY M), AND INSTEAD WE COMPU01637 C* ONLY ONE OF ITS ELEMENTS AT A TIME WHEN REQUIRED. 01638 C* 01639 C* FOR MORE DETAILS CONCERNING THE MATHEMATICAL BASIS OF THIS FUNCTION, 01640 C* SEE: 'INTRODUCTION TO APPROXIMATION THEORY' BY: E.W. CHENEY, 1966, 01641 C* MCGRAW-HILL BOOK COMPANY, NEW YORK. 01642 C* 01643 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*01644 C* 01645 C* USAGE:- THIS FUNCTION WILL BE AUTOMATICALLY CALLED (TO EVALUATE A(I,J01646 C* AS SOON AS THE EXECUTION IN THE CALLING ROUTINE HITS ANY 01647 C* EXPRESSION THAT INVOLVES ANY INDIVIDUAL ELEMENT A(I,J) OF A. 01648 C* 01649 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*01650 C* 01651 C* THE REQUIRED INPUT DATA ARE THE VECTORS OF X & Y LOCAL ORTHOGONAL 01652 C* COORDINATES OF ALL THE DATA POINTS AVAILABLE IN THE INTERPOLATION 01653 C* AREA (E. G. ONE-DEGREE-SQUARE), AND THE DESIRED ORDER OF THE APPROXIM01654 C* ATING POLYNOMIAL.. THIS DATA MUST BE SPECIFIED AS A LABELLED COMMON 01655 C* BLOCK BETWEEN THIS FUNCTION-SUBPROGRAM AND THE MAIN ROUTINE THAT 01656 C* COMPUTES THE X,Y COORDINATES.. 01657 C* 01658 C* THE OUTPUT WILL BE THE NUMERICAL VALUE OF A(I,J) TRANSFERRED TO THE 01659 C* RESPECTIVE EXPRESSION IN THE CALLING ROUTINE . 01660 C* 01661 C* IT SHOULD BE NOTED THAT ALTHOUGH THERE ARE NO LIMITATIONS PUT ON THE 01662 C* VALUE OF THE ORDER OF THE APPROXIMATING POLYNOMIAL, THE NUMBER OF DAT01663 C* POINTS IS NOT TO EXCEED 2000 , OTHERWISE THE CORRESPONDING DIMENSION 01664 C* STATEMENTS SHOULD BE MODIFIED,.. 01665 C* 01666 C***********************************************************************01667 C 01668 DOUBLE PRECISION FUNCTION A(I,J) 01669 DOUBLE PRECISION X,Y,XA,YA 01670 DIMENSION X(500),Y(500) 01671**6 C 01672 C DECLARE EXPLICITLY THE LABELLED COMMON BLOCK BETWEEN THE MAIN-PROGRAM 01673 C AND THIS FUNCTION-SUBPROGRAM 01674 C THIS COMMON BLOCK SHOULD CONTAIN THE VECTORS OF X & Y PLANE LOCAL 01675 C COORDINATES, AND THE DEGREE OF THE APPROXIMATING MIXED-ALGEBRAIC POLYN01676 C 01677 COMMON /XYMAIN/ X,Y,IDEG 01678 NA=IDEG+1 01679**6 K=(J-1)/NA 01680**6 IL2=J-K-(NA-1)*K-1 01681**6 XA=X(I) 01682 YA=Y(I) 01683 C 01684 C CHECK AND RESPOND TO OVERCOME THE PROBLEM OF TERMINATING THE PROGRAM D01685 C TO THE OCCURENCE OF A BASE EQUALS ZERO WITH AN EXPONENT ALSO EQUALS ZE01686 C WHICH WILL RESULT IN AN ERROR-MESSAGE AND STOPPING THE PROGRAM-EXECUTI01687 C 01688 IF(XA.EQ.0.0D 0.AND.K.EQ.0) XA=1.0D 0 01689 IF(YA.EQ.0.0D 0.AND.IL2.EQ.0) YA=1.0D 0 01690**5 A=XA**K*YA**IL2 01691**5 RETURN 01692 END 01693 C 01694 C 01695 C***********************************************************************01696**2 SUBROUTINE SPIN(Q,N,MM,DET,IDEXP,INDEX) 01697**6 C 01698**3 C SUBROUTINE SPIN IS A MATRIX INVERSION ROUTINE FOR SYMMETRIC 01699**3 C POSITIVE-DEFINITE MATRICES. THE MATRIX INVERTED IS THE UPPER 01700**3 C N BY N PORTION OF THE MATRIX Q WHICH IS DIMENSIONED MM BY MM 01701**3 C IN THE CALLING ROUTINE. 01702**3 C 01703**3 C INPUT: 01704**3 C Q - THE MATRIX DIMENSIONED MM BY MM WHICH CONTAINS THE 01705**3 C MATRIX TO BE INVERTED. 01706**3 C 01707**3 C N - THE DIMENSION OF THE ACTUAL PART( UPPER LEFT CORNER) 01708**3 C OF Q WHICH IS TO BE INVERTED. ( N MAY BE EQUAL BUT MUST 01709**3 C NOT BE LARGER THAN MM) . 01710**3 C MM- DIMENSIONED SIZE OF Q IN THE CALLING ROUTINE. 01711**3 C 01712**3 C 01713**3 C OUTPUT: 01714**3 C 01715**3 C Q - THE UPPER LEFT N BY N PORTION CONTAINS THE INVERSE OF 01716**3 C THE INPUT UPPER LEFT N BY N PORTION. 01717**3 C 01718**3 C DET - THE NON-EXPONENT PORTION OF THE DETERMINANT OF THE 01719**3 C INPUT N BY N (UPPER LEFT PORTION OF Q) MATRIX. SEE 01720**3 C IDEXP BELOW. 01721**3 C 01722**3 C IDEXP - THE EXPONENT (OF 10) PART OF THE DETERMINANT DESCRIBED 01723**3 C UNDER DET ABOVE. THUS THE DETERMINANT IS RETURNED IN 01724**3 C TWO PARTS CORRESPONDING TO 01725**3 C DETERMINANT = DET * 10 ** IDEXP . 01726**3 C THIS IS DONE TO AVOID UNDER OR OVERFLOW IN THE 01727**3 C COMPUTATION OF THE DETERMINANT. TO PRINT THE DETERM- 01728**3 C INANT THE USER SHOULD PRINT BOTH NUMBERS AS FOLLOWS; 01729**3 C (FOR EXAMPLE) 01730**3 C PRINT 10,DET,IDEXP 01731**3 C 10 FORMAT(' ','DETERMINANT= ',F7.4,'D',I4) 01732**3 C 01733**3 C R.R. STEEVES 01734**3 C SEPT., 1979 01735**3 C 01736**3 C 01737**3 C********************************************************************* 01738**9 C 01739**9 REAL*8 Q(MM,MM),DET,SUM,DSQRT,DABS,RPART,APART 01740**3 INDEX=0 01741**6 DET=0.D0 01742**3 DO 4 J=1,N 01743**3 DO 4 I=1,J 01744**3 IF(I.EQ.1) GO TO 2 01745**3 M=I-1 01746**3 SUM=0.0D0 01747**3 DO 1 K=1,M 01748**3 1 SUM=SUM+Q(K,I)*Q(K,J) 01749**3 Q(I,J)=Q(I,J)-SUM 01750**3 2 IF(I.EQ.J) GO TO 3 01751**3 Q(I,J)=Q(I,J)/Q(I,I) 01752**3 GO TO 4 01753**3 3 CONTINUE 01754**3 IF(Q(I,I).LE.0.D0) GO TO 11 01755**6 DET=DET+DLOG10(Q(I,I)) 01756**3 Q(I,I)=DSQRT(Q(I,I)) 01757**3 4 CONTINUE 01758**3 IDEXP=DET 01759**3 RPART=DET-IDEXP 01760**3 APART=DABS(RPART) 01761**3 IF(APART.LT.1.D-20) DET=1.D0 01762**3 IF(APART.LT.1.D-20) GO TO 10 01763**3 DET=10**RPART 01764**3 10 CONTINUE 01765**3 DO 7 J=1,N 01766**3 DO 7 I=1,J 01767**3 IF(I.LT.J) GO TO 5 01768**3 Q(J,J)=1.0D0/Q(J,J) 01769**3 GO TO 7 01770**3 5 SUM=0.0D0 01771**3 M=J-1 01772**3 DO 6 K=I,M 01773**3 6 SUM=SUM-Q(I,K)*Q(K,J) 01774**3 Q(I,J)=SUM/Q(J,J) 01775**3 7 CONTINUE 01776**3 DO 9 J=1,N 01777**3 DO 9 I=1,J 01778**3 SUM=0.0D0 01779**3 DO 8 K=J,N 01780**3 8 SUM=SUM+Q(I,K)*Q(J,K) 01781**3 Q(I,J)=SUM 01782**3 IF(I.EQ.J) GO TO 9 01783**3 Q(J,I)=SUM 01784**3 9 CONTINUE 01785**3 RETURN 01786**3 11 INDEX=1 01787**6 RETURN 01788**6 END 01789**3 C 01790**2 C 01791**2 C***********************************************************************01792 C* 01793 C* SUBROUTINE 'G C P L O T' :- ( BY : M.M. NASSAR, JANUARY 1975.) 01794 C* ------------------------ 01795 C* 01796 C* THE PURPOSE OF THIS SUBROUTINE IS TO PLOT01797 C* TWO PERPENDICULAR AXES: X (HORIZONTAL) AND Y (VERTICAL), INTERSECTING01798 C* AT THE CENTRE OF A SQUARE 6 BY 6 INCHES; AND THEN PLOT 'N' NUMBER OF 01799 C* POINTS DEFINED BY THEIR X & Y COORDINATES RELATIVE TO THE CENTRE OF 01800 C* THE SQUARE AND GIVEN IN THE USER'S UNITS. 01801 C* 01802 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*01803 C* 01804 C* NOTE THAT THE TRANSFORMATION OF COORDINATES FROM THE USER UNITS TO TH01805 C* PLOTTING UNITS IS DONE AUTOMATICALLY VIA THE PLOTTING ROUTINES CALLED01806 C* BY THIS SUBROUTINE. FOR MORE DETAILS CONCERNING THE PLOTTING ROUTI01807 C* AVAILABLE AT THE U.N.B. COMPUTING CENTRE, SEE: ''COMPUTER PLOTTING'',01808 C* A USER GUIDE MANUAL BY: U.G. GUJAR, CL # 278, COMPUTING CENTRE, U.N.B01809 C* FREDERICTON, N.B., OCTOBER, 1972. 01810 C* 01811 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*01812 C* 01813 C* USAGE :- CALL GCPLOT(X,Y,N,XMIN,YMIN,XMAX,YMAX) 01814 C* 01815 C*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*01816 C* 01817 C* REQUIRED INPUT DATA :- (VIA THE CALLING ROUTINE) 01818 C* ------------------- 01819 C* 01820 C* X : VECTOR OF X COORDINATES OF THE N-POINTS, IN USER'S UNITS. 01821 C* Y : VECTOR OF Y COORDINATES OF THE N-POINTS, IN USER'S UNITS. 01822 C* N : TOTAL NUMBER OF POINTS TO BE PLOTTED. 01823 C* XMIN : MINIMUM X COORDINATE IN THE USER'S UNITS. 01824 C* YMIN : MINIMUM Y COORDINATE IN THE USER'S UNITS. 01825 C* XMAX : MAXIMUM X COORDINATE IN THE USER'S UNITS. 01826 C* YMAX : MAXIMUM Y COORDINATE IN THE USER'S UNITS. 01827 C* 01828 C* THE OUTPUT WILL BE THE PLOTTED SQUARE, THE TWO PERPENDICULAR AXES, AN01829 C* THE N-POINTS. FINALLY, THIS SUBROUTINE PRINTS-OUT THE TOTAL NUMBER '01830 C* OF PLOTTED POINTS NO LIMITATIONS ON THE NUMBER OF POINTS. 01831 C* 01832 C***********************************************************************01833 C 01834 SUBROUTINE GCPLOT(X,Y,N,XMIN,YMIN,XMAX,YMAX) 01835 IMPLICIT REAL*8(A-H,O-Z) 01836 DIMENSION X(N),Y(N) 01837 XLNGTH=6. 01838 YLNGTH=6. 01839 CALL DEVICE(1627) 01840 CALL AREA(XLNGTH,YLNGTH) 01841 CALL SETPLT(XMIN,YMIN,XMAX,YMAX) 01842 CALL PRNTCH('+') 01843 RECXL=2.0D 0*XMIN 01844 CALL RECT(XMIN,YMIN,RECXL) 01845 CALL AXES(XMIN,YMIN,XMAX,YMAX,0.,0.,-1.,-1.) 01846 CALL PRNTCH('.') 01847 DO 30 K=1,N 01848 CALL NOWPLT(0,X(K),Y(K)) 01849 30 CONTINUE 01850 CALL ENDPLT 01851 WRITE(6,19) N 01852 19 FORMAT(/,5X,'TOTAL NUMBER OF PLOTTED POINTS SURROUNDING THE CENTRA01853 1L-ORIGIN IS =',I6,2X,'POINTS') 01854 RETURN 01855 END 01856 SUBROUTINE SORT(VCLS,NOR,NRES) 01857**4 C***********************************************************************01858**4 C* 01859**4 C* SORT REARRANGES THE ELEMENTS OF VCLS IN ORDER OF INCREASING MAGNITUD01860**4 C* FOR USE IN SEPARATING THESE ELEMENTS INTO CLASS INTERVALS FOR THE CH01861**4 C* SQUARE GOODNESS OF FIT TEST AND HISTOGRAM PLOT. 01862**4 C* 01863**4 C* 01864**4 C* INPUT: 01865**4 C* NRES- NUMBER OF RESIDUALS IN VCLS 01866**4 C* 01867**4 C* 01868**4 C* WRITTEN BY: 01869**4 C* LAURIE PACH, JULY, 1978 01870**4 C* 01871**4 C***********************************************************************01872**4 IMPLICIT REAL*8(A-H,O-Z) 01873**4 DIMENSION VCLS(NOR) 01874**4 NOUT=NRES*(NRES+1)/2 01875**4 NINS=NRES-1 01876**4 DO 12 J=1,NOUT 01877**4 IFLG=0 01878**4 DO 11 I=1,NINS 01879**4 IF(VCLS(I).GT.VCLS(I+1))GOTO10 01880**4 GOTO11 01881**4 10 IFLG=1 01882**4 TEMP=VCLS(I) 01883**4 VCLS(I)=VCLS(I+1) 01884**4 VCLS(I+1)=TEMP 01885**4 11 CONTINUE 01886**4 IF(IFLG.EQ.0)GOTO13 01887**4 12 CONTINUE 01888**4 13 RETURN 01889**4 END 01890**4 SUBROUTINE ILPRNT(NVAL,MAX,PLOTV,RVEC,WINT,N,KK) 01891**4 C***********************************************************************01892**4 C* 01893**4 C* LPRNT CONTROLS THE LINE SPACING FOR THE NORMAL-HISTOGRAM PLOT PRINTI01894**4 C* 01895**4 C* WRITTEN BY: 01896**4 C* LAURIE PACH, JULY, 1978 01897**4 C* 01898**4 C***********************************************************************01899**4 INTEGER WINT,PLOTV,SV,RVEC 01900**4 DIMENSION PLOTV(110),SV(22),RVEC(WINT) 01901**4 DATA SV/' ',' ',' ',' ','R','E','L','A','T','I','V','E',' ','F', 01902**4 @ 'R','E','Q','U','E','N','C','Y'/ 01903**4 B=.3 01904**4 C=.2 01905**4 D=.1 01906**4 DO 2 I=1,100 01907**4 DO 3 N=1,WINT 01908**4 IF(RVEC(N).EQ.MAX)RETURN 01909**4 3 CONTINUE 01910**4 N=N-1 01911**4 IF(NVAL.EQ.MAX)RETURN 01912**4 PRINT101,(PLOTV(L),L=1,110) 01913**4 IF(MAX.EQ.25)PRINT113,B 01914**4 IF(MAX.EQ.17)PRINT113,C 01915**4 IF(MAX.EQ.9)PRINT113,D 01916**4 MAX=MAX-1 01917**4 IF(MAX.GT.32.OR.KK.EQ.23)GOTO2 01918**4 PRINT114,SV(KK) 01919**4 KK=KK+1 01920**4 2 CONTINUE 01921**4 101 FORMAT(' ',6X,110A1) 01922**4 113 FORMAT('+',3X,F3.1,'-') 01923**4 114 FORMAT('+',1X,A1) 01924**4 RETURN 01925**4 END 01926**4 SUBROUTINE PBLANK(PLOTV) 01927**4 C***********************************************************************01928**4 C* 01929**4 C* PBLANK CLEARS (SETS ELEMENTS TO BLANKS) VECTOR PLOTV, WHICH IS 01930**4 C* USED IN SUBROUTINE PLOT. 01931**4 C* 01932**4 C* 01933**4 C* WRITTEN BY: 01934**4 C* LAURIE PACH, JULY, 1978 01935**4 C* 01936**4 C***********************************************************************01937**4 INTEGER PLOTV 01938**4 DIMENSION PLOTV(110) 01939**4 DATA BLNK/' '/ 01940**4 DO 1 I=1,110 01941**4 PLOTV(I)=BLNK 01942**4 1 CONTINUE 01943**4 RETURN 01944**4 END 01945**4 SUBROUTINE PLOT(WINT,HVEC,NB) 01946**5 C***********************************************************************01947**4 C* 01948**4 C* PLOT PLOTS THE STANDARD NORMAL CURVE OVERLAYED WITH THE HISTOGRAM OF01949**4 C* STANDARD RESIDUALS. 01950**4 C* 01951**4 C* 01952**4 C* INPUT: 01953**4 C* WINT- NUMBER OF HISTOGRAM INTERVALS 01954**4 C* HVEC- VECTOR CONTAINING THE NUMBER OF RESIDUALS IN EACH HISTOG01955**4 C* INTERVAL 01956**4 C* 01957**4 C* 01958**4 C* WRITTEN BY: 01959**4 C* LAURIE PACH, JULY, 1978 01960**4 C* 01961**4 C***********************************************************************01962**4 INTEGER WINT,NVEC,RVEC,HVEC,STAR,HLINE,VLINE,PLOTV,PLOTL,PLOTH 01963**4 DIMENSION PLOTV(110),NVEC(53),PLOTL(110),PLOTH(110),RVEC(20), 01964**4 @ HVEC(20),NB(11) 01965**5 DATA STAR,HLINE,VLINE/'.','-','|'/ 01966**4 DATA NVEC/22*0,4*1,2,2,3,4,4,5,6,8,9,10,12,14,16,17,19,21,23,25, 01967**4 @27,28,29,31,31,32,32,32,31/ 01968**4 MAX=50 01969**4 A=.4 01970**4 K=2 01971**4 I=52 01972**4 NUMR=0 01973**4 KK=1 01974**4 NVAL=32 01975**4 NFLG=0 01976**4 IF(WINT.EQ.20)INT=5 01977**4 IF(WINT.EQ.10.OR.WINT.EQ.9)INT=10 01978**4 IF(WINT.EQ.2)INT=50 01979**4 IF(WINT.EQ.3)INT=30 01980**4 IF(WINT.EQ.5)INT=20 01981**4 IF(WINT.EQ.4)INT=25 01982**4 CALL PBLANK(PLOTV) 01983**4 CALL PBLANK(PLOTH) 01984**4 CALL PBLANK(PLOTL) 01985**4 DO 2 JJ=1,WINT 01986**4 NUMR=NUMR+HVEC(JJ) 01987**4 2 CONTINUE 01988**4 DO 33 JJ=1,WINT 01989**4 RVEC(JJ)=(80*HVEC(JJ)/NUMR)+.5 01990**4 33 CONTINUE 01991**4 PRINT103 01992**4 DO 29 JJ=1,50 01993**4 DO 28 N=1,WINT 01994**4 IF(RVEC(N).GE.MAX)GOTO19 01995**4 28 CONTINUE 01996**4 PRINT104 01997**4 IF(N.EQ.(WINT+1))N=WINT 01998**4 IF(MAX.EQ.32)GOTO21 01999**4 MAX=MAX-1 02000**4 29 CONTINUE 02001**4 19 RVEC(N)=0 02002**4 MM=(INT*(N-1))+1 02003**4 IF(WINT.EQ.9.OR.WINT.EQ.3)MM=(INT*(N-1))+6 02004**4 PLOTL(MM)=VLINE 02005**4 III=INT-1 02006**4 DO 12 JJ=1,III 02007**4 PLOTH(JJ+MM)=HLINE 02008**4 12 CONTINUE 02009**4 PLOTL(MM+INT)=VLINE 02010**4 DO 32 N=1,WINT 02011**4 IF(RVEC(N).EQ.MAX.OR.RVEC(N).GT.50)GOTO19 02012**4 32 CONTINUE 02013**4 IF(NFLG.EQ.1)GOTO25 02014**4 PRINT101,(PLOTH(L),L=1,110) 02015**4 CALL PBLANK(PLOTH) 02016**4 CALL ILPRNT(NVAL,MAX,PLOTL,RVEC,WINT,N,KK) 02017**4 IF(RVEC(N).EQ.MAX.AND.MAX.NE.32)GOTO19 02018**4 21 DO 31 L=2,100 02019**4 PLOTV(I)=STAR 02020**4 IF(NVEC(I).NE.NVEC(I-1))GOTO4 02021**4 K=K+1 02022**4 I=I-1 02023**4 31 CONTINUE 02024**4 4 I=I-1 02025**4 NFLG=1 02026**4 IF(RVEC(N).EQ.MAX)GOTO19 02027**4 25 DO 26 L=1,110 02028**4 IF(PLOTH(L).EQ.HLINE)GOTO36 02029**4 26 CONTINUE 02030**4 GOTO38 02031**4 36 DO 39 L=1,110 02032**4 IF(PLOTH(L).EQ.HLINE)PLOTV(L)=PLOTH(L) 02033**4 39 CONTINUE 02034**4 38 CALL PBLANK(PLOTH) 02035**4 11 PRINT102,(PLOTV(L),L=1,110) 02036**4 IF(MAX.EQ.32)PRINT113,A 02037**4 IF(I.EQ.1)GOTO20 02038**4 NVAL=NVEC(I) 02039**4 IF(NVAL.EQ.0)GOTO20 02040**4 CALL PBLANK(PLOTV) 02041**4 DO 37 L=1,110 02042**4 IF(PLOTL(L).EQ.VLINE)PLOTV(L)=PLOTL(L) 02043**4 37 CONTINUE 02044**4 CALL ILPRNT(NVAL,MAX,PLOTL,RVEC,WINT,N,KK) 02045**4 CALL PBLANK(PLOTV) 02046**4 DO 3 L=2,100 02047**4 PLOTV(I)=STAR 02048**4 PLOTV(I+K)=STAR 02049**4 K=K+2 02050**4 IF(I.EQ.1)GOTO11 02051**4 IF(NVEC(I).NE.NVEC(I-1))GOTO4 02052**4 I=I-1 02053**4 3 CONTINUE 02054**4 GOTO4 02055**4 20 PRINT111 02056**4 PRINT112,(NB(IL),IL=1,11) 02057**5 IF(WINT.EQ.20)PRINT107,(HVEC(IL),IL=1,20) 02058**5 IF(WINT.EQ.10)PRINT106,(HVEC(IL),IL=1,10) 02059**5 IF(WINT.EQ.4)PRINT108,(HVEC(IL),IL=1,4) 02060**5 IF(WINT.EQ.2)PRINT109,(HVEC(IL),IL=1,2) 02061**5 IF(WINT.EQ.9)PRINT116,(HVEC(IL),IL=1,9) 02062**5 IF(WINT.EQ.5)PRINT114,(HVEC(IL),IL=1,5) 02063**5 IF(WINT.EQ.3)PRINT115,(HVEC(IL),IL=1,3) 02064**5 101 FORMAT(' ',6X,110A1) 02065**4 102 FORMAT('+',6X,110A1) 02066**4 103 FORMAT('1') 02067**4 104 FORMAT(' ') 02068**4 106 FORMAT(' ',8X,10(I4,6X),/) 02069**4 107 FORMAT(' ',7X,20(I4,1X),/) 02070**4 108 FORMAT(' ',15X,I4,3(21X,I4),/) 02071**4 109 FORMAT(' ',34X,I4,36X,I4,/) 02072**4 111 FORMAT(' ',6X,20('|----'),'|') 02073**4 112 FORMAT(' ',5X,10(I3,7X),I3/) 02074**5 113 FORMAT('+',3X,F3.1,'-') 02075**4 114 FORMAT(' ',13X,I4,4(16X,I4),/) 02076**4 115 FORMAT(' ',23X,3(I4,26X),/) 02077**4 116 FORMAT(' ',14X,9(I4,6X),/) 02078**4 RETURN 02079**4 END 02080**4 C ******************************************************************* 02081**9 C * * 02082**9 C * SUBROUTINE 'FAMNMX' IS DESIGNED TO DETERMINE THE * 02083**9 C * MAXIMUM AND MINIMUM VALUES OF A DATA SERIES. * 02084**9 C * * 02085**9 C * * 02086**9 C * INPUT: * 02087**9 C * * 02088**9 C * AIN - ARRAY OF DATA * 02089*10 C * NAIN - NUMBER OF DATA IN ARRAY * 02090*10 C * * 02091**9 C * OUTPUT: * 02092**9 C * * 02093**9 C * BOUT,AOUT - MAXIMUM AND MINIMUN OF DATA SERIES * 02094*10 C * * 02095**9 C ******************************************************************* 02096**9 C 02097**9 SUBROUTINE FAMNMX(AIN,NAIN,AOUT,BOUT) 02098**5 IMPLICIT REAL*8(A-H,O-Z) 02099**5 DIMENSION AIN(NAIN) 02100**5 AOUT=1000. 02101**5 BOUT=-1000. 02102**5 DO 20 I=1,NAIN 02103**5 IF(AIN(I).LT.AOUT) AOUT=AIN(I) 02104**5 IF(AIN(I).GT.BOUT) BOUT=AIN(I) 02105**5 20 CONTINUE 02106**5 RETURN 02107**5 END 02108**5 C ******************************************************************* 02109**9 C * * 02110**9 C * * 02111**9 C * SUBROUTINE 'GRAV' IS DESIGNED TO PREDICT GRAVITY ANOMALY * 02112**9 C * BY THE METHOD OF LEAST-SQUARES COLLOCATION FOR DETAILS ABOUT * 02113**9 C * THE MATHEMATICAL MODEL SEE REFERENCES IN MAIN PROGRAM 'COLO- * 02114**9 C * CATE'. * 02115**9 C * * 02116**9 C * * 02117**9 C * INPUT DATA: * 02118**9 C * ----------- * 02119**9 C * * 02120**9 C * X, VDG - OBSERVATIONS AND THEIR VARIANCES * 02121*10 C * NDG - NUMBER OF OBSERVATIONS * 02122*10 C * DLT, DLN - POSITION COORDINATES OF EACH OBSERVATION * 02123*10 C * XLAT,XLOG - POSITION COORDINATES OF PREDICTION POINT * 02124*10 C * NBRR - NUMBER OF ELEMENTS IN UPPER TRIANGLE AND DIAGONAL * 02125*10 C * OF COVARIANCE MATRIX (CX) OF OBSERVATIONS * 02126*10 C * IPRINT - PARAMETER SPECIFYING WHETHER (OR NOT) TO PRINT * 02127*10 C * INTERMEDIATE RESULTS * 02128*10 C * * 02129**9 C * * 02130**9 C * OUTPUT: * 02131**9 C * ------- * 02132**9 C * * 02133**9 C * * 02134**9 C * DGP - VALUE OF PREDICTED ANOMALY * 02135*10 C * VAR - VARIANCE OF PREDICTED ANOMALY * 02136*10 C * * 02137*11 C * ELEMENTS OF COVARIANCE MATRICES, VECTORS (IF REQUIRED) : * 02138*10 C * * 02139*10 C * CXX - VECTOR OF ALL ELEMENTS OF UPPER TRIANGLE AND * 02140*10 C * DIAGONAL OF COVARIANCE MATRIX (CX) OF MEASURE- * 02141*10 C * MENTS TAKING THE ROWS SUCCESSIVELY; SEE BELOW: * 02142*10 C * * 02143*10 C * 1 2 3 * 02144*10 C * 4 5 ---> ( 1 2 3 4 5 6 ) * 02145*11 C * 6 * 02146*11 C * CX,CX1,CX2 - COVARIANCE MATRICES OF MEASUREMENTS * 02147*10 C * CX1,CX2 WERE USED TO COPY CX. THEY COULD * 02148*10 C * BE DISCARDED FROM THE SUBROUTINE AND MAIN * 02149*10 C * PROGRAM * 02150*10 C * CSX,CSXA - SIGNAL COVARIANCE MATRIX AND ITS TRANSPOSE * 02151*10 C * ( CSX1 ) CSX1 WAS USED TO COPY CSX. * 02152*10 C * CSXB - SOLUTION MATRIX ( CSXB = CSX * CX(INVERSE) ) * 02153*10 C * GG - AUTO-COVARIANCE OF MEASUREMENTS. * 02154*10 C * * 02155**9 C ******************************************************************* 02156**9 C 02157**9 C 02158**5 C SUBROUTINE GRAV 02159**5 C =============== 02160**5 02161**5 SUBROUTINE GRAV(NDG,VDG,DLT,DLN,XRES,CXX,CSX,CSXA,CSXB, 02162**5 *XLAT,XLOG,DGP,VAR,NBRR,CX,CX2,GG,CX1,CSX1,IPRINT,INDEX) 02163**6 C...PREPARED BY G. LACHAPELLE, TU GRAZ, OCTOBER 75. 02164**5 C...MODIFIED IN FEBRUARY 75. 02165**5 IMPLICIT REAL*8(A-H,O-Z) 02166**5 COMMON /CMCOV/CI(12),CR(51),SIGMA0(300),SIGMA(300), 02167**5 *KI(25),N1,LOCAL 02168**5 DIMENSION DLT(50),DLN(50),XRES(50),CXX(1500),CSX(50),VDG(50), 02169**5 *CSXA(1,50),CSX1(50,1),CSXB(1,50),CX(50,50),CX1(50,50),CX2(50,50), 02170**5 * NIK(50),NEK(50) 02171**5 C 02172**5 10 FORMAT(I4,21F6.1) 02173**6 15 FORMAT(//,4X,4F12.5) 02174**5 16 FORMAT(I4,21F6.1) 02175**6 20 FORMAT(//) 02176**5 30 FORMAT(//,4X,'CXX-MATRIX'/,4X,10('=')/) 02177**5 40 FORMAT(//,4X,'CXX-INVERSE'/,4X,11('=')/) 02178**5 50 FORMAT(//,4X,'DETERMINANT=',F7.4,'D',I4) 02179**5 60 FORMAT(//,4X,'CSX-MATRIX'/,4X,13('-')/) 02180**5 70 FORMAT(//,4X,'CSXB=CSX*CXX(INVERSE)'/,4X,22('-')/) 02181**5 80 FORMAT(//,20X,'NEGATIVE OR ZERO ELEMENT ON DIAGONAL OF CXX-MATRIX:02182**6 * CHECK CAUSE:'//) 02183**6 CR(4)=DSIN(DLT(1)) 02184**5 CR(6)=DCOS(DLT(1)) 02185**5 CR(5)=CR(4) 02186**5 CR(7)=CR(6) 02187**5 CR(8)=0.0D0 02188**5 CR(9)=1.0D0 02189**5 CR(1)=CR(4)*CR(5)+CR(6)*CR(7)*CR(9) 02190**5 GG=CVFUNC(CR(1)) 02191**6 NDG1=NDG-1 02192**5 IV=0 02193**5 DO 120 II=1,NDG1 02194**5 IV=IV+1 02195**5 CXX(IV)=GG+VDG(II) 02196**5 K=II+1 02197**5 DO 120 J=K,NDG 02198**5 IV=IV+1 02199**5 CR(4)=DSIN(DLT(II)) 02200**5 CR(6)=DCOS(DLT(II)) 02201**5 CR(5)=DSIN(DLT(J)) 02202**5 CR(7)=DCOS(DLT(J)) 02203**5 CR(8)=DSIN(DLN(II)-DLN(J)) 02204**5 CR(9)=DCOS(DLN(II)-DLN(J)) 02205**5 IF((DABS(DLT(II)-DLT(J)).LT.1.0D-5).AND.(DABS(DLN(II)-DLN(J)).LT. 02206**5 *1.0D-5)) GO TO 115 02207**5 CR(1)=CR(4)*CR(5)+CR(6)*CR(7)*CR(9) 02208**5 GO TO 120 02209**5 115 CR(1)=1.0D0 02210**5 120 CXX(IV)=CVFUNC(CR(1)) 02211**6 IV=IV+1 02212**5 CXX(IV)=GG+VDG(NDG) 02213**5 C...VECTOR CSX %CSXA # CSX< 02214**5 CR(4)=DSIN(XLAT) 02215**5 CR(6)=DCOS(XLAT) 02216**5 DO 150 II=1,NDG 02217**5 CR(5)=DSIN(DLT(II)) 02218**5 CR(7)=DCOS(DLT(II)) 02219**5 CR(8)=DSIN(DLN(II)-XLOG) 02220**5 CR(9)=DCOS(DLN(II)-XLOG) 02221**5 IF((DABS(XLAT-DLT(II)).LT.1.0D-5).AND.(DABS(XLOG-DLN(II)) 02222**5 *.LT.1.0D-5)) GO TO 140 02223**5 CR(1)=CR(4)*CR(5)+CR(6)*CR(7)*CR(9) 02224**5 GO TO 145 02225**5 140 CR(1)=1.0D0 02226**5 145 CSX(II)=CVFUNC(CR(1)) 02227**6 CSX1(II,1)=CSX(II) 02228**5 150 CSXA(1,II)=CSX(II) 02229**5 N=1 02230**5 DO 162 J=1,NDG 02231**5 CX(1,J)=CXX(N) 02232**5 N=N+1 02233**5 162 CONTINUE 02234**5 DO 157 I=2,NDG 02235**5 DO 157 K=I,NDG 02236**5 CX(I,K)=CXX(N) 02237**5 N=N+1 02238**5 157 CONTINUE 02239**5 DO 158 I=2,NDG 02240**5 I1=I-1 02241**5 DO 158 K=1,I1 02242**5 158 CX(I,K)=CX(K,I) 02243**5 IF(IPRINT.LT.2) GO TO 161 02244**5 WRITE(6,30) 02245**5 DO 159 I=1,NDG 02246**5 WRITE(6,10) I,(CX(I,J),J=1,NDG) 02247**5 159 CONTINUE 02248**5 WRITE(6,20) 02249**5 161 IF(IPRINT.EQ.1) GO TO 141 02250**7 GO TO 160 02251**7 141 CONTINUE 02252**7 WRITE(6,144) 02253**7 DO 142 I=1,NDG 02254**7 WRITE(6,143) I,CX(I,I) 02255**7 142 CONTINUE 02256**7 143 FORMAT(20X,I5,F12.5) 02257**7 144 FORMAT(//,20X,'DIAGONAL ELEMENTS OF CXX-MATRIX'/,20X,31('=')//) 02258**7 160 CONTINUE 02259**7 C... SCALING OF CXX(CX) AGAINST UNDERFLOW OR OVERFLOW DURING INVERSION 02260**5 C... DIVIDE EACH ELEMENT BY S=GG=588.13 GGAL**2. 02261**5 SS1=GG 02262**5 DO 163 I=1,NDG 02263**5 DO 163 J=1,NDG 02264**5 CX(I,J)=CX(I,J)/SS1 02265**5 163 CONTINUE 02266**5 C...SOLUTION 02267**5 CALL SPIN(CX,NDG,50,DET,IDEXP,INDEX) 02268**6 IF(INDEX.EQ.1) GO TO 175 02269**6 C 02270**5 IF(IPRINT.EQ.2) WRITE(6,50) DET,IDEXP 02271**5 C... REFILL THE MATRIX CX AFTER INVERSION: 02272**5 DO 164 I=2,NDG 02273**5 I1=I-1 02274**5 DO 164 J=1,I1 02275**5 CX(I,J)=CX(J,I) 02276**5 164 CONTINUE 02277**5 C... RESCALE THE INVERTED CX MATRIX: 02278**5 DO 165 I=1,NDG 02279**5 DO 165 J=1,NDG 02280**5 CX(I,J)=CX(I,J)/SS1 02281**5 165 CONTINUE 02282**5 C... COMPUTE CSXB=CSXA*CXX(INV): 02283**5 DO 170 I=1,NDG 02284**5 CSXB(1,I)=0.0D0 02285**5 170 CONTINUE 02286**5 DO 171 I=1,NDG 02287**5 DO 171 J=1,NDG 02288**5 CSXB(1,I)=CSXB(1,I)+CSXA(1,J)*CX(J,I) 02289**5 171 CONTINUE 02290**5 168 FORMAT(10X,I5,F9.2) 02291**5 DGP=0.0 02292**5 SOM=0.0 02293**5 DO 200 II=1,NDG 02294**5 DGP=DGP+CSXB(1,II)*XRES(II) 02295**5 SOM=SOM+CSXB(1,II)*CSX1(II,1) 02296**5 200 CONTINUE 02297**5 VAR=GG-SOM 02298**5 IF(IPRINT.LT.2) GO TO 167 02299**5 WRITE(6,40) 02300**5 DO 166 I=1,NDG 02301**5 166 WRITE(6,16) I,(CX(I,J),J=1,NDG) 02302**5 WRITE(6,60) 02303**5 I=1 02304**5 WRITE(6,10) I,(CSXA(1,J),J=1,NDG) 02305**5 WRITE(6,70) 02306**5 WRITE(6,10) I,(CSXB(1,J),J=1,NDG) 02307**5 WRITE(6,15) DGP,GG,SOM,VAR 02308**5 WRITE(6,20) 02309**5 167 CONTINUE 02310**5 GO TO 180 02311**6 175 WRITE(6,80) 02312**6 180 RETURN 02313**6 END 02314**5 C ******************************************************************* 02315**9 C * * 02316**9 C * SUBROUTINE 'PRECOV' IS DESIGNED TO COMPUTE AVERAGE COVAR- * 02317**9 C * IANCES CV AND CORRESPONDING DISTANCES S BETWEEN PAIRS OF * 02318**9 C * (RESIDUAL ANOMALY DATA. CV AND S ARE REQUIRED FOR THE DETER- * 02319**9 C * MINATION OF COEFFICIENTS CI OF THE COVARIANCE FUNCTION EXPRE- * 02320**9 C * SSED AS: * 02321**9 C * * 02322**9 C * CV(S) = C0+C1*S+C2*S**2+C3*S**3 * 02323**9 C * * 02324**9 C * WHICH IS AN ALGEBRAIC POLYNOMIAL IN S. C0, C1, C2, AND C3 ARE * 02325**9 C * THE COEFFICIENTS CI . * 02326**9 C * * 02327**9 C * * 02328**9 C * INPUT DATA: * 02329**9 C * ----------- * 02330**9 C * * 02331**9 C * DELG,VRDG - OBSERVED GRAVITY ANOMALIES AND THEIR VARIANCES * 02332*10 C * DLAT,DLON - POSITIONAL COORDINATES OF EACH OBSERVATION * 02333*10 C * IPR - NUMBER OF OBSERVATIONS * 02334*10 C * * 02335**9 C * * 02336**9 C * OUTPUT * 02337**9 C * ------ * 02338**9 C * * 02339**9 C * CVR - AVERAGE COVARIANCES FOR DISTACCES 0 - 60 MINUTES * 02340*10 C * AT 5 MINUTE INTERVAL * 02341*10 C * ICVR - CORRESPONDING CLASS NUMBERS: 1 --> 0 - 5 * 02342*10 C * 2 --> 6 - 10 ETC. * 02343*10 C * CLASS - CLASS WIDTH * CLASS NUMBER * 02344*10 C * SIGCVR - STANDARD DEVIATION OF COVARIANCES * 02345*10 C * AVGPSI - APPROXIMATE VALUE FOR CLASS * 02346*10 C * * 02347**9 C ******************************************************************* 02348**9 C 02349**9 SUBROUTINE PRECOV(IPR,DLAT,DLON,DELG,VRDG,CVR,SIGCVR,CLASS,AVGPSI,02350**6 *ICVR) 02351**6 C 02352**9 IMPLICIT REAL*8(A-H,O-Z) 02353**6 DIMENSION DLAT(500),DLON(500),DELG(500),VRDG(500),CRR(6),CVR(12), 02354**6 *AVGPSI(12),SIGCVR(12) 02355**6 INTEGER ICVR(12) 02356**6 10 FORMAT(2I2) 02357**6 20 FORMAT(20X,'*** FILE NUMBER',I3,'***'//) 02358**6 25 FORMAT(20X,'DETERMINE LEAST SQUARE S LINE PQRQMETERS'/,32X,10('*')02359**6 *//) 02360**6 30 FORMAT(2F10.5,2F11.5,F8.2,3F10.5) 02361**6 40 FORMAT(24X,'__',10X,'__'/,24X,'DG',10X,'DGSQ',14X,'NRDG'/,17X, 02362**6 *3(F10.2,5X),I10//) 02363**6 50 FORMAT(20X,2F10.5,2(D15.7,I10)) 02364**6 51 FORMAT(20X,F5.2,7X,F5.2,2X,D11.4,3X,D11.4,5X,I10) 02365**6 52 FORMAT(20X,'THEORETICL',2X,'ACTUAL'/,20X,'CLASS',7X,'CLASS',3X, 02366**6 *'COVARIANCE',3X,'STD. DEV',5X,'NR/ANOMALIES'/) 02367**6 60 FORMAT(90X,D15.7,I10) 02368**6 70 FORMAT(//,20X,'B0,B1=',2F12.4,'STANDARD DEV. B1=',F10.5//,20X, 02369**6 *'CORRELATION COEFFICIENT=',F12.7//,20X,'STANDARD DEV. REGR=', 02370**6 *F12.7//) 02371**6 C... 02372**6 CNVRT=11.04656D 00/60.D 00 02373**6 C0=0.D0 02374**6 VARC0=0.D0 02375**6 DEL3=0.D0 02376**6 PI=DARCOS(-1.D0) 02377**6 RHO=PI/180.D0 02378**6 II=IPR 02379**6 DO 100 I=1,II 02380**6 DEL2=DELG(I)**2 02381**6 DEL3=DEL3+DELG(I) 02382**6 C0=C0+DEL2 02383**6 VARC0=VARC0+DEL2*VRDG(I) 02384**6 100 CONTINUE 02385**6 C0CT=DFLOAT(II) 02386**6 C0=C0/C0CT 02387**6 DEL3=DEL3/C0CT 02388**6 VARC0=VARC0/C0CT**2*4.D0 02389**6 WRITE(6,40) DEL3,C0,VARC0,II 02390**6 C... 02391**6 DO 260 I=1,12 02392**6 CVR(I)=0.D0 02393**6 SIGCVR(I)=0.D0 02394**6 AVGPSI(I)=0.D0 02395**6 ICVR(I)=0 02396**6 260 CONTINUE 02397**6 G0=0.D0 02398**6 IG0=0 02399**6 I1=II-1 02400**6 201 DO 250 I=1,I1 02401**6 I2=I+1 02402**6 DO 245 J=I2,II 02403**6 CRR(2)=DSIN(DLAT(I)) 02404**6 CRR(3)=DSIN(DLAT(J)) 02405**6 CRR(4)=DCOS(DLAT(I)) 02406**6 CRR(5)=DCOS(DLAT(J)) 02407**6 CRR(6)=DCOS(DLON(J)-DLON(I)) 02408**6 CRR(1)=CRR(2)*CRR(3)+CRR(4)*CRR(5)*CRR(6) 02409**6 PSI=DARCOS(CRR(1)) 02410**6 PSI=(PSI/RHO)*60.D0 02411**6 ADIS=PSI*CNVRT 02412**6 IF((ADIS.LE.1.).AND.(ADIS.GT.0.)) GO TO 202 02413**6 GO TO 204 02414**6 202 CONTINUE 02415**6 IG0=IG0+1 02416**6 DG0=(DELG(I)-DELG(J))/ADIS 02417**6 DG0=DABS(DG0) 02418**6 G0=G0+DG0 02419**6 204 CONTINUE 02420**6 CALL CLASFY(PSI,IP) 02421**6 IF(IP.EQ.0) GO TO 245 02422**6 CVR(IP)=CVR(IP)+DELG(I)*DELG(J) 02423**6 SIGCVR(IP)=SIGCVR(IP)+DELG(I)**2*VRDG(J)+DELG(J)**2*VRDG(I) 02424**6 AVGPSI(IP)=AVGPSI(IP)+PSI 02425**6 ICVR(IP)=ICVR(IP)+1 02426**6 245 CONTINUE 02427**6 250 CONTINUE 02428**6 WRITE(6,60) G0,IG0 02429**6 IF(IG0.EQ.0) GO TO 252 02430**6 G0CT=DFLOAT(IG0) 02431**6 G0=G0/G0CT 02432**6 WRITE(6,60) G0,IG0 02433**6 252 CONTINUE 02434**6 WRITE(6,52) 02435**6 DO 261 I=1,12 02436**6 IF(ICVR(I).EQ.0) GO TO 261 02437**6 COUNT=DFLOAT(ICVR(I)) 02438**6 CVR(I)=CVR(I)/COUNT 02439**6 AVGPSI(I)=AVGPSI(I)/COUNT 02440**6 SIGCVR(I)=SIGCVR(I)/COUNT**2 02441**6 SCVR=DSQRT(SIGCVR(I)) 02442**6 CLASS=I*5.D0 02443**6 WRITE(6,51) CLASS,AVGPSI(I),CVR(I),SCVR,ICVR(I) 02444**6 261 CONTINUE 02445**6 RETURN 02446**6 END 02447**6 C ******************************************************************* 02448**9 C * * 02449**9 C * SUBROUTINE 'CLASFY' IS DESIGNED TO PLACE AN OBSERVATION * 02450**9 C * FROM A SERIES INTO A CLASS , GIVEN CLASS WIDTH= 5 MINUTES * 02451**9 C * * 02452**9 C * * 02453**9 C * INPUT DATA: * 02454**9 C * ----------- * 02455**9 C * * 02456**9 C * PP - SPHERICAL DISTANCE ASSOCIATED WITH THE OBSERVATION* 02457*10 C * * 02458**9 C * OUTPUT: * 02459**9 C * ------- * 02460**9 C * * 02461**9 C * ICL - CLASS NUMBER: 1 --> 0 - 5 * 02462*10 C * 2 --> 6 - 10 * 02463**9 C * 10 -->46 -50 ETC. * 02464**9 C * * 02465**9 C ******************************************************************* 02466**9 C 02467**9 SUBROUTINE CLASFY(PP,ICL) 02468**6 IMPLICIT REAL*8(A-H,O-Z) 02469**6 ICL=0 02470**6 DCL=5.D0 02471**6 CL1=0.D0 02472**6 CL2=5.D0 02473**6 DO 50 JI=1,12 02474**6 IF((PP.GT.CL1).AND.(PP.LE.CL2)) GO TO 55 02475**6 CL1=CL2 02476**6 CL2=CL2+DCL 02477**6 GO TO 50 02478**6 55 CONTINUE 02479**6 ICL=JI 02480**6 GO TO 60 02481**6 50 CONTINUE 02482**6 60 RETURN 02483**6 END 02484**6 C ******************************************************************* 02485**9 C * FUNCTION-SUBPROGRAM 'CVFUNC' IS DESIGNED TO EVALUATE THE * 02486**9 C * COVARIANCE FUNCTION EXPRESSED AS: * 02487**9 C * * 02488**9 C * CV(S) = C0 + C1*S + C2*S**2 + C3*S**3 * 02489**9 C * * 02490**9 C * AN ALGEBRAIC POLYNOMIAL IN S * 02491**9 C * * 02492**9 C * THE COEFFICIENTS C0,C1,C2,C3 ARE SPECIFIED THROUGH A DATA * 02493**9 C * STATEMENT. IF A NEW SET OF COEFFICIENTS ARE TO BE USED , THE * 02494**9 C * DATA STATEMENT MUST BE ADJUSTED ACCORDINGLY. * 02495**9 C * * 02496**9 C * * 02497**9 C * INPUT: * 02498*10 C * * 02499*10 C * PSIR - SPHERICAL DISTANCE FOR WHICH COVARIANCE * 02500*10 C * IS REQUIRED * 02501*10 C * * 02502*10 C * * 02503**9 C * OUTPUT: * 02504*10 C * * 02505*10 C * - COVARIANCE * 02506*10 C * * 02507**9 C ******************************************************************* 02508**9 C 02509**9 C 02510**6 DOUBLE PRECISION FUNCTION CVFUNC(PSIR) 02511**6 C 02512**6 IMPLICIT REAL*8(A-H,O-Z) 02513**6 DIMENSION C(3) 02514**6 DATA PI/3.1415926535D0/,C0,C(1),C(2),C(3)/111.9147D0,-625.6191D0, 02515**6 *1018.4520D0,-497.8900D0/ 02516**6 RHO=PI/180.D0 02517**6 PSI=DARCOS(PSIR) 02518**6 S=(PSI/RHO) 02519**6 CS=C0 02520**6 DO 10 I=1,3 02521**6 CS=CS+C(I)*S**I 02522**6 10 CONTINUE 02523**6 CVFUNC=CS 02524**6 RETURN 02525**6 END 02526**6 C ******************************************************************* 02527**9 C * * 02528**9 C * SUBROUTINE 'CHISQ' IS DESIGNED TO CARRY OUT CHI-SQUARE * 02529**9 C * GOODNESS OF FIT TEST ON A DATA SERIES. * 02530**9 C * * 02531**9 C * * 02532**9 C * INPUT DATA: * 02533**9 C * ----------- * 02534**9 C * * 02535**9 C * DG4 - ARRAY OF STANDARDIZED OBSERVATIONS * 02536*10 C * NV - NUMBER OF OBSERVATIONS * 02537*10 C * NCHI - NUMBER OF CLASSES * 02538*10 C * SEG - DELIMITERS OF CLASSES * 02539*10 C * CHI - CHI-SQUARE VALUE FROM TABLES * 02540*10 C * NCHI - NUMBER OF KNOWN POPULATION PARAMETERS * 02541*10 C * * 02542**9 C * * 02543**9 C * OUTPUT: * 02544**9 C * ------- * 02545**9 C * * 02546**9 C * CHISQD - COMPUTED CHI-SQUARE VALUE * 02547*10 C * NDOF - DEGREE OF FREEDOM * 02548*10 C * ICHI - COUNTER FOR OBS. IN EACH CLASS * 02549*10 C * * 02550**9 C ******************************************************************* 02551**9 C 02552**9 SUBROUTINE CHISQ(DG7,NV,SEG,ICHI,NCHI,CHISQD,CHI) 02553**8 IMPLICIT REAL*8(A-H,O-Z) 02554**7 DIMENSION DG7(NV),SEG(NCHI) 02555**8 INTEGER ICHI(NCHI),NV,ICOT,NSEG,NSEG1 02556**7 NSEG=NCHI-1 02557**7 DO 10 K=1,NCHI 02558**7 10 ICHI(K)=0 02559**7 J=1 02560**7 ICOT=0 02561**7 DO 20 I=1,NV 02562**7 IF(DG7(I).GE.SEG(J)) ICHI(J)=ICHI(J)+1 02563**8 20 CONTINUE 02564**7 ICOT=ICOT+ICHI(J) 02565**7 J=J+1 02566**7 JM1=J-1 02567**7 25 IF(J.GT.NSEG) GO TO 35 02568**7 DO 30 I=1,NV 02569**7 IF((DG7(I).LT.SEG(JM1)).AND.(DG7(I).GE.SEG(J))) ICHI(J)=ICHI(J)+1 02570**8 30 CONTINUE 02571**7 ICOT=ICOT+ICHI(J) 02572**7 IF(ICOT.EQ.NV) GO TO 100 02573**7 JM1=J 02574**7 J=J+1 02575**7 GO TO 25 02576**7 35 CONTINUE 02577**7 ICHI(J)=NV-ICOT 02578**7 100 EXPD=NV/DFLOAT(NCHI) 02579**7 CHISQD=0.D0 02580**7 DO 102 K=1,NCHI 02581**7 IF(ICHI(K).EQ.0) GO TO 101 02582**7 CHISQD=CHISQD+((DFLOAT(ICHI(K))-EXPD)**2)/EXPD 02583**7 GO TO 102 02584**7 101 CHISQD=CHISQD+EXPD 02585**7 102 CONTINUE 02586**7 C 02587**7 NDOF=NCHI-2 02588**7 WRITE(6,128) NDOF,CHISQD,CHI 02589**7 128 FORMAT(///,50X,'CHI-SQUARE GOODNESS OF FIT TEST:'//,50X,'DEGREE OF02590**7 * FREEDOM=',I6/,42X,'COMPUTED CHI-SQUARE VALUE=',F12.5/,39X,'CHI-SQ02591**7 *UARE VALUE FROM TABLES=',F12.5) 02592**7 RETURN 02593**7 END 02594**7 C ******************************************************************* 02595**9 C * * 02596**9 C * SUBROUTINE 'AMNSDR' IS DESIGNED TO COMPUTE THE MEAN * 02597**9 C * AND STANDARD DEVIATION OF A SERIES OF OBSERVATIONS. * 02598**9 C * * 02599**9 C * * 02600**9 C * INPUT: * 02601**9 C * * 02602**9 C * RVAR - ARRAY OF OBSERVATIONS * 02603*10 C * NRVAR - NUMBER OF OBSERVATIONS * 02604*10 C * * 02605**9 C * * 02606**9 C * OUTPUT: * 02607**9 C * * 02608**9 C * * 02609**9 C * AMN - MEAN OF OBSERVATIONS * 02610*10 C * SDR - STANDARD DEVIATION OF THE MEAN * 02611*10 C * * 02612**9 C ******************************************************************* 02613**9 C 02614**9 SUBROUTINE AMNSDR(RVAR,NRVAR,AMN,SDR) 02615**7 C 02616**7 IMPLICIT REAL*8(A-H,O-Z) 02617**7 DIMENSION RVAR(NRVAR) 02618**7 RVSUM=0.D0 02619**7 RVSQSM=0.D0 02620**7 DO 10 I=1,NRVAR 02621**7 RVSUM=RVSUM+RVAR(I) 02622**7 RVSQSM=RVSQSM+RVAR(I)**2 02623**7 10 CONTINUE 02624**7 ANRV=DFLOAT(NRVAR) 02625**7 AMN=RVSUM/ANRV 02626**7 SDRSQ=(RVSQSM*ANRV-RVSUM**2)/(ANRV*(ANRV-1.)) 02627**7 SDR=DSQRT(SDRSQ) 02628**7 RETURN 02629**7 END 02630**7 C ******************************************************************* 02631**9 C * * 02632**9 C * SUBROUTINE 'HSTPLT' IS DESIGNED TO PLOT THE HISTOGRAM * 02633**9 C * OF A SERIES OF OBSERVATIONS. IT ASSEMLES THE OBSERVATIONS * 02634**9 C * INTO THEIR RESPECTIVE CLASSES AND CALLS SUBROUTINE 'PLOT' * 02635**9 C * TO PLOT THE HISTOGRAM. * 02636**9 C * * 02637**9 C * * 02638**9 C * INPUT DATA: DDD - ARRAY OF OBSERVATIONS * 02639*10 C * NDDD- NUMBER OF OBSERVATIONS. * 02640*10 C * * 02641**9 C * OUTPUT: - HISTOGRAM PLOT. * 02642**9 C * * 02643**9 C ******************************************************************* 02644**9 C 02645**9 SUBROUTINE HSTPLT(DDD,NDDD,NCLASS,HVEC,NB) 02646*10 IMPLICIT REAL*8(A-H,O-Z) 02647**7 INTEGER HVEC 02648**7 DIMENSION DDD(NDDD),HVEC(20),NB(11) 02649**7 NCLASS=20 02650**7 ACLASS=DFLOAT(NCLASS) 02651**7 DO 120 J=1,NCLASS 02652**7 120 HVEC(J)=0 02653**7 20 CALL FAMNMX(DDD,NDDD,GMN,GMX) 02654**7 GM1=DABS(GMN)+0.5 02655**7 GN=DABS(GMX)+0.5 02656**7 GR=GM1/GN 02657**7 IF(GR.GE.1) GO TO 559 02658**7 XX=-GN 02659**7 AXX=2.*GN/ACLASS 02660**7 GO TO 560 02661**7 559 CONTINUE 02662**7 XX=-GM1 02663**7 AXX=2.*GM1/ACLASS 02664**7 560 CONTINUE 02665**7 XX1=XX+AXX 02666**7 NB(1)=XX 02667**7 NC=NCLASS/2+1 02668**7 DO 561 I=2,NC 02669**7 IM1=I-1 02670**7 NB(I)=XX+IM1*AXX*2. 02671**7 561 CONTINUE 02672**7 K=1 02673**7 121 CONTINUE 02674**7 DO 122 J=1,NDDD 02675**7 IF((DDD(J).LT.XX).OR.(DDD(J).GT.XX1)) GO TO 122 02676**7 HVEC(K)=HVEC(K)+1 02677**7 122 CONTINUE 02678**7 K=K+1 02679**7 XX=XX1 02680**7 XX1=XX+AXX 02681**7 IF(K.GT.NCLASS) GO TO 123 02682**7 GO TO 121 02683**7 123 CONTINUE 02684**7 CALL PLOT(NCLASS,HVEC,NB) 02685**8 RETURN 02686**7 END 02687**7 C ******************************************************************* 02688**9 C * SUBROUTINE 'TONMN' IS DESIGNED TO PERFORM THE STUDENT'S T * 02689**9 C * TEST ON THE MEAN OF A SERIES OF OBSERVATIONS. * 02690**9 C * * 02691**9 C * * 02692**9 C * INPUT: RMN - MEAN VALUE * 02693*10 C * ST - STANDARD DEVIATION * 02694*10 C * NST - NUMBER OF OBSERVATIONS * 02695*10 C * NDEX - PARAMETER SPECIFYING WHETHER(OR NOT) * 02696*10 C * TO USE FACTOR FOR STANDARD NORMAL CURVE * 02697*10 C * * 02698**9 C * OUTPUT: - MESSAGE INDICATING WHETHER TEST * 02699**9 C * PASSED OR FAILED. * 02700**9 C * * 02701**9 C ******************************************************************* 02702**9 C 02703**9 SUBROUTINE TONMN(RMN,ST,NST,NDEX) 02704**7 IMPLICIT REAL*8(A-H,O-Z) 02705**7 TDFALF=1.98 02706**7 ANRM=1.96D0 02707**7 AST=DFLOAT(NST) 02708**7 AMULT=ST/DSQRT(AST) 02709**7 IF(NDEX.EQ.0) PROD=AMULT*ANRM 02710**7 PROD=AMULT*TDFALF 02711**7 EX1=RMN-PROD 02712**7 EX2=RMN+PROD 02713**7 IX=NST-1 02714**7 WRITE(6,124) IX,EX1,EX2 02715**7 124 FORMAT(///,50X,'STUDENTS-T TEST ON MEAN DIFFERENCE:'//,50X,'DEGREE02716**7 * OF FREEDOM=',I6/,51X,'CONFIDENCE LEVEL= 95%'//,50X,F10.5,' < 0 02717**7 *< ',F10.5) 02718**7 IF((EX1.LT.0.).AND.(EX2.GT.0.)) GO TO 125 02719**7 WRITE(6,126) 02720**7 GO TO 99 02721**7 125 WRITE(6,127) 02722**7 126 FORMAT(//,50X,'TEST FAILS'/,50X,10('*')//) 02723**7 127 FORMAT(//,50X,'==========> TEST PASSES <=========='//) 02724**7 99 RETURN 02725**7 END 02726**7 C ******************************************************************* 02727**9 C * SUBROUTINE 'WTDAVG' IS DESIGNED TO COMPUTE A WEIGHTTD * 02728**9 C * AVERAGE OF A SERIES OF OBSERVATIONS, USING AS WEIGHT, * 02729**9 C * DISTANCE DWT RAISED TO SOME POWER IDWT. * 02730**9 C * * 02731**9 C * * 02732**9 C * INPUT DATA: * 02733**9 C * ----------- * 02734**9 C * * 02735**9 C * DELG - ARRAY OF OBSERVATIONS * 02736*10 C * VDELG - ARRAY OF VARIANCES OF OBSERVATIONS * 02737*10 C * DWT - CORRESPONDING ARRAY OF DISTANCES FOR WEIGHTING * 02738*10 C * NDLG - NUMBER OF OBSERVATIONS * 02739*10 C * IDWT - POWER FOR ASSIGNING WEIGHTS * 02740*10 C * TVAL - TEST VALUE USED TO AVOID DIVISION BY ZERO * 02741*10 C * * 02742**9 C * OUTPUT: * 02743**9 C * ------- * 02744**9 C * * 02745**9 C * * 02746**9 C * DELGP - WEIGHTED MEAN * 02747*10 C * VDELGP - STANDARD DEVIATION * 02748*10 C * * 02749**9 C ******************************************************************* 02750**9 C 02751**9 SUBROUTINE WTDAVG(NDLG,DELG,VDELG,DWT,IDWT,TVAL,DELGP,VDELGP) 02752**7 IMPLICIT REAL*8(A-H,O-Z) 02753**7 DIMENSION DELG(NDLG),VDELG(NDLG),DWT(NDLG) 02754**7 C 02755**7 SMWT=0. 02756**7 WTDSM=0. 02757**7 WDSQSM=0. 02758**7 C 02759**7 DO 100 I=1,NDLG 02760**7 IF(DWT(I).EQ.0) GO TO 100 02761**7 WTINV=DWT(I)**IDWT 02762**7 IF(WTINV.LE.TVAL) GO TO 10 02763**7 AWT=1.D00/WTINV 02764**7 GO TO 12 02765**7 10 AWT=1.0D 16 02766**7 12 CONTINUE 02767**7 AWTSQ=AWT*AWT 02768**7 SMWT=SMWT+AWT 02769**7 WTDSM=WTDSM+DELG(I)*AWT 02770**7 WDSQSM=WDSQSM+AWTSQ*VDELG(I) 02771**7 100 CONTINUE 02772**7 SMWTSQ=SMWT*SMWT 02773**7 DELGP=WTDSM/SMWT 02774**7 VDELGP=WDSQSM/SMWTSQ 02775**7 RETURN 02776**7 END 02777**7 C ******************************************************************* 02778**9 C * * 02779**9 C * SUBROUTINE 'REJECT' IS DESIGNED TO PERFORM THE TEST OF * 02780**9 C * OUTLIERS ON A SERIES OF OBSERVATIONS. * 02781**9 C * INPUT: * 02782**9 C * * 02783**9 C * OBS - ARRAY OF OBSERVATIONS * 02784*10 C * RMS - STANDARD DEVIATION * 02785*10 C * IOBS - NUMBER OF OBSERVATIONS * 02786*10 C * NLOC1 - MAXIMUM NUMBER OF OBSERVATIONS DECLARED * 02787*10 C * IN THE MAIN PROGRAM * 02788*10 C * RJCT - FACTOR (MULTIPLIER) * 02789*10 C * * 02790**9 C * OUTPUT: * 02791**9 C * * 02792**9 C * ILOC - LOCATIONS OF THE REJECTED OBSERVATIONS IN THE * 02793*10 C * SERIES * 02794*10 C * NLOC - NUMBER OF REJECTED OBSERVATIONS * 02795*10 C * * 02796**9 C ******************************************************************* 02797**9 C 02798**9 SUBROUTINE REJECT(OBS,IOBS,RMS,RJCT,ILOC,NLOC,NLOC1) 02799**7 IMPLICIT REAL*8(A-H,O-Z) 02800**7 DIMENSION OBS(IOBS),ILOC(NLOC1) 02801**7 J=0 02802**7 RJCT=RJCT*RMS 02803**8 Y1=-RJCT 02804**8 Y2=RJCT 02805**8 DO 10 I=1,IOBS 02806**7 IF((OBS(I).GT.Y1).AND.(OBS(I).LT.Y2)) GO TO 10 02807**8 J=J+1 02808**7 ILOC(J)=I 02809**7 10 CONTINUE 02810**7 NLOC=J 02811**7 RETURN 02812**7 END 02813**7