C PROGRAM SWSPEC C C THIS PROGRAM EVALUATES THE SPECTRUM AND C AUTOCOVARIANCE FUNCTION OF DIFFERENT TIME C SERIES CONSTRUCTED FROM THE FILE OF C SWISS LEVELLING DATA C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION T(500), F(500), DAT(1), C(2), P(100), PER(1), CO(100), $ S(500), A(500,10), D(500), PINDEX(10), X(100), COV(100), OM(500) DIMENSION INDEX(500) DIMENSION CMEAN(100),XNEW(100) C REAL*8 HEAD(8) C ABS(XX) = DABS(XX) FLOAT(II)= DFLOAT(II) C C DEFINING I-O UNITS C IRE=5 IPR = 6 C REWIND IRE C C C READ HEADER LINE OF FILE C READ(IRE,1001) (HEAD(LL),LL=1,8) 1001 FORMAT(8A8) WRITE(IPR,1002) (HEAD(LL),LL=1,8) 1002 FORMAT(1H1, 10X, 8A8) C C PREPARE ARGUMENTS FOR ANALYSIS C DO 5 J=1,10 5 PINDEX(J)=0. C C I= 0 500 CONTINUE I= I + 1 READ(IRE,1003,END=9999)DD,H,DH,DS,T1,DAT1,ST1,STA1, $ OBS1,T2,DAT2,ST2,STA2,OBS2,RL 1003 FORMAT(13X,F6.2,F6.0,F10.5,F5.0,2(1X,2F3.0,2F2.0,F3.0),F6.0) D(I) = DD A(I,1) = RL PINDEX(1) = PINDEX(1) +RL A(I,2) = DS A(I,3) = H A(I,4) = DH IF(PINDEX(1).EQ.0.00) GO TO 6 C*** A(I,5) = DH / DS PINDEX(6) = PINDEX(6) + T1 6 A(I,6) = T2 - T1 A(I,7) = 0.5* (T1+T2) PINDEX(8) = PINDEX(8) + DAT1 A(I,8) = DAT2-DAT1 A(I,9) = 0.5* (DAT1+DAT2) GO TO 500 9999 CONTINUE N = I - 1 C IF(PINDEX(1).EQ.0.0) GO TO 7 PINDEX(2) = 1.0 PINDEX(5) = 1.0 7 IF(PINDEX(6).NE.0.0) PINDEX(7)= 1.0 IF(PINDEX(8).NE.0.0) PINDEX(9)= 1.0 PINDEX(3)= 1.0 PINDEX(4)= 1.0 C C THIS IS THE END OF READING IN C ------------------------------------------------------------- NC = 100 NP = 100 C*** M = 9 M = 1 C C FOLLOWS THE MAIN PART OF PROGRAMME C THAT ORDERS THE TIME SERIES, PRODUCES C THE LEAST SQUARES SPECTRUM, PRODUCES THE C COVARIANCE FUNCTION AND PLOTS C DO 600 J= 1, M WRITE(IPR,1000) J IF(PINDEX(J).NE.0.0) GO TO 8 WRITE(IPR,1020) GO TO 600 8 CONTINUE DO 700 I= 1, N T(I) = A(I,J) 700 F(I) = D(I) C CALL SORTD(T,500,1,N,1,INDEX,1,F) RAT= T(N)-T(1) C C C MAKING SURE T(I) INCREASES MONOTONICALLY C DEL = .0001*RAT DINC = .000001*RAT DO 750 I=2,N IF(ABS(T(I)-T(I-1)).LE.DEL) T(I)= T(I-1)+DINC 750 CONTINUE C----------------------------------------------------------------------- DAT(1) = T(1) C CALL FPLOT(T,F,N,DAT,1,IPR) C C THE N-WAY ANALYSIS MAY TAP THE DATA HERE FOR C ORDERED SERIES C C C LEAST SQUARES SPECTRUM C PER(1) = 0.0 DO 800 I=1, NP 800 P(I) = FLOAT(NP-1)/(FLOAT(NP-I)/(10000*RAT)+ *FLOAT(I-1)/(0.03*RAT)) NDAT = 1 LT = 1 NPER = 0 NBASE = 0 IB = 1 NK = 2 POL = -2.0 CALL SPEC(T,F,N,FNORM,NK,DAT,NDAT,LT,PER,NPER,NBASE, *C,P,S,NP,IB,POL) C VAR = FNORM/(N-1) WRITE(IPR,1050) VAR C CALL AMPL(2,DAT,1,1,PER,0,0,C,IPR) C CALL SPLOT(P,S,NP,1,IPR) C DO 10 I=1,NC 10 X(I) = RAT*FLOAT(I)/FLOAT(NC) C C ACTUAL AUTOCOVARIANCE FUNCTION C CALL ACCOV(T,F,N,RAT,COV,NC,VAR) WRITE(IPR,1050) VAR WRITE(IPR,1060) NSTEP=10 CZERO=1.0 DAT(1) = 0.0 CALL FPLOT(X,COV,NC,DAT,1,IPR) CALL SMOOTH(X,COV,NC,CZERO,NSTEP,CMEAN,XNEW,NCMEAN) WRITE(IPR,1070) CALL FPLOT(XNEW,CMEAN,NCMEAN,DAT,1,IPR) C C COSINE TRANSFORM OF SPECTRUM C CALL COVAR(S,NP,RAT,X,CZERO,CO,NC) WRITE(IPR,1030) CALL FPLOT(X,CO,NC,DAT,1,IPR) CALL SMOOTH(X,CO,NC,CZERO,NSTEP,CMEAN,XNEW,NCMEAN) WRITE(IPR,1040) CALL FPLOT(XNEW,CMEAN,NCMEAN,DAT,1,IPR) C C 600 CONTINUE C C 1000 FORMAT(1H1,9X,15HANALYSIS FOR ,I2,13H-TH ARGUMENT ///) 1020 FORMAT(9X,35HNO DATA AVAILABLE FOR THIS ANALYSIS ,///) 1030 FORMAT(///,9X,36HFOLLOWS COSINE TRANSFORM OF SPECTRUM,///) 1040 FORMAT(///,9X,45HFOLLOWS SMOOTHED COSINE TRANSFORM OF SPECTRUM///) 1050 FORMAT(///,9X,38HTHE VARIANCE OF THE ANALYSED SERIES IS, *F10.4,///) 1060 FORMAT(///,9X,38HFOLLOWS ACTUAL AUTOCOVARIANCE FUNCTION,///) 1070 FORMAT(///,9X,47HFOLLOWS SMOOTHED ACTUAL AUTOCOVARIANCE FUNCTION,/ &//) C STOP END SUBROUTINE COVAR(S, NP, RAT, X, C0, C, NC) C C THIS SUBROUTINE COMPUTES AUTOCOVARIANCE FUNCTION C OF TIME SERIES DEFINED ON (TMIN,TMAX), RAT WIDE, C GIVEN ITS SPECTRUM S. C C ********************************************************************** C IMPLICIT REAL*8(A-H,O-Z) DIMENSION S(NP), C(NC), X(NC) COS(XX) = DCOS(XX) FLOAT(II) = DFLOAT(II) NPM1 = NP -1 C0= 0.5* S(1) DO 3 J=2, NPM1 C0= C0 + S(J) 3 CONTINUE C0 = C0 + 0.5 * S(NP) C DO 1 I=1,NC C(I)= 0.5*S(1)* COS(.00062832*X(I)/RAT) C DO 2 J=2, NPM1 OM = (0.0003*FLOAT(NP)+99.997*FLOAT(J)-100.0)/ * (RAT*3.0*(FLOAT(NP-1))) C(I)= C(I) + S(J)*COS(6.2831852*X(I)*OM) 2 CONTINUE C C(I)= (C(I)+0.5*S(NP)*COS(2094.395*X(I)/RAT)) C C(I)= C(I)/C0 1 CONTINUE C0=1.0 RETURN END SUBROUTINE ACCOV(T,F,N,RAT,BIN,NB,VAR) C C THIS SUBROUTINE COMPUTES THE ACTUAL ISOTROPIC AUTOVARIANCE C OF TIME SERIES F DEFINED ON (T(1),...,T(N)) C AND PRODUCES BIN(1),...,BIN(NB) VALUES C FOR EQUIDISTANT ARGUMENT RANGING FROM 0 TO RAT IMPLICIT REAL*8(A-H,O-Z) DIMENSION F(N),T(N),BIN(NB) FLOAT(II) = DFLOAT(II) INT(XX) = IDINT(XX) DO 1 I=1, NB 1 BIN(I) = 0.0 VAR = 0.0 DO 2 I=1,N VAR = VAR+F(I)*F(I) L = I + 1 DO 2 J=L,N K= INT(FLOAT(NB)*(T(J)-T(I))/RAT+0.9999999) BIN(K) = BIN(K) +F(I)*F(J) 2 CONTINUE DO 3 I=1,NB 3 BIN(I) = BIN(I)/VAR VAR = VAR/(N-1) RETURN END DOUBLE PRECISION FUNCTION BASE(I,T,DAT,NDAT,LT,PER,NPER) REAL*8 COS,DAT,EXP,PER,PI,SIN,T INTEGER I,IND,LT,NDAT,NPER BASE 3 DIMENSION DAT(1),PER(1) BASE 4 DATA PI/3.14159/ C BASE 6 C FUNCTION: BASE COMPUTES KNOWN CONSTITUENT FUNCTIONAL BASE 7 C VALUES. BASE 8 C BASE 9 C ARGUMENTS: I = INDEX OF KNOWN CONSTITUENT TO BE COMPUTED BASE 10 C T = TIME AT WHICH KNOWN CONSTITUENT COMPUTED BASE 11 C DAT(NDAT) = INPUT TIMES NEW DATUM BEGINS BASE 12 C LT = INPUT LINEAR TREND SWITCH (1 = INCLUDED) BASE 13 C PER(NPER) = INPUT FORCED PERIODS BASE 14 C BASE 15 C EXTERNALS: COS,EXP,SIN BASE 16 C BASE 17 C LIMITATION: USER MUST SUPPLY CODING TO COMPUTE EACH BASE 18 C USER-DEFINED CONSTITUENT. AS AN EXAMPLE, BASE 19 C THIS VERSION CONTAINS THE EXPONENTIAL BASE 20 C TREND EXP(-T/25). BASE 21 C BASE 22 C DATUM BIAS BASE 23 COS(X)=DCOS(X) EXP(X)=DEXP(X) SIN(X)=DSIN(X) IF(I .GT. NDAT) GO TO 5 BASE 24 BASE = 1. BASE 25 IF(I .EQ. NDAT .AND. T .GE. DAT(I)) RETURN BASE 26 IF(I .LT. NDAT .AND. T .GE. DAT(I) BASE 27 $ .AND. T .LT. DAT(I+1)) RETURN BASE 28 BASE = 0. BASE 29 RETURN BASE 30 C LINEAR TREND BASE 31 5 IF(I .GT. NDAT + LT) GO TO 10 BASE 32 BASE = T BASE 33 RETURN BASE 34 C FORCED PERIODS BASE 35 10 IF(I .GT. NDAT + LT + 2 * NPER) GO TO 20 BASE 36 IND = (I - NDAT - LT + 1) / 2 BASE 37 IF(I - NDAT - LT .EQ. IND * 2) GO TO 15 BASE 38 BASE = COS(2. * PI * T / PER(IND)) BASE 39 RETURN BASE 40 15 BASE = SIN(2. * PI * T / PER(IND)) BASE 41 RETURN BASE 42 C EXPONENTIAL TREND BASE 43 20 IF(I .GT. NDAT + LT + 2 * NPER + 1) GO TO 25 BASE 44 BASE = EXP(-T / 25.) BASE 45 RETURN BASE 46 C ADD ADDITIONAL USER-DEFINED FUNCTIONS HERE BASE 47 25 BASE = 0. BASE 48 RETURN BASE 49 END BASE 50 SUBROUTINE CHOLS(A,IRDA,NA) CHOL 1 REAL*8 A,EPS,ARG,ROUND,SQRT,SUM INTEGER I,IRDA,J,J2,K,NA CHOL 3 DIMENSION A(IRDA,NA) CHOL 4 DATA ROUND /500./ CHOL 5 C CHOL 6 C FUNCTION: CHOLS INVERTS MATRIX A IN PLACE CHOL 7 C USING CHOLESKY DECOMPOSITION CHOL 8 C CHOL 9 C ARGUMENTS: A(IRDA,NA) = ARRAY CONTAINING POSITIVE DEFINITE CHOL 10 C SYMMETRIC INPUT MATRIX, WITH ROW DIMENSION CHOL 11 C IRDA. THE INPUT MATRIX SIZE IS (NA,NA) CHOL 12 C AND IS INVERTED IN PLACE, DESTROYING THE CHOL 13 C INPUT, RETURNING THE INVERSE. CHOL 14 C CHOL 15 C EXTERNALS: EPS,ERROR,SQRT CHOL 16 C CHOL 17 C ERROR CONDITIONS: CHOL 18 C 101 = FATAL. DIMENSION OF A .LT. 1 CHOL 19 C 102 = FATAL. NEGATIVE SQUARE ROOT. A PROBABLY SINGULAR. CHOL 20 C 103 = FATAL. DIAGONAL ELEMENT OF CHOLESKI DECOMPOSITION CHOL 21 C NEGLIGIBLY SMALL COMPARED TO DIAGONAL ELEMENT OF A. CHOL 22 C A PROBABLY SINGULAR. CHOL 23 C ("NEGLIGIBLY SMALL" MEANS LESS THAN EPS*ROUND, CHOL 24 C WHERE EPS IS THE SMALLEST NUMBER SO THAT CHOL 25 C 1. + EPS .GT. 1., AND ROUND ACCOUNTS FOR CHOL 26 C ACCUMULATED ROUNDOFF) CHOL 27 C CHOL 28 SQRT(X)=DSQRT(X) C MATRIX DIMENSION CHECK CHOL 29 IF(NA .LT. 1) CALL ERROR(101) CHOL 30 C INVERSION OF 1X1 MATRIX CHOL 31 IF(NA .GT. 1) GO TO 5 CHOL 32 A(1,1) = 1. / A(1,1) CHOL 33 RETURN CHOL 34 C CHOLESKI DECOMPOSITION OF INPUT MATRIX CHOL 35 5 A(1,1) = SQRT(A(1,1)) CHOL 36 DO 10 I = 2,NA CHOL 37 10 A(I,1) = A(I,1) / A(1,1) CHOL 38 DO 30 J = 2,NA CHOL 39 SUM = 0. CHOL 40 DO 15 K = 2,J CHOL 41 15 SUM = SUM + A(J,K-1) ** 2 CHOL 42 IF(A(J,J) .LT. SUM) CALL ERROR(102) CHOL 43 SUM = SQRT(A(J,J) - SUM) CHOL 44 IF(SUM/A(J,J) .LT. EPS(ARG)*ROUND) CALL ERROR(103) CHOL 45 A(J,J) = SUM CHOL 46 IF(J .EQ. NA) GO TO 30 CHOL 47 J2 = J + 1 CHOL 48 DO 25 I = J2,NA CHOL 49 SUM = 0. CHOL 50 DO 20 K = 2,J CHOL 51 20 SUM = SUM + A(I,K-1) * A(J,K-1) CHOL 52 25 A(I,J) = (A(I,J) - SUM) / A(J,J) CHOL 53 30 CONTINUE CHOL 54 C INVERSION OF LOWER TRIANGULAR MATRIX CHOL 55 DO 35 I = 1,NA CHOL 56 35 A(I,I) = 1. / A(I,I) CHOL 57 DO 45 J = 2,NA CHOL 58 DO 45 I = J,NA CHOL 59 SUM = 0. CHOL 60 DO 40 K = J,I CHOL 61 40 SUM = SUM + A(I,K-1) * A(K-1,J-1) CHOL 62 45 A(I,J-1) = - A(I,I) * SUM CHOL 63 C CONSTRUCTION OF INVERSE OF INPUT MATRIX CHOL 64 DO 65 J = 1,NA CHOL 65 IF(J .EQ. 1) GO TO 55 CHOL 66 DO 50 I = 2,J CHOL 67 50 A(I-1,J) = A(J,I-1) CHOL 68 55 DO 65 I = J,NA CHOL 69 SUM = 0. CHOL 70 DO 60 K = I,NA CHOL 71 60 SUM = SUM + A(K,I) * A(K,J) CHOL 72 65 A(I,J) = SUM CHOL 73 RETURN CHOL 74 END CHOL 75 DOUBLE PRECISION FUNCTION EPS(ARG) REAL*8 ARG C EPS 3 C FUNCTION: EPS SETS FUNCTIONAL VALUE EPS AND ARGUMENT ARG EPS 4 C BOTH EQUAL TO THE SMALLEST NUMBER SO THAT EPS 5 C 1. + EPS .GT. 1. EPS 6 C EPS 7 EPS = 1. EPS 8 10 EPS = EPS / 2. EPS 9 IF((1. + EPS) - 1. .EQ. EPS) GO TO 10 EPS 10 EPS = EPS * 2. EPS 11 ARG = EPS EPS 12 RETURN EPS 13 END EPS 14 SUBROUTINE ERROR(IER) ERRO 1 INTEGER IER,IPR ERRO 2 DATA IPR /6/ ERRO 3 C ERRO 4 C FUNCTION: ERROR DETECTS WHETHER ERROR IS WARNING ERRO 5 C OR FATAL, AND PRINTS MESSAGE ERRO 6 C ERRO 7 C ARGUMENT: IER = ERROR INDEX ERRO 8 C WARNINGS HAVE INDICES 1 - 99 ERRO 9 C FATAL ERRORS HAVE INDICES 100 AND OVER. ERRO 10 C ERRO 11 IF(IER .GT. 100) GO TO 10 ERRO 12 WRITE(IPR,1001) IER ERRO 13 RETURN ERRO 14 10 WRITE(IPR,1002) IER ERRO 15 STOP ERRO 16 1001 FORMAT(11H ***WARNING,I5) ERRO 17 1002 FORMAT(15H ***FATAL ERROR,I5) ERRO 18 END ERRO 19 SUBROUTINE RESID(T,F,NF, RESI 1 $ NK,DAT,NDAT,LT,PER,NPER,NBASE, RESI 2 $ A,B,C,NKDIM) RESI 3 REAL*8 T,F,DAT,PER,A,B,C,FUNC,BASE INTEGER I,J,K,LT,NBASE,NDAT,NF,NPER,NK,NKDIM RESI 5 DIMENSION A(NKDIM,1),B(1),C(1),DAT(1),F(1), RESI 6 $ PER(1),T(1) RESI 7 C RESI 8 C FUNCTION: RESID COMPUTES THE RESIDUAL TIME SERIES RESI 9 C AFTER REMOVING THE KNOWN CONSTITUENTS RESI 10 C RESI 11 C ARGUMENTS: T(NF) = INPUT TIME SERIES TIMES RESI 12 C F(NF) = INPUT TIME SERIES VALUES RESI 13 C = OUTPUT RESIDUAL TIME SERIES VALUES RESI 14 C NK = INPUT TOTAL NUMBER OF KNOWN CONSTITUENTS RESI 15 C DAT(NDAT) = INPUT TIMES NEW DATUM BEGINS RESI 16 C LT = INPUT LINEAR TREND SWITCH (1 = INCLUDED) RESI 17 C PER(NPER) = INPUT FORCED PERIODS RESI 18 C NBASE = NUMBER OF USER-DEFINED CONSTITUENTS RESI 19 C A(NKDIM,NK) = OUTPUT ARRAY CONTAINING NORMAL RESI 20 C EQUATION COEFFICIENT MATRIX RESI 21 C B(NK) = OUTPUT NORMAL EQUATION KNOWN VECTOR RESI 22 C C(NK) = OUTPUT NORMAL EQUATION UNKNOWN VECTOR RESI 23 C RESI 24 C EXTERNALS: BASE,CHOLS RESI 25 C RESI 26 C SUMMARY: RESI 27 C CLEAR NORMAL EQUATION ARRAYS RESI 28 C CONSTRUCT NORMAL EQUATIONS FOR KNOWN CONSTITUENTS RESI 29 C INVERT NORMAL EQUATION MATRIX USING CHOLESKY ALGORITHM RESI 30 C COMPUTE SOLUTION TO NORMAL EQUATIONS RESI 31 C COMPUTE RESIDUAL TIME SERIES RESI 32 C RESI 33 DO 5 I = 1,NK RESI 34 B(I) = 0. RESI 35 DO 5 J = 1,NK RESI 36 5 A(I,J) = 0. RESI 37 DO 10 I = 1,NF RESI 38 DO 10 J = 1,NK RESI 39 FUNC = BASE(J,T(I),DAT,NDAT,LT,PER,NPER) RESI 40 B(J) = B(J) + FUNC * F(I) RESI 41 DO 10 K = J,NK RESI 42 10 A(K,J) = A(K,J) + FUNC * RESI 43 $ BASE(K,T(I),DAT,NDAT,LT,PER,NPER) RESI 44 CALL CHOLS(A,NKDIM,NK) RESI 45 DO 15 I = 1,NK RESI 46 C(I) = 0. RESI 47 DO 15 J = 1,NK RESI 48 15 C(I) = C(I) + A(I,J) * B(J) RESI 49 DO 20 I = 1,NF RESI 50 DO 20 J = 1,NK RESI 51 20 F(I) = F(I) - C(J) * RESI 52 $ BASE(J,T(I),DAT,NDAT,LT,PER,NPER) RESI 53 RETURN RESI 54 END RESI 55 SUBROUTINE AMPL(NK,DAT,NDAT,LT,PER,NPER,NBASE,C,IPR) AMPL 1 REAL*8 C,DAT,PER INTEGER I,IPR,K,K1,LT,NBASE,NDAT,NK,NPER AMPL 3 DIMENSION C(1),DAT(1),PER(1) AMPL 4 C AMPL 5 C FUNCTION: AMPL LISTS PRELIMINARY AMPLITUDES OF KNOWN AMPL 6 C CONSTITUENTS AMPL 7 C AMPL 8 C ARGUMENTS: NK = TOTAL NUMBER OF KNOWN CONSTITUENTS AMPL 9 C DAT(NDAT) = TIMES NEW DATUM BEGINS AMPL 10 C LT = LINEAR TREND SWITCH (1 = INCLUDED) AMPL 11 C PER(NPER) = FORCED PERIODS AMPL 12 C NBASE = NUMBER OF USER-DEFINED CONSTITUENTS AMPL 13 C C(NK) = PRELIMINARY AMPLITUDES OF KNOWN AMPL 14 C CONSTITUENTS AMPL 15 C IPR = UNIT NUMBER FOR OUTPUT AMPL 16 C AMPL 17 WRITE(IPR,1001) NDAT,LT,NPER,NPER,NBASE,NK AMPL 18 IF(NDAT .GE. 1) WRITE(IPR,1002) (K,C(K),K=1,NDAT) AMPL 19 K = NDAT + 1 AMPL 20 IF(LT .EQ. 1) WRITE(IPR,1003) K,C(K) AMPL 21 IF(NPER .EQ. 0) GO TO 10 AMPL 22 DO 5 I = 1,NPER AMPL 23 K = NDAT + LT + 2 * I - 1 AMPL 24 K1 = K + 1 AMPL 25 5 WRITE(IPR,1004) K,K1,PER(I),C(K),C(K1) AMPL 26 10 IF(NBASE .EQ. 0) RETURN AMPL 27 DO 15 I = 1,NBASE AMPL 28 K = NDAT + LT + 2 * NPER + I AMPL 29 15 WRITE(IPR,1005) K,C(K) AMPL 30 RETURN AMPL 31 1001 FORMAT(33H1 SOLUTION FOR KNOWN CONSTITUENTS/19X, AMPL 32 $ 5HDATUM,8X,6HLINEAR,8X,6HFORCED,8X,6HCOSINE,10X, AMPL 33 $ 4HSINE,10X,4HUSER,5X,5HTOTAL/20X,4HBIAS,9X,5HTREND, AMPL 34 $ 8X,6HPERIOD,35X,7HDEFINED//2X,7HNUMBER ,2I14,14X, AMPL 35 $ 3I14,I9//2X,9HAMPLITUDE//) AMPL 36 1002 FORMAT(7X,I3,E14.6) AMPL 37 1003 FORMAT(7X,I3,14X,E14.6) AMPL 38 1004 FORMAT(3X,I3,1H-,I3,28X,F14.6,2E14.6) AMPL 39 1005 FORMAT(7X,I3,70X,E14.6) AMPL 40 END AMPL 41 SUBROUTINE FPLOT(T,F,NF,DAT,NDAT,IPR) FPLO 1 REAL*8 AMAX1,AMIN1,DAT,F,FMAX,FMIN,T INTEGER I,IBLANK,IDAT,IFIX,IPLOT,IPR,ISTAR, FPLO 3 $ KF,NDAT,NF FPLO 4 DIMENSION DAT(1),F(1),IPLOT(100),T(1) FPLO 5 DATA IBLANK /1H /, FPLO 6 $ ISTAR /1H*/ FPLO 7 C FPLO 8 C FUNCTION: FPLOT PLOTS A TIME SERIES F FPLO 9 C DETECTING TIMES OF NEW DATUM BIASES FPLO 10 C FPLO 11 C ARGUMENTS: T(NF) = INPUT TIME SERIES TIMES FPLO 12 C F(NF) = INPUT TIME SERIES VALUES FPLO 13 C DAT(NDAT) = INPUT TIMES NEW DATUM BEGINS FPLO 14 C IPR = UNIT NUMBER FOR OUTPUT FPLO 15 C FPLO 16 C EXTERNALS: AMAX1,AMIN1,IFIX FPLO 17 C FPLO 18 C SUMMARY: FPLO 19 C INITIALIZE PLOT ARRAY FPLO 20 C COMPUTE MAXIMUM AND MINIMUM VALUES IN F FPLO 21 C SCAN THROUGH TIME SERIES FPLO 22 C CHECKING FOR NEW DATUM FPLO 23 C PLOTTING TIME SERIES VALUES FPLO 24 C FPLO 25 AMAX1(X1,X2)=DMAX1(X1,X2) AMIN(X1,X2)=DMIN1(X1,X2) WRITE(IPR,1001) FPLO 26 DO 5 I = 1,100 FPLO 27 5 IPLOT(I) = IBLANK FPLO 28 FMIN = F(1) FPLO 29 FMAX = F(1) FPLO 30 DO 10 I = 2,NF FPLO 31 FMIN = AMIN1(FMIN,F(I)) FPLO 32 10 FMAX = AMAX1(FMAX,F(I)) FPLO 33 IDAT = 1 FPLO 34 DO 20 I = 1,NF FPLO 35 IF(NDAT .EQ. 0 .OR. IDAT .GT. NDAT) GO TO 15 FPLO 36 IF(DAT(IDAT) .GT. T(I)) GO TO 15 FPLO 37 WRITE(IPR,1002) IDAT FPLO 38 IDAT = IDAT + 1 FPLO 39 15 KF= 1 + (99.* (F(I)-FMIN)/ (FMAX-FMIN)) IF(KF .LT. 1 ) KF = 1 FPLO 41 IF(KF .GT. 100) KF = 100 FPLO 42 IPLOT(KF) = ISTAR FPLO 43 WRITE(IPR,1003) I,T(I),F(I),IPLOT FPLO 44 20 IPLOT(KF) = IBLANK FPLO 45 RETURN FPLO 46 1001 FORMAT(1H1,9X,11HTIME SERIES// FPLO 47 $ 5X,1HI,8X,4HT(I),8X,4HF(I)/) FPLO 48 1002 FORMAT(30X,38(1H-),18HBEGINNING OF DATUM,I5,39(1H-)) FPLO 49 1003 FORMAT(2X,I4,2E12.4,100A1) FPLO 50 END FPLO 51 SUBROUTINE SPLOT(P,S,NW,IB,IPR) SPLO 1 REAL*8 P,S INTEGER I,IB,IBLANK,IFIX,IPLOT,IPR,ISTAR,KS,NW SPLO 3 DIMENSION IPLOT(100),P(1),S(1) SPLO 4 DATA IBLANK /1H /, SPLO 5 $ ISTAR /1H*/ SPLO 6 C SPLO 7 C FUNCTION: SPLOT PLOTS A SPECTRUM S SPLO 8 C SPLO 9 C ARGUMENTS: P(NW) = INPUT SPECTRAL PERIODS SPLO 10 C S(NW) = INPUT SPECTRAL VALUES SPLO 11 C IPR = UNIT NUMBER FOR OUTPUT SPLO 12 C SPLO 13 C EXTERNALS: IFIX SPLO 14 C SPLO 15 C SUMMARY: SPLO 16 C INITIALIZE PLOT ARRAY SPLO 17 C SCAN THROUGH SPECTRUM SPLO 18 C PLOTTING SPECTRAL VALUES (PERCENTAGE VARIANCES) SPLO 19 C SPLO 20 WRITE(IPR,1001) IB,NW,P(1),P(NW) SPLO 21 DO 5 I = 1,100 SPLO 22 5 IPLOT(I) = IBLANK SPLO 23 DO 10 I = 1,NW SPLO 24 KS=S(I) IF(KS .LT. 1) KS = 1 SPLO 26 IF(KS .GT. 100) KS = 100 SPLO 27 IPLOT(KS) = ISTAR SPLO 28 WRITE(IPR,1002) P(I),S(I),IPLOT SPLO 29 10 IPLOT(KS) = IBLANK SPLO 30 RETURN SPLO 31 1001 FORMAT(1H1,10X,13HSPECTRAL BAND,I5/ SPLO 32 $ 10X,I5,24H SPECTRAL VALUES BETWEEN,2F14.6/ SPLO 33 $ 10X,6HPERIOD,6X,19HPERCENTAGE VARIANCE/) SPLO 34 1002 FORMAT(5X,F10.5,F12.3,100A1) SPLO 35 END SPLO 36 SUBROUTINE SPEC(T,F,NF,FNORM, $ NK,DAT,NDAT,LT,PER,NPER,NBASE,C, SPEC 2 $ P,S,NW,IB,D) REAL*8 A,ABS,AMAX1,B,BASE,C,CC,COS,COSWT, $ CS,DAT,DET,EPS,EPSARG,F,FLOAT,FCOS, SPEC 5 $ FMAX,FNORM,FSIN,FUNC,OMEGA,P,PER,PI, SPEC 6 $ ROUND,S,SIN,SINWT,SQRT,SS,T,U,UAU, SPEC 7 $ UAV,V,VAV,WT,D INTEGER I,IB,NKDIM,J,K,LT,NBASE,NDAT, SPEC 9 $ NF,NK,NPER,NW SPEC 10 DIMENSION A(15,15),B(15),C(1),DAT(1),F(1),P(1), SPEC 11 $ PER(1),S(1),T(1),U(15),V(15) SPEC 12 DATA PI /3.14159265359/, SPEC 13 $ ROUND /100000./, SPEC 14 $ NKDIM /15/ SPEC 15 AMAX1(X1,X2)=DMAX1(X1,X2) COS(X)=DCOS(X) FLOAT(N)=DFLOAT(N) SIN(X)=DSIN(X) SQRT(X)=DSQRT(X) ABS(X)=DABS(X) C SPEC 16 C FUNCTION: SPEC COMPUTES DIRECT (FOR D=-2) OR INVERSE (FOR D=2) C LEAST SQUARES SPECTRUM OF C AN UNEQUALLY SPACED TIME SERIES SPEC 18 C AFTER SUPPRESSING KNOWN CONSTITUENTS SPEC 19 C SPEC 20 C ARGUMENTS: SPEC 21 C SPECIFYING THE INPUT TIME SERIES SPEC 22 C T(NF) = INPUT TIME SERIES TIMES SPEC 23 C F(NF) = INPUT TIME SERIES VALUES SPEC 24 C = OUTPUT RESIDUAL TIME SERIES VALUES SPEC 25 C FNORM = OUTPUT QUADRATIC NORM OF RESIDUAL F SPEC 26 C SPEC 27 C SPECIFYING THE KNOWN CONSTITUENTS SPEC 28 C NK = INPUT TOTAL NUMBER OF KNOWN CONSTITUENTS SPEC 29 C DAT(NDAT) = INPUT TIMES NEW DATUM BEGINS SPEC 30 C LT = INPUT LINEAR TREND SWITCH (1 = USE TREND) SPEC 31 C (0 = DO NOT USE) SPEC 32 C PER(NPER) = INPUT FORCED PERIODS SPEC 33 C NBASE = INPUT NUMBER OF USER-DEFINED CONSTITUENTS SPEC 34 C C(NK) = OUTPUT PRELIMINARY AMPLITUDES OF KNOWN SPEC 35 C CONSTITUENTS SPEC 36 C SPEC 37 C SPECIFYING THE OUTPUT SPECTRUM SPEC 38 C P(NW) = INPUT SPECTRAL PERIODS SPEC 39 C S(NW) = OUTPUT SPECTRAL VALUES SPEC 40 C IB = INPUT SPECTRAL BAND LABEL SPEC 41 C SPEC 42 C EXTERNALS: ABS,AMAX1,BASE,COS,EPS,ERROR,FLOAT,RESID, SPEC 43 C SIN,SQRT SPEC 44 C SPEC 45 C ERROR CONDITIONS: SPEC 46 C 1 = WARNING. ARGUMENT NDAT .LT. 0. (SET TO 0.) SPEC 47 C 2 = WARNING. ARGUMENT LT NOT 0 OR 1. (SET TO 0.) SPEC 48 C 3 = WARNING. ARGUMENT NPER .LT. 0. (SET TO 0.) SPEC 49 C 4 = WARNING. ARGUMENT NBASE .LT. 0. (SET TO 0.) SPEC 50 C 5 = WARNING. ARGUMENT NK .NE. NDAT+LT+2*NPER+NBASE. SPEC 51 C (SET TO NDAT + LT + 2 * NPER + NBASE.) SPEC 52 C 104 = FATAL. LESS THAN 3 TIME SERIES VALUES INPUT. SPEC 53 C 105 = FATAL. T ELEMENT VALUES NOT MONOTONIC INCREASING SPEC 54 C 106 = FATAL. NK TOO LARGE FOR DIMENSIONS OF A,B,U,V SPEC 55 C (LIMITATION NO. 2 BELOW) SPEC 56 C 107 = FATAL. DAT(1) .NE. T(1). (REQUIREMENT NO. 2 BELOW) SPEC 57 C 108 = FATAL. RESIDUAL TIME SERIES CONSISTS OF ROUNDOFF SPEC 58 C SPEC 59 C CALLING ROUTINE REQUIREMENTS: SPEC 60 C 1. WHEN NO KNOWN CONSTITUENTS ARE TO BE SUPPRESSED, THE SPEC 61 C CALLING ROUTINE MUST PASS ZERO VALUES FOR NK,NDAT, SPEC 62 C LT,NPER AND NBASE. SPEC 63 C 2. WHEN NDAT .GT. 0, THE CALLING ROUTINE MUST SET SPEC 64 C DAT(1) = T(1) SPEC 65 C 3. THE CALLING ROUTINE MUST SET SPEC 66 C NK = NDAT + LT + 2 * NPER + NBASE. SPEC 67 C 4. WHEN NBASE .GT. 0, THE USER MUST SUPPLY CODING IN SPEC 68 C FUNCTION BASE TO COMPUTE EACH USER-DEFINED SPEC 69 C CONSTITUENT. SPEC 70 C 5. ON INITIAL CALL, CALLING ROUTINE MUST SET IB = 1 TO SPEC 71 C COMPUTE RESIDUAL TIME SERIES. MANY SPECTRAL BANDS SPEC 72 C FOR THE SAME SPECTRUM CAN THEN BE COMPUTED BY SPEC 73 C SETTING IB .NE. 1, AND CALLING REPEATEDLY, SPEC 74 C CHANGING ONLY P(NW). SPEC 75 C 6. CALLING ROUTINE MUST DIMENSION ARGUMENT ARRAYS .GE. SPEC 76 C T(NF),F(NF),DAT(NDAT),PER(NPER),C(NK),P(NW),S(NW). SPEC 77 C 7. T ELEMENT VALUES ARE UNRESTRICTED AS TO SPACING, BUT SPEC 78 C MUST MONOTONICALLY INCREASE. P,DAT AND PER ELEMENT SPEC 79 C VALUES MUST BE IN THE SAME UNITS AS T. SPEC 80 C SPEC 81 C LIMITATIONS: SPEC 82 C 1. WHEN CALLED WITH IB = 1, AND NK .GT. 0, THE CONTENTS SPEC 83 C OF THE TIME SERIES F IS REPLACED BY THE RESIDUAL SPEC 84 C TIME SERIES VALUES. SPEC 85 C 2. WHEN NK .GT. NKDIM, A,B,U AND V MUST BE REDIMENSIONED SPEC 86 C .GE. NK, AND NKDIM CHANGED TO THE NEW DIMENSION. SPEC 87 C SPEC 88 IF(IB .NE. 1) GO TO 65 SPEC 89 C SPEC 90 C PROCESS INPUT ARGUMENTS SPEC 91 C CHECK NF .GE. 3 SPEC 92 C CHECK T INCREASES MONOTONICALLY SPEC 93 C COMPUTE FMAX = MAXIMUM ABSOLUTE VALUE IN F SPEC 94 C CHECK VALUES OF NDAT,LT,NPER,NBASE,AND NK SPEC 95 C CHECK DAT(1) .EQ. T(1) SPEC 96 C SPEC 97 IF(NF .LT. 3) CALL ERROR(104) SPEC 98 DO 5 I = 2,NF SPEC 99 IF(T(I) .LE. T(I-1)) CALL ERROR(105) SPEC 100 5 CONTINUE SPEC 101 FMAX = ABS(F(1)) SPEC 102 DO 10 I = 2,NF SPEC 103 10 FMAX = AMAX1(FMAX,ABS(F(I))) SPEC 104 IF(NDAT .GE. 0) GO TO 15 SPEC 105 CALL ERROR(1) SPEC 106 NDAT = 0 SPEC 107 15 IF(LT .EQ. 0 .OR. LT .EQ. 1) GO TO 20 SPEC 108 CALL ERROR(2) SPEC 109 LT = 0 SPEC 110 20 IF(NPER .GE. 0) GO TO 25 SPEC 111 CALL ERROR(3) SPEC 112 NPER = 0 SPEC 113 25 IF(NBASE .GE. 0) GO TO 30 SPEC 114 CALL ERROR(4) SPEC 115 NBASE = 0 SPEC 116 30 IF(NK .EQ. NDAT + LT + 2 * NPER + NBASE) GO TO 35 SPEC 117 CALL ERROR(5) SPEC 118 NK = NDAT + LT + 2 * NPER + NBASE SPEC 119 35 IF(NK .GT. NKDIM) CALL ERROR(106) SPEC 120 IF(NDAT .GE. 1 .AND. DAT(1) .NE. T(1)) CALL ERROR(107) SPEC 121 C SPEC 122 C SUPPRESS KNOWN CONSTITUENTS SPEC 123 C REPLACE F WITH RESIDUAL TIME SERIES SPEC 124 C COMPUTE QUADRATIC NORM OF F SPEC 125 C TERMINATE IF RMS VALUE OF RESIDUAL F IS LESS SPEC 126 C THAN EPS * FMAX * ROUND, WHERE SPEC 127 C EPS = EPSARG = SMALLEST NUMBER SO 1 + EPS .GT. 1 SPEC 128 C FMAX = MAXIMUM ABSOLUTE VALUE OF ORIGINAL F SPEC 129 C ROUND ACCOUNTS FOR ACCUMULATED ROUNDOFF IN SPEC 130 C COMPUTING RESIDUAL F SPEC 131 C SPEC 132 IF(NK .GT. 0) CALL RESID(T,F,NF, SPEC 133 $ NK,DAT,NDAT,LT,PER,NPER,NBASE, SPEC 134 $ A,B,C,NKDIM) SPEC 135 FNORM = 0. SPEC 136 DO 60 I = 1,NF SPEC 137 60 FNORM = FNORM + F(I)**2 SPEC 138 IF(SQRT(FNORM/FLOAT(NF)) .LT. SPEC 139 $ EPS(EPSARG)*FMAX*ROUND) CALL ERROR(108) SPEC 140 C SPEC 141 C FOR EACH SPECTRAL PERIOD P(I),COMPUTE SPECTRAL VALUE S(I) SPEC 142 C COMPUTE SCALAR PRODUCTS FCOS,FSIN,CC,CS,SS,U,V SPEC 143 C COMPUTE BILINEAR FORMS UAU,UAV,VAV SPEC 144 C COMPUTE PERCENTAGE VARIANCE S SPEC 145 C SPEC 146 65 DO 130 I = 1,NW SPEC 147 OMEGA = 2. * PI / P(I) SPEC 148 FCOS = 0. SPEC 149 FSIN = 0. SPEC 150 CC = 0. SPEC 151 CS = 0. SPEC 152 SS = 0. SPEC 153 IF(NK .EQ. 0) GO TO 75 SPEC 154 DO 70 J = 1,NK SPEC 155 U(J) = 0. SPEC 156 70 V(J) = 0. SPEC 157 75 DO 85 J = 1,NF SPEC 158 WT = OMEGA * T(J) SPEC 159 COSWT = COS(WT) SPEC 160 SINWT = SIN(WT) SPEC 161 FCOS = FCOS + F(J) * COSWT SPEC 162 FSIN = FSIN + F(J) * SINWT SPEC 163 CC = CC + COSWT * COSWT SPEC 164 CS = CS + COSWT * SINWT SPEC 165 SS = SS + SINWT * SINWT SPEC 166 IF(NK .EQ. 0) GO TO 85 SPEC 167 DO 80 K = 1,NK SPEC 168 FUNC = BASE(K,T(J),DAT,NDAT,LT,PER,NPER) SPEC 169 U(K) = U(K) + FUNC * COSWT SPEC 170 80 V(K) = V(K) + FUNC * SINWT SPEC 171 85 CONTINUE SPEC 172 UAU = 0. SPEC 174 UAV = 0. SPEC 175 VAV = 0. SPEC 176 IF(NK .EQ. 0) GO TO 125 SPEC 173 DO 120 J = 1,NK SPEC 177 DO 120 K = 1,NK SPEC 178 UAU = UAU + U(J) * A(J,K) * U(K) SPEC 179 UAV = UAV + U(J) * A(J,K) * V(K) SPEC 180 120 VAV = VAV + V(J) * A(J,K) * V(K) SPEC 181 125 S(I) = 0. SPEC 182 DET = (CC-UAU) * (SS-VAV) - (CS-UAV)**2 SPEC 183 IF(ABS(DET) .LT. EPSARG) GO TO 130 SPEC 184 S(I) = 100. * ((SS - VAV) * FCOS * FCOS - SPEC 185 $ D * (CS - UAV) * FCOS * FSIN + $ (CC - UAU) * FSIN * FSIN) / SPEC 187 $ (DET * FNORM) SPEC 188 130 CONTINUE SPEC 189 RETURN SPEC 190 END SPEC 191 SUBROUTINE SMOOTH(X,COV,NC,CZERO,NSTEP,CMEAN,XNEW,NCMEAN) IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(NC),COV(NC),CMEAN(NC),XNEW(NC) Z=1./NSTEP XNEW(1)=0.D0 CMEAN(1)=CZERO NN=(NC-NSTEP+1)+1 DO 1 I=2,NN CMEAN(I)=0.D0 L=I-1 K=(I+NSTEP-1)-1 DO 2 J=L,K CMEAN(I)=CMEAN(I)+COV(J) 2 CONTINUE CMEAN(I)= Z* CMEAN(I) XNEW(I)=X(I-1) 1 CONTINUE NCMEAN=NC J=NN+1 DO 3 I=J,NCMEAN XNEW(I)=X(I-1) CMEAN(I)=0.D0 3 CONTINUE RETURN END