//LEVEL JOB ' ','VERN VOGT',MSGCLASS=S /*SERVICE DEFERRED /*JOBPARM L=15,S=10,R=512 // EXEC WATFIV,RG=512K //SYSIN DD * $JOB WATFIV VOGT/U,PAGES=30 C MODIFIED BY V. VOGT SEPTEMBER 1982 C MATHEMATICAL ADJUSTMENT SECTION GEODETIC SURVEY OF CANADA C C INPUT OF DATA CARD 1 COLS 1-5 NUMBER OF SEPARATE NETS TO BE ADJUSTED CARD 2 COLS 1-72 DESCRIPTION OF JOB TO BE RUN (MAY BE LEFT BLANK) CARD 3 COLS 1-5 NUMBER OF LINKS (RIGHT JUSTIFY) C COLS 6-10 TOTAL NUMBER OF JUNCTIONS (RIGHT JUSTIFY) C COLS 11-15 NUMBER OF FIXED JUNCTIONS (RIGHT JUSTIFY) C COLS 16-20 CODE TO PRINT NORMAL EQUATIONS, ETC.ICODE = 1 CARD 4 COLS 2-10 STATION NAME OR NUMBER (JUNCTION DATA) C (TYPE) COLS 11-20 APPROXIMATE OR FIXED ELEVATION C NOTE: LIST FIXED JUNCTIONS LAST CARD 5 COLS 3-10 STATION NAME 'FROM' (LINK DATA) C (TYPE) COLS 13-20 STATION NAME 'TO' C NOTE THE NAME FIELDS MUST ALIGN COLUMNWISE WITH THE NAME FIELDS OF CARD TYPE 3 C COLS 21-30 ELEVATION DIFFERENCE OBSERVED 'TO - FROM' C COLS 31-40 INVERSE OF WEIGHT CARD 6 COLS 1-5 CONFIDENCE LEVEL FOR STATISTICAL TESTS-E.G.95.0% C PROGRAM LEVELOB C LIMITS--TOTAL OF 300 JUNCS.--250 ALLOWED FREE--ALLOWED 500 LINKS C ALL FIXED JUNCTIONS MUST BE NUMBERED LAST C B -NORMAL EQUATIONS C J -LINE NO C K -JUNCTION NUMBER C V -RESIDUALS C W -WEIGHT MATRIX C X -OBSERVED HEIGHT DIFFERENCES C EA -MATRIX OF ESTIMATED AND ADJUSTED ELEVS. OF JUNCTION POINTS C XF -VECTOR OF ADJUSTED HEIGHT DIFFERENCES C NAF -NUMBER OF FIXED STATIONS C NFS -NUMBER OF FREE STATIONS C ALPH -CONFIDENCE LEVEL OF STATISTICAL TESTING C NOBS -NUMBER OF OBSERVATIONS C SNAM -VECTOR CONTAINING NAMES OF STATIONS C NFRM -VECTOR CONTAINING STATION NUMBERS 'FROM' C NTOS -VECTOR CONTAINING STATION NUMBERS 'TO' IMPLICIT REAL*8(A-H,O-Z) DIMENSION W(500),D(250),EA(300),X(500),B(31375),C(300),SNAM(300), &NFRM(500),NTOS(500),V(500),XF(500),NTO(2),NFR(2),TITLE(10) COMMON /AA/B,D,C,NFS,NFSS,I1/BB/W,NFRM,NTOS,V/CC/SNAM 1001 FORMAT(10A8) 1002 FORMAT('1',T5,80('*'),/,T5,'*',78X,'*',/,T5,'*',3X,9A8,3X, &'*',/,T5,'*',78X,'*',/,T5,80('*'),////) 1003 FORMAT (4I5) 1004 FORMAT (T9,'INPUT DATA',//,25('*'),' FREE STATIONS ',25('*'), &///T10,'POINT',11X,'STATION',10X,'APPROXIMATE',/,T10,'NUMBER', &11X,'NAME',13X,'ELEVATION') 1005 FORMAT(A10,F10.4) 1006 FORMAT (//,25('*'),' FIXED STATIONS ',25('*'),///, &T10,'POINT',11X,'STATION',14X,'FIXED',/,T10,'NUMBER',11X, &'NAME',14X,'ELEVATION') 1007 FORMAT(T8,I5,12X,A8,2(11X,F10.4)) 1008 FORMAT(2A10,2F10.4) 1009 FORMAT (2X,13F8.4) 1010 FORMAT (100X, F10.3) 1011 FORMAT(///2X,'SOLUTION TO NORMALS',/) 1012 FORMAT(///,T25,'ADJUSTED ELEVATIONS',/,T25,19('-'),//,T10,'POINT', &12X,'STATION',12X,'ADJUSTED',16X,'STD.',/,T10,'NUMBER',12X,'NAME', &14X,'ELEVATION',15X,'DEV.',/) 1013 FORMAT(////,T5,'LINE',6X,'FROM',10X,' TO ',8X,'OBS. ELEV.',5X, &'ADJ. ELEV.',7X,' RESIDUAL',6X,'WEIGHT',5X,'LINK',/,T5,' # ', &7X,'STA.',10X,'STA.',8X,' DIFF. ',7X,' DIFF. ',9X, &'(OBS-ADJ)',17X,'S.D.',/) 1014 FORMAT(T4,I3,2(5X,A8),3(6X,F10.4),8X,F4.2,5X,F7.5) 1015 FORMAT(////12X,'MEAN RESIDUAL',17X,F10.4,/,9X,'STD.ERROR OF OBS. O &F UNIT WT.',4X, F10.4) I0=5 I1=6 READ(I0,1003) NNETS DO 999 NA = 1,NNETS C READ TITLE OF JOB READ (I0,1001) (TITLE(I),I=1,9) WRITE (I1,1002) (TITLE(I),I=1,9) READ (I0,1003) NOBS, NK, NAF, ICODE NFS = NK - NAF C NFS IS RESULTING DIMENSION OF B(N,K) WRITE (I1,1004) DO 10 J = 1,NK READ (I0,1005) SNAM(J), EA(J) IF(J.EQ.NFS+1)WRITE(I1,1006) 10 WRITE(I1,1007) J,SNAM(J), EA(J) WS = 0.0D0 DO 20 K = 1, NOBS READ(I0,1008)STO,SFR,X(K),W(K) CALL SEARCH(I,STO,NK,I1) CALL SEARCH(J,SFR,NK,I1) NFRM(K) = I NTOS (K) = J W(K) = 1.0D0 / W(K) 20 WS = WS + W(K) DO 80 K = 1, NOBS I = NFRM(K) J = NTOS (K) 80 V(K) = X(K) + EA(I) - EA(J) DO 50 I = 1, NFS C(I) = 0.0D0 NFSS = (NFS*NFS+NFS)/2 C NFSS IS THE DIMENSION OF THE SINGLE SUBSCRIPT ARRAY B--NORMAL EQNS. DO 50 J = 1,NFSS 50 B(J) = 0.0D0 DO 62 K = 1,NOBS I = NFRM(K) J = NTOS(K) II = I*NFS-(I*I-I)/2 - NFS + I JJ= J*NFS -(J*J-J)/2 - NFS + J IF(I.GT.NFS)GO TO 59 58 B(II) = W(K) + B(II) C(I) = -W(K)*V(K) + C(I) 59 IF(J.GT.NFS) GO TO 62 63 B(JJ) = W(K) +B(JJ) C(J) = W(K) * V(K) + C(J) IF(I.GT.NFS.OR.J .GT.NFS) GO TO 62 60 IF (J.GT.I) GO TO 31 30 IJ = J*NFS-(J*J-J)/2 - NFS + I GO TO 32 31 IJ=I*NFS-(I*I-I)/2-NFS+J 32 B(IJ) = -W(K) + B(IJ) 62 CONTINUE NFST = NFS JB = 1 JE = NFST IC = 1 C PRINT DESIGN MATRIX IF DESIRED 54 IF(ICODE.EQ.1)WRITE(I1,1009)(B(J),J=JB,JE) IF(ICODE.EQ.1)WRITE(I1,1010) C(IC) IF(JE .EQ.NFSS)GO TO 72 53 NFST = NFST-1 JB = JE +1 JE = JE + NFST IC = IC +1 GO TO 54 72 CALL BORDER IF(ICODE.EQ.1)WRITE (I1,1011) IF(ICODE.EQ.1)WRITE (I1,1009) (D(I), I = 1, NFS) DO 90 K = 1, NFS EA(K) = EA(K) + D(K) 90 CONTINUE WRITE (I1,1012) SUMS = 0.0D0 SUM = 0.0D0 DO 120 K = 1, NOBS I = NFRM(K) J = NTOS(K) XF(K) = EA(J) - EA(I) V(K) = X(K) - XF(K) SUM = SUM + V(K)* W(K) 120 SUMS = SUMS + V(K)**2 * W(K) XMEAN = SUM / WS S = DSQRT (SUMS / (NOBS-NFS)) DO 100 K = 1, NK IF(K .GT. NFS)GO TO 170 180 C(K) = S * DSQRT (C(K)) GO TO 100 170 C(K) = 0.0D0 100 WRITE(I1,1007) K,SNAM(K),EA(K),C(K) WRITE (I1,1013) DO 150 K = 1, NOBS I = NFRM(K) J = NTOS(K) IF(I .GT. J) GO TO 190 200 IA = I JA = J GO TO 210 190 IA = J JA = I 210 II = IA*NFS -(IA*IA-IA)/2 - NFS + IA IJ = II - IA + JA JJ = JA*NFS-(JA*JA-JA)/2 - NFS + JA IF(IA .GT. NFS .OR. JA .GT. NFS) GO TO 220 221 GO TO 250 220 IF(IA .GT. NFS) GO TO 230 GO TO 240 230 SD = 0.0D0 GO TO 150 240 SD = S * DSQRT (B(II)) GO TO 150 250 SD = S*DSQRT (B(II) + B(JJ) - 2.0D0*B(IJ)) 150 WRITE(I1,1014) K,SNAM(I), SNAM(J),X(K), XF(K),V(K),W(K),SD WRITE(I1,1015) XMEAN, S C PERFORM STATISTICAL TESTING READ(I0,1009) ALPH CALL STATS(NOBS,NFS,ALPH,S,I1) 999 CONTINUE STOP END SUBROUTINE BORDER C INVERSE MATRIX OBTAINED BY METHOD OF BORDERING IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(31375),B(250),X(250),C(300) COMMON /AA/A,B,C,NFS,NFSS,I1 C ARRAYS--A=NORMALS,B=TEMPORY,X=TEMPORARY,C=CONSTANTS C A COMES OUT AS A(INVERSE), C COMES BACK AS BM VARIANCE FACTORS C ARRAY B RETURNS AS SOLUTION C COMPUTE THE INITIAL 2X2 INVERSE I = 2 J = 2 IJ = I*NFS -(I*I-I)/2 -NFS + J F = 1.0D0/(A(1)*A(IJ) - A(2) * A(2)) E = A(1) A(1) = F*A(IJ) A(2) = -F*A(2) A(IJ) = F*E C TRANSFER NEXT COLUMN OF NORMALS INTO ARRAY B J = 3 C ARRAY B WILL BE N12 COLUMN + N22 ELEMENT 100 JL = J - 1 DO 10 I = 1, JL IJ = I*NFS -(I*I-I)/2 - NFS + J 10 B(I) = A(IJ) C ARRAY X WILL BE N11(INVERSE).BY.N12 ARRAY DO 20 I = 1, JL X(I) = 0.0 DO 20 K = 1,JL IF(I.GT.K)GO TO 22 21 IK = I*NFS -(I*I-I)/2 - NFS + K GO TO 20 22 IK = K*NFS -(K*K-K)/2 - NFS + I 20 X(I) = A(IK)*B(K) + X(I) JJ=J*NFS-(J*J-J)/2 - NFS + J TINV = A(JJ) DO 30 I = 1,JL 30 TINV = TINV -B(I)*X(I) IF(TINV) 40,50,40 40 A(JJ) = 1.0D0/TINV TINV = A(JJ) C ARRAY B WILL BE -CT(INVERSE) ARRAY 50 DO 60 I=1,JL 60 B(I) = -TINV*X(I) DO 70 I = 1,JL DO 70 K = 1, JL IF(I.GT.K)GO TO 70 71 IK = I*NFS - (I*I-I)/2 -NFS + K A(IK) = A(IK) -B(I) * X(K) 70 CONTINUE DO 80 I = 1, JL IJ = I*NFS - (I*I-I)/2 - NFS + J 80 A(IJ) = B(I) J = J+1 IF(J .GT. NFS) GO TO 90 GO TO 100 C COMPLETED INVERSE 90 DO 200 I = 1, NFS B(I) = 0.0D0 DO 200 J = 1,NFS IF(I.GT.J) GO TO 201 202 IJ = I*NFS -(I*I-I)/2- NFS + J GO TO 200 201 IJ = J*NFS - (J*J-J)/2 - NFS + I 200 B(I) = B(I) +A(IJ)*C(J) DO 300 I=1,NFS II = I*NFS -(I*I-I)/2 - NFS + I 300 C(I) = A(II) WRITE(I1,25) 25 FORMAT(//2X,'COVARIANCE MATRIX DIVIDED BY VAR. OF OBS. OF UNIT WT. &',/) DO 500 J = 1, NFS DO 400 I = 1, J IJ = I*NFS -(I*I-I)/2 - NFS + J 400 X(I) = A(IJ) WRITE(I1,5) J 500 WRITE(I1,15) (X(K), K = 1, J) 5 FORMAT(1X,'ROW/COL NO.=', I5) 15 FORMAT(1X,10F12.5) RETURN END SUBROUTINE SEARCH(IND,ANUM,N,IW) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION SNAM(300) ; COMMON/CC/SNAM DO 30 J = 1, N IF(SNAM(J) - ANUM)30,20,30 20 IND = J RETURN 30 CONTINUE C NO SUCH STATION IN LIST WRITE(IW,5) ANUM 5 FORMAT(10X,10('-'),' NO SUCH STATION AS',3X,A10,10('-')) STOP END SUBROUTINE STATS(NOBS,NFS,ALPH,VF,I1) C SUBROUTINE TO COMPUTE THE APOSTERIORI VARIANCE FACTOR, AND PERFORM C THE CHI-SQUARE TEST ON IT, AND PERFORM THE TAU-MAX TEST ON THE C RESIDUALS C TESTS DONE TO THE 'ALPHA' % CONFIDENCE LEVEL IMPLICIT REAL*8(A-H,O-Z) COMMON /BB/W,NFRM,NTOS,V/CC/SNAM REAL*4 DF,ALPH1,ALPH2,CLO,CHI DIMENSION W(500),V(500),NFRM(500),NTOS(500),SNAM(300) C COMPUTE THE ESTIMATED VARIANCE FACTOR NDF = NOBS-NFS DF = FLOAT(NDF) WRITE(I1,1) NOBS,NFS,NDF VF = VF**2 WRITE(I1,2)VF C PERFORM THE CHI-SQUARE TEST ALPH2 = (1.0D0 - ALPH/100.D0)/2.0D0 ALPH1 = 1.0D0 - ALPH2 CALL MDCHI(ALPH1,DF,CLO,IER) CALL MDCHI(ALPH2,DF,CHI,IER) CLO = DF*VF/CLO CHI = DF*VF/CHI WRITE(I1,3) CLO,CHI IF(CLO.GT.1.D0.OR.CHI.LT.1.D0)WRITE(I1,4)ALPH IF(CLO.LE.1.D0.AND.CHI.GE.1.D0)WRITE(I1,5)ALPH C COMPUTE STANDARDIZED RESIDUALS USING POPE'S APPROXIMATION C AND PERFORM TAU-MAX TEST ON THE RESIDUALS WRITE(I1,6)ALPH ALPH = 1.0D0 - ALPH/100.D0 CALL TAURE(NOBS,NDF,ALPH,CTAU) NREJ = 0 POP = DSQRT(DBLE(DF/NOBS)) DO 20 I=1,NOBS SIGV = POP * DSQRT(VF/W(I)) CRIT = SIGV * CTAU C PERFORM TAU-MAX TEST IF(DABS(V(I)).GE.CRIT)THEN DO NREJ = NREJ+1 WRITE(I1,7) SNAM(NFRM(I)),SNAM(NTOS(I)),V(I),SIGV,CRIT END IF 20 CONTINUE NPC = NREJ*100/NOBS WRITE(I1,8)NREJ,NPC 1 FORMAT('1',T48,'STATISTICS SUMMARY',/,T48,18('-'),//,20X,'NUMBER O &F OBSERVATIONS =',I12,/,20X,'NUMBER OF FREE STATIONS =',I11,/,20X, &'NUMBER OF DEGREES OF FREEDOM =',I6) 2 FORMAT(////,34X,'ESTIMATED VARIANCE FACTOR =',D15.7) 3 FORMAT(//,36X,'CHI-SQUARE TEST ON THE VARIANCE FACTOR',/,T37, &38('-'),//,28X,F12.6,6X,' < 1.000000 < ',F12.6,5X,'?') 4 FORMAT(//,24X,'TEST ON THE VARIANCE FACTOR AT THE',F5.1,' % CONFID &ENCE LEVEL FAILS',/,T84,5('-')) 5 FORMAT(//,24X,'TEST ON THE VARIANCE FACTOR AT THE',F5.1,' % CONFID &ENCE LEVEL PASSES',/,T84,6('-')) 6 FORMAT(///,24X,'SUMMARY OF REJECTION OF RESIDUALS AT THE',F5.1, &' % CONFIDENCE LEVEL',/,34X,'(TAU-MAX TEST USING POPE S APPROXIMAT &ION)'//,24X,' FROM ',10X,' TO ',10X,'RESIDUAL',10X, &'STD.DEV.',10X,'CRITICAL',/,T79,'RESIDUAL',11X,'VALUE',//) 7 FORMAT(24X,2(A8,10X),2(F8.4,10X),F8.6) 8 FORMAT(//,24X,I5,' RESIDUALS (',I3,' % OF THE OBSERVATIONS) WERE &FLAGGED FOR REJECTION',//////) RETURN END SUBROUTINE TAURE( NT,NU,ALPH,CRTAU ) C COMPUTES THE REJECTION LEVEL FOR NORMALISED RESIDUALS FOR A GIVEN NU C OBSERVATIONS , DEGREES OF FREEDOM AND DESIRED LEVEL OF TYPE I ERROR C PARAMETERS C NT - NUMBER OF OBSERVATIONS C NU - DEGREES OF FREEDOM C ALPH - DESIRED PROBABILITY OF TYPE I ERROR C CRTAU - CRITICAL VALUE ( MAX-TAU ) PRODUCED BY THE SUBROUTINE C REFERENCE : A.J. POPE (1976) - "THE STATISTICS OF RESIDUALS AND THE C DETECTION OF OUTLIERS" , U.S. DEPT. OF COMMERCE , NOAA T C REPORT NOS 65 NGS 1. IMPLICIT REAL*8(A-H,O-Z) DATA PI/ 3.1415926535898 / PD = 2. /PI S = 1. WNU = NU U = WNU -1. IF( U.EQ.0. ) GO TO 72 IF ( ALPH.EQ.0. ) GO TO 72 IF ( ALPH.LT.1. ) GO TO 10 CRTAU = 0. C RETURN C 10 Q = NT IF ( ALPH.GT.0.5 ) GO TO 19 AA = ALPH / Q DELT = AA DO 18 I = 1,100 XI = I DELT = DELT * ALPH * (( XI*Q - 1.)/(( XI+1.)*Q)) IF ( DELT.LT.1.D-14 ) GO TO 20 18 AA = AA + DELT 19 AA = 1. - (1.-ALPH)**(1./Q) 20 P = 1. - AA IF(U.EQ.1. ) GO TO 71 F = 1.3862943611199 - 2.*DLOG(AA) G = DSQRT(F) X = G - (2.515517 + 0.802853*G + 0.010328*F) & / (1. + 1.432788*G + F*(0.189269 + 0.001308*G)) Y = X*X A = X*(1. + Y)/4. B = X*(3. + Y*(16. + 5.*Y))/96. C = X*(-15. + Y*(17. + Y*(19. + 3.*Y)))/384. E = X*(-945. + Y*(-1920. + Y*(1482. + Y*(776. + 79.*Y))))/92160. V = 1./U T = X + V*(A + V*(B + V*(C + E*V))) S = T/DSQRT(U + T*T) UM = U - 1. DELL = 1. DO 70 M = 1,50 H = 1. - S*S R = 0.0 IF ( DMOD(U,2.D0).EQ.0.0 ) GO TO 49 DD = DSQRT(H) N = 0.5*UM DO 45 I = 1,N Z = 2*I R = R + DD D = DD 45 DD = DD * H * (Z/(Z+1.)) R = PD*(R*S + DARSIN(S)) D = PD*D*UM GO TO 61 49 DD = 1. N = 0.5*U DO 55 I = 1,N Z = 2*I R = R + DD D = DD 55 DD = DD*H*((Z-1.)/Z) R = R*S D = D*UM 61 CONTINUE DEL = (P-R)/D IF( DABS( DEL/DELL ) .GT.0.5) GO TO 72 DELL = DEL S = S + DEL IF( DABS(DEL) .LT. 1.D-8 ) GO TO 72 70 CONTINUE GO TO 72 71 S =DSIN(P/PD) 72 CRTAU = S*DSQRT(WNU) RETURN END $ENTRY