PROGRAM TSPEC TSPE 001 IMPLICIT DOUBLE PRECISION (A-H,O-Z) TSPE 002 CHARACTER *2 EQORUN TSPE 003 CHARACTER *5 MODE TSPE 004 DIMENSION DAT(3), FF(500), PER(5), T(500) TSPE 005 DATA IPR /6/ TSPE 006 C TSPE 007 C FUNCTION: TSPEC CALLS TIMSER TO GENERATE TEST TIME SERIES TSPE 008 C AND CALLS DRIVER TO COMPUTE SEVERAL TSPE 009 C LEAST SQUARES SPECTRA OF THE TEST TSPE 010 C TIME SERIES TSPE 011 C TSPE 012 C UNIT NUMBER: IPR = 6 = LISTING OF INPUT AND OUTPUT TSPE 013 C TSPE 014 C EXTERNALS: DRIVER,FPLOT,TIMSER TSPE 015 C TSPE 016 C SUMMARY: TSPE 017 C CALL TIMSER TSPE 018 C CALL FPLOT TO PLOT INPUT TIME SERIES TSPE 019 C CALL DRIVER ADDING DATUM BIAS, LINEAR TREND, USER-DEFINED TSPE 020 C CONSTITUENTS AND FORCED FREQUENCIES SIMULTANEOUSLY (MODE=BATCH)TSPE 021 C CALL DRIVER WITH NO KNOWN CONSTITUENTS TSPE 022 C CALL DRIVER ADDING DATUM BIAS CONSTITUENTS TSPE 023 C CALL DRIVER ADDING LINEAR TREND CONSTITUENT TSPE 024 C CALL DRIVER ADDING USER-DEFINED CONSTITUENTS TSPE 025 C CALL DRIVER ADDING ONE FORCED FREQUENCY AT A TIME TSPE 026 CALL TIMSER(T, FF, NF, DAT, MDAT, MT, PER, MPER, MBASE, MODE, TSPE 027 $ EQORUN, PCENT, CLEVEL) TSPE 028 CALL FPLOT(T, FF, NF, DAT, MDAT, EQORUN, IPR) TSPE 029 IF (MODE .EQ. 'SQNTL') GO TO 1 TSPE 030 IF (MODE .EQ. 'BATCH') TSPE 031 $CALL DRIVER(T, FF, NF, DAT, MDAT, MT, PER, MPER, MBASE, IPR, TSPE 032 $ EQORUN, PCENT, CLEVEL) TSPE 033 STOP TSPE 034 1 NDAT = 0 TSPE 035 LT = 0 TSPE 036 NBASE = 0 TSPE 037 NPER = 0 TSPE 038 CALL DRIVER(T, FF, NF, DAT, NDAT, LT, PER, NPER, NBASE, IPR, TSPE 039 $ EQORUN, PCENT, CLEVEL) TSPE 040 IF(MDAT .EQ. 0) GO TO 5 TSPE 041 NDAT = MDAT TSPE 042 CALL DRIVER(T, FF, NF, DAT, NDAT, LT, PER, NPER, NBASE, IPR, TSPE 043 $ EQORUN, PCENT, CLEVEL) TSPE 044 5 IF(MT .EQ. 0) GO TO 10 TSPE 045 LT = MT TSPE 046 CALL DRIVER(T, FF, NF, DAT, NDAT, LT, PER, NPER, NBASE, IPR, TSPE 047 $ EQORUN, PCENT, CLEVEL) TSPE 048 10 IF(MBASE .EQ. 0) GO TO 15 TSPE 049 NBASE = MBASE TSPE 050 CALL DRIVER(T, FF, NF, DAT, NDAT, LT, PER, NPER, NBASE, IPR, TSPE 051 $ EQORUN, PCENT, CLEVEL) TSPE 052 15 IF(MPER .EQ. 0) GO TO 25 TSPE 053 DO 20 NPER = 1,MPER TSPE 054 20 CALL DRIVER(T, FF, NF, DAT, NDAT, LT, PER, NPER, NBASE, IPR, TSPE 055 $ EQORUN, PCENT, CLEVEL) TSPE 056 25 STOP TSPE 057 END TSPE 058 C*********************************************************************** SUBROUTINE AMPL(A, NF, NK, FNORM, DAT, NDAT, LT, PER, NPER, AMPL 001 $ NBASE, C, IPR) AMPL 002 IMPLICIT DOUBLE PRECISION (A-H,O-Z) AMPL 003 DIMENSION A(100,100), C(1), DAT(1), PER(1) AMPL 004 DATA PI/3.141592653589793D0/ AMPL 005 C AMPL 006 C FUNCTION: AMPL LISTS PRELIMINARY COSINE AND SINE AMPL 007 C COEFFICIENTS. AMPL 008 C AMPL 009 C CALLED FROM: DRIVER AMPL 010 C AMPL 011 C ARGUMENTS: A = INVERTED MATRIX OF NORMAL EQUATIONS AMPL 012 C NK = TOTAL NUMBER OF KNOWN CONSTITUENTS AMPL 013 C DAT(NDAT) = TIMES NEW DATUM BEGINS AMPL 014 C LT = LINEAR TREND SWITCH (1 = INCLUDED) AMPL 015 C PER(NPER) = FORCED PERIODS AMPL 016 C NBASE = NUMBER OF USER-DEFINED CONSTITUENTS AMPL 017 C C(NK) = PRELIMINARY AMPLITUDES OF KNOWN AMPL 018 C CONSTITUENTS AMPL 019 C IPR = UNIT NUMBER FOR OUTPUT AMPL 020 C AMPL 021 C EXTERNALS: DSQRT, DATA2, DMOD AMPL 022 WRITE(IPR,1001) NDAT,LT,NPER,NPER,NPER,NBASE AMPL 023 IF(NDAT .GE. 1) WRITE(IPR,1002) (K,C(K),K=1,NDAT) AMPL 024 K = NDAT + 1 AMPL 025 IF(LT .EQ. 1) WRITE(IPR,1003) K,C(K) AMPL 026 IF(NPER .EQ. 0) GO TO 10 AMPL 027 DO 5 I = 1, NPER AMPL 028 K = NDAT + LT + 2 * I - 1 AMPL 029 K1 = K + 1 AMPL 030 AMP = DSQRT(C(K)*C(K) + C(K1)*C(K1)) AMPL 031 GPL = DATAN2(C(K1),C(K)) * 180.D0 / PI AMPL 032 GPL = DMOD(GPL + 360.D0, 360.D0) AMPL 033 ESTSD = DSQRT(FNORM / (NF - NK)) AMPL 034 SIGAMP = DSQRT(C(K) * C(K) * A(K,K) + C(K1) * C(K1) * AMPL 035 $ A(K1,K1) + 2.D0 * C(K) * C(K1) * A(K,K1)) * AMPL 036 $ ESTSD / AMP AMPL 037 SIGPL = DSQRT(C(K1) * C(K1) * A(K,K) + AMPL 038 $ C(K) * C(K) * A(K1,K1) - AMPL 039 $ 2.D0 * C(K) * C(K1) * A(K,K1)) * AMPL 040 $ (ESTSD / AMP) * (180.D0 / PI) AMPL 041 WRITE(IPR,1004) K,K1,PER(I),C(K),C(K1),AMP,SIGAMP,GPL,SIGPL AMPL 042 5 CONTINUE AMPL 043 10 IF(NBASE .EQ. 0) RETURN AMPL 044 DO 15 I = 1,NBASE AMPL 045 K = NDAT + LT + 2 * NPER + I AMPL 046 15 WRITE(IPR,1005) K, C(K) AMPL 047 RETURN AMPL 048 1001 FORMAT(1H1,2X,31HSOLUTION FOR KNOWN CONSTITUENTS,///, AMPL 049 $ 14X,5HDATUM,4X,6HLINEAR,4X,6HFORCED,4X,6HCOSINE,6X, AMPL 050 $ 4HSINE,8X,4HUSER,/,15X,4HBIAS,5X,5HTREND, AMPL 051 $ 4X,6HPERIOD,4X,4HTERM,8X,4HTERM,5X,7HDEFINED,2X, AMPL 052 $ 9HAMPLITUDE,3X,7H(SIGMA),4X,5HPHASE,2X,7H(SIGMA),//, AMPL 053 $ 2X,6HNUMBER,5X,I5,5(7X,I3),//) AMPL 054 1002 FORMAT(7X,I3,E11.3) AMPL 055 1003 FORMAT(7X,I3,10X,E11.3) AMPL 056 1004 FORMAT(3X,I3,1H-,I3,18X,F11.3,2E11.3,10X,E10.3,1X,1H(,E10.3, AMPL 057 $ 1H),1X,F6.2,1X,1H(,F6.2,1H)) AMPL 058 1005 FORMAT(7X,I3,51X,E11.3) AMPL 059 END AMPL 060 C*********************************************************************** DOUBLE PRECISION FUNCTION BASE(I, T, DAT, NDAT, LT, PER, NPER) BASE 001 IMPLICIT DOUBLE PRECISION (A-H,O-Z) BASE 002 DIMENSION DAT(1), PER(1) BASE 003 DATA PI/3.141592653589793D0/ BASE 004 C BASE 005 C FUNCTION: BASE COMPUTES KNOWN CONSTITUENT FUNCTIONAL BASE 006 C VALUES. BASE 007 C BASE 008 C CALLED FROM: RESID BASE 009 C BASE 010 C ARGUMENTS: I = INDEX OF KNOWN CONSTITUENT TO BE COMPUTED BASE 011 C T = TIME AT WHICH KNOWN CONSTITUENT COMPUTED BASE 012 C DAT(NDAT) = INPUT TIMES NEW DATUM BEGINS BASE 013 C LT = INPUT LINEAR TREND SWITCH (1 = INCLUDED) BASE 014 C PER(NPER) = INPUT FORCED PERIODS BASE 015 C BASE 016 C EXTERNALS: DCOS, DEXP, DSIN BASE 017 C BASE 018 C LIMITATION: USER MUST SUPPLY CODING TO COMPUTE EACH BASE 019 C USER-DEFINED CONSTITUENT. AS AN EXAMPLE, BASE 020 C THIS VERSION CONTAINS THE EXPONENTIAL BASE 021 C TREND EXP(-T/25). BASE 022 C BASE 023 C DATUM BIAS BASE 024 IF(I .GT. NDAT) GO TO 5 BASE 025 BASE = 1.0D0 BASE 026 IF(I .EQ. NDAT .AND. T .GE. DAT(I)) RETURN BASE 027 IF(I .LT. NDAT .AND. T .GE. DAT(I) BASE 028 $ .AND. T .LT. DAT(I+1)) RETURN BASE 029 BASE = 0.0D0 BASE 030 RETURN BASE 031 C BASE 032 C LINEAR TREND BASE 033 5 IF(I .GT. NDAT + LT) GO TO 10 BASE 034 BASE = T BASE 035 RETURN BASE 036 C BASE 037 C FORCED PERIODS BASE 038 10 IF(I .GT. NDAT + LT + 2 * NPER) GO TO 20 BASE 039 IND = (I - NDAT - LT + 1) / 2 BASE 040 IF(I - NDAT - LT .EQ. IND * 2) GO TO 15 BASE 041 BASE = DCOS(2.D0 * PI * T / PER(IND)) BASE 042 RETURN BASE 043 15 BASE = DSIN(2.D0 * PI * T / PER(IND)) BASE 044 RETURN BASE 045 C BASE 046 C EXPONENTIAL TREND BASE 047 20 IF(I .GT. NDAT + LT + 2 * NPER + 1) GO TO 25 BASE 048 BASE = DEXP(-T / 25.D0) BASE 049 RETURN BASE 050 C BASE 051 C ADD ADDITIONAL USER-DEFINED FUNCTIONS HERE BASE 052 25 BASE = 0.D0 BASE 053 RETURN BASE 054 END BASE 055 C*********************************************************************** SUBROUTINE CHOLS(A, IRDA, NA) CHOL 001 IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHOL 002 DIMENSION A(IRDA, NA) CHOL 003 DATA ROUND /500.D0/ CHOL 004 C CHOL 005 C FUNCTION: CHOLS INVERTS MATRIX A IN PLACE CHOL 006 C USING CHOLESKY DECOMPOSITION CHOL 007 C CHOL 008 C CALLED FROM: RESID CHOL 009 C CHOL 010 C ARGUMENTS: A(IRDA,NA) = ARRAY CONTAINING POSITIVE DEFINITE CHOL 011 C SYMMETRIC INPUT MATRIX, WITH ROW DIMENSION CHOL 012 C IRDA. THE INPUT MATRIX SIZE IS (NA,NA) CHOL 013 C AND IS INVERTED IN PLACE, DESTROYING THE CHOL 014 C INPUT, RETURNING THE INVERSE. CHOL 015 C CHOL 016 C EXTERNALS: EPS, ERROR, DSQRT CHOL 017 C CHOL 018 C ERROR CONDITIONS: CHOL 019 C 101 = FATAL. DIMENSION OF A .LT. 1 CHOL 020 C 102 = FATAL. NEGATIVE SQUARE ROOT. A PROBABLY SINGULAR. CHOL 021 C 103 = FATAL. DIAGONAL ELEMENT OF CHOLESKI DECOMPOSITION CHOL 022 C NEGLIGIBLY SMALL COMPARED TO DIAGONAL ELEMENT OF A. CHOL 023 C A PROBABLY SINGULAR. CHOL 024 C ("NEGLIGIBLY SMALL" MEANS LESS THAN EPS*ROUND, CHOL 025 C WHERE EPS IS THE SMALLEST NUMBER SO THAT CHOL 026 C 1. + EPS .GT. 1., AND ROUND ACCOUNTS FOR CHOL 027 C ACCUMULATED ROUNDOFF) CHOL 028 C CHOL 029 C MATRIX DIMENSION CHECK CHOL 030 IF(NA .LT. 1) CALL ERROR(101) CHOL 031 C CHOL 032 C INVERSION OF 1X1 MATRIX CHOL 033 IF(NA .GT. 1) GO TO 5 CHOL 034 A(1,1) = 1.0D0 / A(1,1) CHOL 035 RETURN CHOL 036 C CHOL 037 C CHOLESKI DECOMPOSITION OF INPUT MATRIX CHOL 038 5 A(1,1) = DSQRT(A(1,1)) CHOL 039 DO 10 I = 2, NA CHOL 040 10 A(I,1) = A(I,1) / A(1,1) CHOL 041 DO 30 J = 2, NA CHOL 042 SUM = 0.0D0 CHOL 043 DO 15 K = 2, J CHOL 044 15 SUM = SUM + A(J,K-1) ** 2 CHOL 045 IF(A(J,J) .LT. SUM) CALL ERROR(102) CHOL 046 SUM = DSQRT(A(J,J) - SUM) CHOL 047 IF(SUM/A(J,J) .LT. EPS(ARG)*ROUND) CALL ERROR(103) CHOL 048 A(J,J) = SUM CHOL 049 IF(J .EQ. NA) GO TO 30 CHOL 050 J2 = J + 1 CHOL 051 DO 25 I = J2, NA CHOL 052 SUM = 0.0D0 CHOL 053 DO 20 K = 2, J CHOL 054 20 SUM = SUM + A(I,K-1) * A(J,K-1) CHOL 055 25 A(I,J) = (A(I,J) - SUM) / A(J,J) CHOL 056 30 CONTINUE CHOL 057 C CHOL 058 C INVERSION OF LOWER TRIANGULAR MATRIX CHOL 059 DO 35 I = 1, NA CHOL 060 35 A(I,I) = 1.0D0 / A(I,I) CHOL 061 DO 45 J = 2, NA CHOL 062 DO 45 I = J, NA CHOL 063 SUM = 0.0D0 CHOL 064 DO 40 K = J, I CHOL 065 40 SUM = SUM + A(I,K-1) * A(K-1,J-1) CHOL 066 45 A(I,J-1) = - A(I,I) * SUM CHOL 067 C CHOL 068 C CONSTRUCTION OF INVERSE OF INPUT MATRIX CHOL 069 DO 65 J = 1, NA CHOL 070 IF(J .EQ. 1) GO TO 55 CHOL 071 DO 50 I = 2, J CHOL 072 50 A(I-1,J) = A(J,I-1) CHOL 073 55 DO 65 I = J, NA CHOL 074 SUM = 0.0D0 CHOL 075 DO 60 K = I, NA CHOL 076 60 SUM = SUM + A(K,I) * A(K,J) CHOL 077 65 A(I,J) = SUM CHOL 078 RETURN CHOL 079 END CHOL 080 C*********************************************************************** SUBROUTINE COVAR(FNORM, NF, NK, A, C, PCENT, CLEVEL, IPR) COVA 001 IMPLICIT DOUBLE PRECISION (A-H,O-Z) COVA 002 DIMENSION A(100,1), COV(1,1), C(1) COVA 003 C COVA 004 C FUNCTION: COVAR COMPUTES THE VARIANCE-COVARIANCE MATRIX COVA 005 C OF THE UNKNOWN CONSTITUENTS, THE CORRELATION COVA 006 C MATRIX AND PRINTS RESULTS. COVA 007 C COVA 008 C CALLED FROM: DRIVER COVA 009 C COVA 010 C EXTERNALS: DSQRT, DABS COVA 011 C COVA 012 IF (NK .LE. 0) RETURN COVA 013 C COVA 014 C COMPUTE THE STANDARD DEVIATIONS STD, CHECK IF COVA 015 C EXCEED PCENT*C(I) AND PRINT ALL OUTSTANDING COVA 016 C STANDARD DEVIATIONS COVA 017 SIGMA2 = FNORM / (NF - NK) COVA 018 WRITE (IPR,1000) PCENT COVA 019 NSTD = 0 COVA 020 DO 25 I = 1, NK COVA 021 STD = DSQRT(SIGMA2 * A(I,I)) COVA 022 IF (STD .LT. DABS(PCENT*C(I) / 100.0D0)) GO TO 25 COVA 023 WRITE (IPR,1001) I, STD COVA 024 NSTD = NSTD + 1 COVA 025 25 CONTINUE COVA 026 IF(NSTD .EQ. 0) WRITE (IPR, 1004) COVA 027 C COVA 028 C CHECK IF ANY CORRELATION EXCEEDS CLEVEL AND PRINT COVA 029 C ALL OUTSTANDING CORRELATIONS COVA 030 WRITE(IPR,1002) CLEVEL COVA 031 NLEVEL = 0 COVA 032 DO 35 I = 1, NK COVA 033 DO 30 J = 1, NK COVA 034 IF(I .GE. J) GO TO 30 COVA 035 COR = A(I,J) / DSQRT(A(I,I) * A(J,J)) COVA 036 IF(DABS(COR) .LT. CLEVEL) GO TO 30 COVA 037 WRITE (IPR,1003) I, J, COR COVA 038 NLEVEL = NLEVEL + 1 COVA 039 30 CONTINUE COVA 040 35 CONTINUE COVA 041 IF(NLEVEL .EQ. 0) WRITE (IPR, 1004) COVA 042 1000 FORMAT(1H1, 5X, COVA 043 $ 53HOUTSTANDING STANDARD DEVIATIONS OF KNOWN CONSTITUENTS COVA 044 $ ,/,5X,12H(LARGER THAN, 1X, F5.1, 1X, COVA 045 $ 25H% OF ESTIMATED MAGNITUDE),//, 5X, 6HNUMBER, COVA 046 $ 2X, 18HSTANDARD DEVIATION,/) COVA 047 1001 FORMAT(6X, I3, 9X, E9.3) COVA 048 1002 FORMAT(1H1, 5X, COVA 049 $ 51HOUTSTANDING CORRELATIONS BETWEEN KNOWN CONSTITUENTS,/, COVA 050 $ 5X, 30H(LARGER IN ABSOLUTE VALUE THAN, 1X, F4.2,1H),//, COVA 051 $ 6X,6HNUMBER, 6X, 11HCORRELATION,/) COVA 052 1003 FORMAT(5X, I3, 1H-, I3, 5X, F11.8) COVA 053 1004 FORMAT(10X, 14HNONE WAS FOUND) COVA 054 END COVA 055 C*********************************************************************** SUBROUTINE DRIVER(T, FF, NF, DAT, NDAT, LT, PER, NPER, NBASE, IPR,DRIV 001 $ EQORUN, PCENT, CLEVEL) DRIV 002 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DRIV 003 CHARACTER *2 EQORUN DRIV 004 DIMENSION A(100,100), C(100), DAT(1), F(2000), FF(1), DRIV 005 $ P(2000), PER(1), S(2000), T(1) DRIV 006 DATA PL/ 200.D0/, DRIV 007 $ PS/ 2.D0/, DRIV 008 $ NW/125/, DRIV 009 $ IB/ 1/ DRIV 010 C DRIV 011 C FUNCTION: DRIVER CALLS SPECEQ OR SPECUN TO COMPUTE A DRIV 012 C LEAST SQUARES SPECTRUM (P,S) FOR THE INPUT DRIV 013 C TIME SERIES (T,F). DRIV 014 C DRIV 015 C CALLED FROM: TSPEC DRIV 016 C DRIV 017 C ARGUMENTS: T(NF) = INPUT TIME SERIES TIMES DRIV 018 C FF(NF) = INPUT TIME SERIES VALUES DRIV 019 C DAT(NDAT) = INPUT TIMES NEW DATUM BEGINS DRIV 020 C LT = INPUT LINEAR TREND SWITCH (1 = INCLUDED) DRIV 021 C PER(NPER) = INPUT FORCED PERIODS DRIV 022 C NBASE = NUMBER OF USER-DEFINED CONSTITUENTS DRIV 023 C IPR = UNIT NUMBER FOR OUTPUT DRIV 024 C EQORUN = FLAG FOR EQUALLY OR UNEQUALLY SPACED SERIES DRIV 025 C PCENT = PERCENTAGE LEVEL FOR DETECTING OUTSTANDING DRIV 026 C STANDARD DEVIATIONS OF UNKNOWNS DRIV 027 C CLEVEL = CRITICAL LEVEL FOR DETECTING OUTSTANDING DRIV 028 C CORRELATIONS IN THE SOLUTION DRIV 029 C DRIV 030 C EXTERNALS: AMPL, DFLOAT, SPECEQ, SPECUN, SPLOT, ERROR(108) DRIV 031 C DRIV 032 C SUMMARY: DRIV 033 C COMPUTE SPECTRAL PERIODS, P DRIV 034 C PL = LONGEST PERIOD IN P DRIV 035 C PS = SHORTEST PERIOD IN P DRIV 036 C NW = NUMBER OF PERIODS IN P DRIV 037 C COPY VECTOR F (MODIFIED BY SPECEQ AND SPECUN) DRIV 038 C COMPUTE NK = TOTAL NUMBER OF KNOWN CONSTITUENTS DRIV 039 C CALL SPECEQ OR SPECUN TO COMPUTE SPECTRUM DRIV 040 C CALL AMPL TO LIST KNOWN CONSTITUENT AMPLITUDES DRIV 041 C CALL COVAR TO LIST ALL OUTSTANDING STANDARD DEVIATIONS DRIV 042 C AND CORRELATIONS DRIV 043 C LIST RESIDUAL TIME SERIES AND ITS QUADRATIC NORM DRIV 044 C COMPUTE RS = MEAN SPECTRAL VALUE FOR WHITE NOISE DRIV 045 C COMPUTE RS95 = CRITICAL PERCENTAGE VARIANCE AT 95% DRIV 046 C CONFIDENCE LEVEL FOR DETECTING STATISTICALLY DRIV 047 C SIGNIFICANT PEAKS IN THE SPECTRUM DRIV 048 C PLOT SPECTRUM DRIV 049 DO 5 I = 1,NW DRIV 050 5 P(I) = DFLOAT(NW-1)/(DFLOAT(NW-I)/PL + DFLOAT(I-1)/PS) DRIV 051 DO 10 I = 1,NF DRIV 052 10 F(I) = FF(I) DRIV 053 NK = NDAT + LT + NBASE + 2 * NPER DRIV 054 IF(EQORUN. EQ. 'EQ') DRIV 055 $CALL SPECEQ(T, F, NF, FNORM, DRIV 056 $ NK, DAT, NDAT, LT, PER, NPER, NBASE, A, C, DRIV 057 $ P, S, NW, IB, ICRIT) DRIV 058 IF(EQORUN .EQ. 'UN') DRIV 059 $CALL SPECUN(T, F, NF, FNORM, DRIV 060 $ NK, DAT, NDAT, LT, PER, NPER, NBASE, A, C, DRIV 061 $ P, S, NW, IB, ICRIT) DRIV 062 CALL AMPL(A, NF, NK, FNORM, DAT, NDAT, LT, PER, NPER, DRIV 063 $ NBASE, C, IPR) DRIV 064 CALL COVAR(FNORM, NF, NK, A, C, PCENT, CLEVEL, IPR) DRIV 065 RS = 200.0D0 / (NF - NK) DRIV 066 RS95 = 100.0D0/(1.0D0/(0.05D0**(-2.0D0/(NF-NK-2))-1)+1) DRIV 067 WRITE(IPR,1001) (F(I),I=1,NF) DRIV 068 WRITE(IPR,1002) FNORM DRIV 069 IF(ICRIT .EQ. 0) CALL ERROR(108) DRIV 070 CALL SPLOT(P, S, NW, RS, RS95, IB, IPR) DRIV 071 RETURN DRIV 072 1001 FORMAT(1H1,9X,20HRESIDUAL TIME SERIES//110(11E10.2/)) DRIV 073 1002 FORMAT(9X,35HRESIDUAL TIME SERIES QUADRATIC NORM,E15.5) DRIV 074 END DRIV 075 C*********************************************************************** DOUBLE PRECISION FUNCTION EPS(ARG) EPS 001 IMPLICIT DOUBLE PRECISION (A-H,O-Z) EPS 002 C EPS 003 C FUNCTION: EPS SETS FUNCTIONAL VALUE EPS AND ARGUMENT ARG EPS 004 C BOTH EQUAL TO THE SMALLEST NUMBER SO THAT EPS 005 C 1. + EPS .GT. 1. EPS 006 C EPS 007 C CALLED FROM: CHOLS, SPECUN, SPECEQ EPS 008 C EPS 009 EPS = 1.0D0 EPS 010 10 EPS = EPS / 2.0D0 EPS 011 IF ((1.0D0 + EPS) - 1.0D0 .EQ. EPS) GO TO 10 EPS 012 EPS = EPS * 2.0D0 EPS 013 ARG = EPS EPS 014 RETURN EPS 015 END EPS 016 C*********************************************************************** SUBROUTINE ERROR(IER) ERRO 001 INTEGER IER,IPR ERRO 002 DATA IPR /6/ ERRO 003 C ERRO 004 C FUNCTION: ERROR DETECTS WHETHER ERROR IS WARNING ERRO 005 C OR FATAL, AND PRINTS MESSAGE ERRO 006 C ERRO 007 C CALLED FROM: DRIVER, CHOLS, SPECUN, SPECEQ ERRO 008 C ERRO 009 C ARGUMENT: IER = ERROR INDEX ERRO 010 C WARNINGS HAVE INDICES 1 - 99 ERRO 011 C FATAL ERRORS HAVE INDICES 100 AND OVER. ERRO 012 C ERRO 013 IF(IER .GT. 100) GO TO 10 ERRO 014 WRITE(IPR,1001) IER ERRO 015 RETURN ERRO 016 10 WRITE(IPR,1002) IER ERRO 017 STOP ERRO 018 1001 FORMAT(11H ***WARNING,I5) ERRO 019 1002 FORMAT(15H ***FATAL ERROR,I5) ERRO 020 END ERRO 021 C*********************************************************************** SUBROUTINE FPLOT(T, F, NF, DAT, NDAT, EQORUN, IPR) FPLO 001 IMPLICIT DOUBLE PRECISION (A-H,O-Z) FPLO 002 CHARACTER *1 IBLANK, ISTAR, IPLOT FPLO 003 CHARACTER *2 EQORUN FPLO 004 DIMENSION DAT(1), F(1), IPLOT(100), T(500) FPLO 005 DATA IBLANK /' '/, FPLO 006 $ ISTAR /'*'/ FPLO 007 C FPLO 008 C FUNCTION: FPLOT PLOTS TIME SERIES F FPLO 009 C DETECTING TIMES OF NEW DATUM BIASES FPLO 010 C FPLO 011 C CALLED FROM: TSPEC FPLO 012 C FPLO 013 C ARGUMENTS: T(NF) = INPUT TIME SERIES TIMES FPLO 014 C F(NF) = INPUT TIME SERIES VALUES FPLO 015 C DAT(NDAT) = INPUT TIMES NEW DATUM BEGINS FPLO 016 C IPR = UNIT NUMBER FOR OUTPUT FPLO 017 C FPLO 018 C EXTERNALS: DMAX1,DMIN1,IFIX FPLO 019 C FPLO 020 C SUMMARY: FPLO 021 C INITIALIZE PLOT ARRAY FPLO 022 C COMPUTE MAXIMUM AND MINIMUM VALUES IN F FPLO 023 C SCAN THROUGH TIME SERIES FPLO 024 C CHECKING FOR NEW DATUM FPLO 025 C PLOTTING TIME SERIES VALUES FPLO 026 C FPLO 027 WRITE(IPR,1001) FPLO 028 DO 5 I = 1,100 FPLO 029 5 IPLOT(I) = IBLANK FPLO 030 FMIN = F(1) FPLO 031 FMAX = F(1) FPLO 032 DO 10 I = 2,NF FPLO 033 FMIN = DMIN1(FMIN,F(I)) FPLO 034 10 FMAX = DMAX1(FMAX,F(I)) FPLO 035 IDAT = 1 FPLO 036 STEP = T(2) - T(1) FPLO 037 NGAP = 0 FPLO 038 DO 20 I = 1,NF FPLO 039 IF (I .EQ. 1 .OR. EQORUN .EQ. 'UN') GO TO 12 FPLO 040 IF (T(I)-T(I-1) .LT. 1.5D0*STEP) GO TO 12 FPLO 041 NGAP = NGAP + 1 FPLO 042 NPNT = (T(I) - T(I-1)) / STEP - 1 FPLO 043 WRITE (IPR, 1004) NGAP, NPNT FPLO 044 12 IF(NDAT .EQ. 0 .OR. IDAT .GT. NDAT) GO TO 15 FPLO 045 IF(DAT(IDAT) .GT. T(I)) GO TO 15 FPLO 046 WRITE(IPR,1002) IDAT FPLO 047 IDAT = IDAT + 1 FPLO 048 15 KF= 1 + (99.0D0 * (F(I) - FMIN) / (FMAX - FMIN)) FPLO 049 IF(KF .LT. 1 ) KF = 1 FPLO 050 IF(KF .GT. 100) KF = 100 FPLO 051 IPLOT(KF) = ISTAR FPLO 052 WRITE(IPR,1003) I,T(I),F(I),IPLOT FPLO 053 20 IPLOT(KF) = IBLANK FPLO 054 RETURN FPLO 055 1001 FORMAT(1H1,9X,11HTIME SERIES// FPLO 056 $ 5X,1HI,8X,4HT(I),8X,4HF(I)/) FPLO 057 1002 FORMAT(/,30X,38(1H-),18HBEGINNING OF DATUM,I5,39(1H-),/) FPLO 058 1003 FORMAT(2X,I4,2E12.4,100A1) FPLO 059 1004 FORMAT(/,30X,37(1H-),5HGAP #,I4,1X,2HOF,I6,1X,6HPOINTS,38(1H-)) FPLO 060 END FPLO 061 C*********************************************************************** SUBROUTINE RESID(T, F, NF, RESI 001 $ NK, DAT, NDAT, LT, PER, NPER, NBASE, RESI 002 $ A, B, C, NKDIM) RESI 003 IMPLICIT DOUBLE PRECISION (A-H,O-Z) RESI 004 DIMENSION A(NKDIM,1), B(1), C(1), DAT(1), F(1), RESI 005 $ PER(1), T(1) RESI 006 C RESI 007 C FUNCTION: RESID COMPUTES THE RESIDUAL TIME SERIES RESI 008 C AFTER REMOVING THE KNOWN CONSTITUENTS RESI 009 C RESI 010 C CALLED FROM: SPECUN, SPECEQ RESI 011 C RESI 012 C ARGUMENTS: T(NF) = INPUT TIME SERIES TIMES RESI 013 C F(NF) = INPUT TIME SERIES VALUES RESI 014 C = OUTPUT RESIDUAL TIME SERIES VALUES RESI 015 C NK = INPUT TOTAL NUMBER OF KNOWN CONSTITUENTS RESI 016 C DAT(NDAT) = INPUT TIMES NEW DATUM BEGINS RESI 017 C LT = INPUT LINEAR TREND SWITCH (1 = INCLUDED) RESI 018 C PER(NPER) = INPUT FORCED PERIODS RESI 019 C NBASE = NUMBER OF USER-DEFINED CONSTITUENTS RESI 020 C A(NKDIM,NK) = OUTPUT ARRAY CONTAINING NORMAL RESI 021 C EQUATION COEFFICIENT MATRIX RESI 022 C B(NK) = OUTPUT NORMAL EQUATION KNOWN VECTOR RESI 023 C C(NK) = OUTPUT NORMAL EQUATION UNKNOWN VECTOR RESI 024 C RESI 025 C EXTERNALS: BASE, CHOLS RESI 026 C RESI 027 C SUMMARY: RESI 028 C CLEAR NORMAL EQUATION ARRAYS RESI 029 C CONSTRUCT NORMAL EQUATIONS FOR KNOWN CONSTITUENTS RESI 030 C INVERT NORMAL EQUATION MATRIX USING CHOLESKY ALGORITHM RESI 031 C COMPUTE SOLUTION TO NORMAL EQUATIONS RESI 032 C COMPUTE RESIDUAL TIME SERIES RESI 033 C RESI 034 DO 5 I = 1, NK RESI 035 B(I) = 0.0D0 RESI 036 DO 5 J = 1, NK RESI 037 5 A(I,J) = 0.0D0 RESI 038 DO 10 I = 1, NF RESI 039 DO 10 J = 1, NK RESI 040 FUNC = BASE(J, T(I), DAT, NDAT, LT, PER, NPER) RESI 041 B(J) = B(J) + FUNC * F(I) RESI 042 DO 10 K = J, NK RESI 043 10 A(K,J) = A(K,J) + FUNC * RESI 044 $ BASE(K, T(I), DAT, NDAT, LT, PER, NPER) RESI 045 CALL CHOLS(A, NKDIM, NK) RESI 046 DO 15 I = 1, NK RESI 047 C(I) = 0.0D0 RESI 048 DO 15 J = 1, NK RESI 049 15 C(I) = C(I) + A(I,J) * B(J) RESI 050 DO 20 I = 1, NF RESI 051 DO 20 J = 1, NK RESI 052 20 F(I) = F(I) - C(J) * RESI 053 $ BASE(J, T(I), DAT, NDAT, LT, PER, NPER) RESI 054 RETURN RESI 055 END RESI 056 C*********************************************************************** SUBROUTINE SPLOT(P, S, NW, RS, RS95, IB, IPR) SPLO 001 IMPLICIT DOUBLE PRECISION (A-H,O-Z) SPLO 002 CHARACTER *1 IBLANK, ISTAR, IPLOT SPLO 003 DIMENSION IPLOT(100), P(1), S(1) SPLO 004 DATA IBLANK /' '/, SPLO 005 $ ISTAR /'*'/ SPLO 006 C SPLO 007 C FUNCTION: SPLOT PLOTS SPECTRUM S SPLO 008 C SPLO 009 C CALLED FROM: DRIVER SPLO 010 C SPLO 011 C ARGUMENTS: P(NW) = INPUT SPECTRAL PERIODS SPLO 012 C S(NW) = INPUT SPECTRAL VALUES SPLO 013 C IPR = UNIT NUMBER FOR OUTPUT SPLO 014 C SPLO 015 C EXTERNALS: IFIX SPLO 016 C SPLO 017 C SUMMARY: SPLO 018 C INITIALIZE PLOT ARRAY SPLO 019 C SCAN THROUGH SPECTRUM SPLO 020 C PLOTTING SPECTRAL VALUES (PERCENTAGE VARIANCES) SPLO 021 C SPLO 022 WRITE(IPR,1001) IB, NW, P(1), P(NW), RS, RS95 SPLO 023 DO 5 I = 1,100 SPLO 024 5 IPLOT(I) = IBLANK SPLO 025 DO 10 I = 1, NW SPLO 026 KS=S(I) SPLO 027 IF(KS .LT. 1) KS = 1 SPLO 028 IF(KS .GT. 100) KS = 100 SPLO 029 IPLOT(KS) = ISTAR SPLO 030 WRITE(IPR,1002) P(I), S(I), IPLOT SPLO 031 10 IPLOT(KS) = IBLANK SPLO 032 RETURN SPLO 033 1001 FORMAT(1H1,10X,13HSPECTRAL BAND,I5,//, SPLO 034 $ 10X,I5,24H SPECTRAL VALUES BETWEEN,2F14.6,//, SPLO 035 $ 12X,35HMEAN SPECTRAL VALUE FOR WHITE NOISE,F11.2,//, SPLO 036 $ 12X,35HCRITICAL RERCENTAGE VARIANCE AT 95%,/, SPLO 037 $ 12X,35HCONFIDENCE LEVEL FOR DETECTING,/, SPLO 038 $ 12X,35HSIGNIFICANT PEAKS IN THE SPECTRUM,F11.2,//, SPLO 039 $ 10X,6HPERIOD,6X,19HPERCENTAGE VARIANCE/) SPLO 040 1002 FORMAT(5X,F10.5,F12.3,100A1) SPLO 041 END SPLO 042 C*********************************************************************** SUBROUTINE TIMSER(T, F, NF, DAT, NDAT, LT, PER, NPER, NBASE, TIMS 001 $ MODE, EQORUN, PCENT, CLEVEL) TIMS 002 IMPLICIT DOUBLE PRECISION (A-H,O-Z) TIMS 003 CHARACTER *2 EQORUN, EQ, UN TIMS 004 CHARACTER *5 MODE, SQNTL, BATCH TIMS 005 DIMENSION X(5), Y(5), Z(5), DAT(3), F(500), IVL(4), TIMS 006 $ NB(4), NE(4), P(5), PER(5), T(500) TIMS 007 DATA X / 0.500D0, 1.000D0, 0.000D0, 0.500D0, -0.250D0/, TIMS 008 $ Y / 1.000D0, 0.500D0, 1.000D0, -0.500D0, 0.000D0/, TIMS 009 $ Z / 1.000D0, -1.000D0, 3.000D0, 0.010D0, 3.000D0/, TIMS 010 $ P / 2.759D0, 3.636D0, 5.714D0, 40.000D0, 16.000D0/, TIMS 011 $ NB / 1, 201, 281, 471/, TIMS 012 $ NE /100, 250, 400, 500/, TIMS 013 $ IVL/ 1, 2, 3, 3/, TIMS 014 $ NIVL/ 4/, TIMS 015 $ STEP/ 0.1D0/, TIMS 016 $ PI / 3.141592653589793D0/ TIMS 017 DATA EQ / 'EQ' /, TIMS 018 $ UN / 'UN' /, TIMS 019 $ SQNTL / 'SQNTL' /, TIMS 020 $ BATCH / 'BATCH' / TIMS 021 C TIMS 022 C TIMS 023 C FUNCTION: TIMSER GENERATES A TEST TIME SERIES WITH TIMS 024 C 3 DATUM BIASES, A LINEAR TREND, TIMS 025 C SIN/COSINE TERMS FOR 5 FREQUENCIES, TIMS 026 C AND AN EXPONENTIAL TREND TIMS 027 C TIMS 028 C CALLED FROM: TSPEC TIMS 029 C TIMS 030 C F(NF) = TEST TIME SERIES VALUES TIMS 031 C DAT(NDAT) = TIMES NEW DATUM BEGINS TIMS 032 C LT = LINEAR TREND SWITCH (1 = INCLUDED) TIMS 033 C PER(NPER) = PERIODS FOR TRIGONOMETRIC TERMS TIMS 034 C NBASE = NUMBER OF USER-DEFINED CONSTITUENTS TIMS 035 C EQORUN = EQUAL OR UNEQUAL SPACED TIME SERIES TIMS 036 C EQORUN = EQ: EQUAL SPACED TIME SERIES TIMS 037 C (SUBROUTINE SPECEQ IS USED) TIMS 038 C EQORUN = UN: UNEQUAL SPACED TIME SERIES TIMS 039 C (SUBROUTINE SPECUN IS USED) TIMS 040 C MODE = BATCH OR SEQUENTIAL FORCING OF UNKNOWNS TIMS 041 C MODE = SQNTL: SEQUENTIAL SOLUTION TIMS 042 C MODE = BATCH: BATCH SOLUTION TIMS 043 C PCENT = PERCENTAGE LEVEL FOR DETECTING TIMS 044 C OUTSTANDING STANDARD DEVIATIONS OF UNKNOWNS TIMS 045 C CLEVEL = CRITICAL LEVEL FOR DETECTING OUTSTANDING TIMS 046 C CORRELATIONS IN THE SOLUTION TIMS 047 C TIMS 048 C EXTERNALS: DCOS,DEXP,DFLOAT,DSIN TIMS 049 C TIMS 050 NDAT = 3 TIMS 051 LT = 1 TIMS 052 NPER = 5 TIMS 053 NBASE = 1 TIMS 054 EQORUN = EQ MODE = SQNTL TIMS 056 PCENT = 25.0D0 TIMS 057 CLEVEL = 0.50D0 TIMS 058 DO 5 I = 1,NPER TIMS 059 5 PER(I) = P(I) TIMS 060 C TIMS 061 C EQUAL SPACING RUN TIMS 062 DT1 = 0.D0 TIMS 063 DT2 = 1.D0 TIMS 064 C TIMS 065 C UNEQUAL SPACING RUN TIMS 066 C DT1 = 0.5D0 TIMS 067 C DT2 = 50.0D0 TIMS 068 NF = 0 TIMS 069 DO 10 K = 1,NIVL TIMS 070 I1 = NB(K) TIMS 071 I2 = NE(K) TIMS 072 DO 10 I = I1,I2 TIMS 073 NF = NF + 1 TIMS 074 T(NF) = DFLOAT(I) * STEP + TIMS 075 $ DT1 * DSIN(PI * DFLOAT(I) / DT2) TIMS 076 F(NF) = Z(IVL(K)) + Z(NDAT+1) * T(NF) + TIMS 077 $ Z(NDAT+2) * DEXP(-T(NF) / 25.D0) TIMS 078 DO 10 J = 1,NPER TIMS 079 10 F(NF) = F(NF) + X(J) * DCOS(2.D0 * PI * T(NF) / P(J)) TIMS 080 $ + Y(J) * DSIN(2.D0 * PI * T(NF) / P(J)) TIMS 081 DAT(1) = T(1) TIMS 082 DAT(2) = T(101) TIMS 083 DAT(3) = T(151) TIMS 084 RETURN TIMS 085 END TIMS 086 C*********************************************************************** SUBROUTINE SPECEQ(T, F, NF, FNORM, SPCQ 001 $ NK, DAT, NDAT, LT, PER, NPER, NBASE, A, C, SPCQ 002 $ P, S, NW, IB, ICRIT) SPCQ 003 IMPLICIT DOUBLE PRECISION (A-H,O-Z) SPCQ 004 DIMENSION A(100,100), B(100), C(1), CLP(60,50), SPCQ 005 $ CNP(60,50), CP(50), DAT(1), F(1), IVL(60), P(1), SPCQ 006 $ PER(1), S(1), SLP(60,50), SNP(60,50), SP(50), SPCQ 007 $ SPMQ(50), SPPQ(50), T(500), U(100), V(100), SPCQ 008 $ XL(60), XN(60) SPCQ 009 DATA PI/3.141592653589793D0/, SPCQ 010 $ ROUND /100000./, SPCQ 011 $ NKDIM /100/, SPCQ 012 $ NPERDM /50/, SPCQ 013 $ IVLDIM /60/ SPCQ 014 C SPCQ 015 C FUNCTION: SPECEQ COMPUTES THE LEAST SQUARES SPECTRUM OF SPCQ 016 C A PIECEWISE EQUALLY SPACED TIME SERIES SPCQ 017 C AFTER SUPPRESSING KNOWN CONSTITUENTS SPCQ 018 C SPCQ 019 C CALLED FROM: DRIVER SPCQ 020 C SPCQ 021 C ARGUMENTS: SPCQ 022 C SPECIFYING THE INPUT TIME SERIES SPCQ 023 C T(NF) = INPUT TIME SERIES TIMES SPCQ 024 C F(NF) = INPUT TIME SERIES VALUES SPCQ 025 C = OUTPUT RESIDUAL TIME SERIES VALUES SPCQ 026 C FNORM = OUTPUT QUADRATIC NORM OF RESIDUAL F SPCQ 027 C SPCQ 028 C SPECIFYING THE KNOWN CONSTITUENTS SPCQ 029 C NK = INPUT TOTAL NUMBER OF KNOWN CONSTITUENTS SPCQ 030 C DAT(NDAT) = INPUT TIMES NEW DATUM BEGINS SPCQ 031 C LT = INPUT LINEAR TREND SWITCH (1 = USE TREND) SPCQ 032 C (0 = DO NOT USE) SPCQ 033 C PER(NPER) = INPUT FORCED PERIODS SPCQ 034 C NBASE = INPUT NUMBER OF USER-DEFINED CONSTITUENTS SPCQ 035 C A(NKDIM,NKDIM) = OUTPUT NORMAL EQUATION MATRIX RESULTING SPCQ 036 C FROM SUPPRESSION OF KNOWN CONSTITUENTS SPCQ 037 C C(NK) = OUTPUT PRELIMINARY AMPLITUDES OF KNOWN SPCQ 038 C CONSTITUENTS SPCQ 039 C SPCQ 040 C SPECIFYING THE OUTPUT SPECTRUM SPCQ 041 C P(NW) = INPUT SPECTRAL PERIODS SPCQ 042 C S(NW) = OUTPUT SPECTRAL VALUES SPCQ 043 C IB = INPUT SPECTRAL BAND LABEL SPCQ 044 C ICRIT = ROUNDOFF FLAG SPCQ 045 C (1 = OK. CONTINUE ANALYSIS) SPCQ 046 C (0 = RESIDUAL TIME SERIES CONSISTS ONLY OF ROUNDOFF) SPCQ 047 C SPCQ 048 C EXTERNALS: DABS, DMAX1, BASE, DCOS, EPS, ERROR, DFLOAT, RESID, DSIGN, SPCQ 049 C DSIN, DSQRT SPCQ 050 C SPCQ 051 C ERROR CONDITIONS: SPCQ 052 C 1 = WARNING. ARGUMENT NDAT .LT. 0. (SET TO 0.) SPCQ 053 C 2 = WARNING. ARGUMENT LT NOT 0 OR 1. (SET TO 0.) SPCQ 054 C 3 = WARNING. ARGUMENT NPER .LT. 0. (SET TO 0.) SPCQ 055 C 4 = WARNING. ARGUMENT NBASE .LT. 0. (SET TO 0.) SPCQ 056 C 5 = WARNING. ARGUMENT NK .NE. NDAT+LT+2*NPER+NBASE. SPCQ 057 C (SET TO NDAT + LT + 2 * NPER + NBASE.) SPCQ 058 C 104 = FATAL. LESS THAN 3 TIME SERIES VALUES INPUT. SPCQ 059 C 105 = FATAL. T ELEMENT VALUES NOT MONOTONIC INCREASING SPCQ 060 C 106 = FATAL. NK TOO LARGE FOR DIMENSIONS OF A,B,U,V SPCQ 061 C (LIMITATION NO. 2 BELOW) SPCQ 062 C 107 = FATAL. DAT(1) .NE. T(1). (REQUIREMENT NO. 2 BELOW) SPCQ 063 C 108 = FATAL. RESIDUAL TIME SERIES CONSISTS OF ROUNDOFF SPCQ 064 C (NOW CALLED IN DRIVER) SPCQ 065 C 109 = FATAL. NPER TOO LARGE FOR DIMENSIONS OF SPCQ 066 C CLP,CNP,CP,SLP,SNP,SP,SPMQ,SPPQ. SPCQ 067 C 110 = FATAL. NIVL TOO LARGE FOR DIMENSIONS OF SPCQ 068 C CLP,CNP,IVL,SLP,SNP,XL,XN. SPCQ 069 C 111 = FATAL. PER CONTAINS FORCED PERIOD .LT. 2. * STEP. SPCQ 070 C 112 = FATAL. P CONTAINS SPECTRAL PERIOD .LT. 2. * STEP. SPCQ 071 C SPCQ 072 C CALLING ROUTINE REQUIREMENTS: SPCQ 073 C 1. WHEN NO KNOWN CONSTITUENTS ARE TO BE SUPPRESSED, THE SPCQ 074 C CALLING ROUTINE MUST PASS ZERO VALUES FOR NK,NDAT, SPCQ 075 C LT,NPER AND NBASE. SPCQ 076 C 2. WHEN NDAT .GT. 0, THE CALLING ROUTINE MUST SET SPCQ 077 C DAT(1) = T(1) SPCQ 078 C 3. THE CALLING ROUTINE MUST SET SPCQ 079 C NK = NDAT + LT + 2 * NPER + NBASE. SPCQ 080 C 4. WHEN NBASE .GT. 0, THE USER MUST SUPPLY CODING IN SPCQ 081 C FUNCTION BASE TO COMPUTE EACH USER-DEFINED SPCQ 082 C CONSTITUENT. SPCQ 083 C 5. ON INITIAL CALL, CALLING ROUTINE MUST SET IB = 1 TO SPCQ 084 C COMPUTE RESIDUAL TIME SERIES. MANY SPECTRAL BANDS SPCQ 085 C FOR THE SAME SPECTRUM CAN THEN BE COMPUTED BY SPCQ 086 C SETTING IB .NE. 1, AND CALLING REPEATEDLY, SPCQ 087 C CHANGING ONLY P(NW). SPCQ 088 C 6. CALLING ROUTINE MUST DIMENSION ARGUMENT ARRAYS .GE. SPCQ 089 C T(NF),F(NF),DAT(NDAT),PER(NPER),C(NK),P(NW),S(NW). SPCQ 090 C 7. T ELEMENT VALUES MUST CONSIST OF NIVL SUBINTERVALS, SPCQ 091 C EACH SUBINTERVAL CONTAINING EQUALLY SPACED DATA SPCQ 092 C SEPARATED BY A TIME INCREMENT STEP COMMON TO ALL SPCQ 093 C SUBINTERVALS. THE GAPS BETWEEN SUBINTERVALS NEED SPCQ 094 C NOT BE INTEGRAL MULTIPLES OF STEP. THE FIRST DATA SPCQ 095 C POINT MUST NOT BE ISOLATED (MUST NOT BE FOLLOWED SPCQ 096 C BY A GAP). T ELEMENT VALUES MUST INCREASE SPCQ 097 C MONOTONICALLY. DAT, P AND PER ELEMENT VALUES MUST SPCQ 098 C BE IN THE SAME UNITS AS T. SPCQ 099 C 8. THE FORCED PERIODS IN PER AND SPECTRAL PERIODS IN P SPCQ 100 C MUST BE SHORTER THAN 2 * STEP, EQUIVALENT TO THE SPCQ 101 C MAXIMUM INTERVAL USED IN FOURIER ANALYSIS. SPCQ 102 C SPCQ 103 C LIMITATIONS: SPCQ 104 C 1. WHEN CALLED WITH IB = 1, AND NK .GT. 0, THE CONTENTS SPCQ 105 C OF THE TIME SERIES F IS REPLACED BY THE RESIDUAL SPCQ 106 C TIME SERIES VALUES. SPCQ 107 C 2. WHEN NK .GT. NKDIM, A,B,U AND V MUST BE REDIMENSIONED SPCQ 108 C .GE. NK, AND NKDIM CHANGED TO THE NEW DIMENSION. SPCQ 109 C 3. WHEN NPER .GT. NPERDM THEN SP,CP,SPMQ,SPPQ AND THE SPCQ 110 C SECOND INDEX OF CLP,SLP,CNP,SNP MUST BE SPCQ 111 C REDIMENSIONED .GE. NPER, AND NPERDM CHANGED TO THE SPCQ 112 C NEW DIMENSION. SPCQ 113 C 4. WHEN NIVL (= NDAT + NUMBER OF GAPS IN T BETWEEN DATUM SPCQ 114 C CHANGES) .GT. IVLDIM THEN XN,XL,IVL AND THE FIRST SPCQ 115 C INDEX OF CLP,SLP,CNP,SNP MUST BE REDIMENSIONED .GE. SPCQ 116 C NIVL, AND IVLDIM CHANGED TO THE NEW DIMENSION. SPCQ 117 C SPCQ 118 IF(IB .NE. 1) GO TO 65 SPCQ 119 C SPCQ 120 C PROCESS INPUT ARGUMENTS SPCQ 121 C CHECK NF .GE. 3 SPCQ 122 C CHECK T INCREASES MONOTONICALLY SPCQ 123 C COMPUTE FMAX = MAXIMUM ABSOLUTE VALUE IN F SPCQ 124 C CHECK VALUES OF NDAT,LT,NPER,NBASE,AND NK SPCQ 125 C CHECK DAT(1) .EQ. T(1) SPCQ 126 IF(NF .LT. 3) CALL ERROR(104) SPCQ 127 DO 5 I = 2,NF SPCQ 128 IF(T(I) .LE. T(I-1)) CALL ERROR(105) SPCQ 129 5 CONTINUE SPCQ 130 FMAX = DABS(F(1)) SPCQ 131 DO 10 I = 2,NF SPCQ 132 10 FMAX = DMAX1(FMAX,DABS(F(I))) SPCQ 133 IF(NDAT .GE. 0) GO TO 15 SPCQ 134 CALL ERROR(1) SPCQ 135 NDAT = 0 SPCQ 136 15 IF(LT .EQ. 0 .OR. LT .EQ. 1 ) GO TO 20 SPCQ 137 CALL ERROR(2) SPCQ 138 LT = 0 SPCQ 139 20 IF(NPER .GE. 0) GO TO 25 SPCQ 140 CALL ERROR(3) SPCQ 141 NPER = 0 SPCQ 142 25 IF(NBASE .GE. 0) GO TO 30 SPCQ 143 CALL ERROR(4) SPCQ 144 NBASE = 0 SPCQ 145 30 IF(NK .EQ. NDAT + LT + 2 * NPER + NBASE) GO TO 35 SPCQ 146 CALL ERROR(5) SPCQ 147 NK = NDAT + LT + 2 * NPER + NBASE SPCQ 148 35 IF(NK .GT. NKDIM) CALL ERROR(106) SPCQ 149 IF(NDAT .GE. 1 .AND. DAT(1) .NE. T(1)) CALL ERROR(107) SPCQ 150 IF(NPER .GT. NPERDM) CALL ERROR(109) SPCQ 151 EPSARG = EPS(ARG) SPCQ 152 C SPCQ 153 C COMPUTE STEPSIZE IN T SPCQ 154 STEP = T(2) - T(1) SPCQ 155 C SPCQ 156 C COMPUTE CRITICAL VALUE OF STEP FOR DETECTING GAPS IN T SPCQ 157 STEP1 = 1.5D0 * STEP SPCQ 158 C SPCQ 159 C INITIALIZE ARGUMENTS IDAT, NIVL, NGAP, TA SPCQ 160 IDAT = 1 SPCQ 161 NIVL = 0 SPCQ 162 NGAP = 0 SPCQ 163 TA = T(1) SPCQ 164 C SPCQ 165 C FIND SUBINTERVAL BOUNDARIES (GAPS OR NEW DATUM SHIFTS) IN T SPCQ 166 DO 45 N = 2,NF SPCQ 167 NEWIVL = 0 SPCQ 168 IF(N .NE. NF) GO TO 39 SPCQ 169 NEWIVL = NGAP SPCQ 170 IF(NEWIVL .EQ. 0) GO TO 45 SPCQ 171 GO TO 42 SPCQ 172 C SPCQ 173 C CHECK IF THERE IS GAP AND DATUM SHIFT IN T SPCQ 174 39 IF((T(N) - T(N-1)) .GT. STEP1 .AND. SPCQ 175 $ DABS(DAT(IDAT+1)-T(N)) .LT. EPSARG*ROUND) GO TO 40 SPCQ 176 C SPCQ 177 C CHECK IF THERE IS ONLY GAP IN T SPCQ 178 IF((T(N) - T(N-1)) .GT. STEP1) GO TO 41 SPCQ 179 C SPCQ 180 C CHECK IF THERE IS ONLY DATUM SHIFT IN T SPCQ 181 IF (DABS(DAT(IDAT+1)-T(N)) .LT. EPSARG*ROUND) GO TO 40 SPCQ 182 IF (NEWIVL .EQ. 0) GO TO 45 SPCQ 183 GO TO 42 SPCQ 184 C SPCQ 185 C COMPUTE NEWIVL, IDAT AND NGAP WHEN THERE IS DATUM SHIFT REGARDLESS SPCQ 186 C OF PRESENCE OF GAP IN T SPCQ 187 40 NEWIVL = NGAP + 1 SPCQ 188 DAT(IDAT+1) = T(N) SPCQ 189 IDAT = IDAT + 1 SPCQ 190 GO TO 42 SPCQ 191 C SPCQ 192 C COMPUTE NEWIVL, IDAT AND NGAP WHEN THERE IS ONLY GAP IN T SPCQ 193 41 NEWIVL = NGAP + 1 SPCQ 194 NGAP = NGAP + 1 SPCQ 195 IF(NEWIVL .EQ. 0) GO TO 45 SPCQ 196 C SPCQ 197 C COMPUTE XL, XN, IVL FOR EACH SUBINTERVAL IN T SPCQ 198 42 TB = T(N-1) SPCQ 199 IF (N .EQ. NF) TB = T(NF) SPCQ 200 NIVL = NIVL + 1 SPCQ 201 IF(NIVL .GT. IVLDIM) CALL ERROR(110) SPCQ 202 XN(NIVL) = 1.0D0 + (TB - TA) / STEP SPCQ 203 XL(NIVL) = (TB + TA) / STEP SPCQ 204 IVL(NIVL) = IDAT - NEWIVL + NGAP SPCQ 205 TA = T(N) SPCQ 206 45 CONTINUE SPCQ 207 IF (NPER.LE.0) GO TO 52 SPCQ 208 C SPCQ 209 C COMPUTE TRIG(PK), TRIG(XN*PK), TRIG(XL*PK) FOR EACH SUBINTERVAL IN T SPCQ 210 DO 50 I = 1,NPER SPCQ 211 IF(PER(I) .LT. 2.0D0 * STEP) CALL ERROR(111) SPCQ 212 PK = PI * STEP / PER(I) SPCQ 213 SP(I) = DSIN(PK) SPCQ 214 CP(I) = DCOS(PK) SPCQ 215 DO 50 J = 1,NIVL SPCQ 216 XNPK = XN(J) * PK SPCQ 217 XLPK = XL(J) * PK SPCQ 218 SNP(J,I) = DSIN(XNPK) SPCQ 219 CNP(J,I) = DCOS(XNPK) SPCQ 220 SLP(J,I) = DSIN(XLPK) SPCQ 221 CLP(J,I) = DCOS(XLPK) SPCQ 222 50 CONTINUE SPCQ 223 C SPCQ 224 C CHECK VALUES IF P .GE. 2*STEP SPCQ 225 52 DO 55 I = 1,NW SPCQ 226 IF(P(I) .LT. 2.0D0 * STEP) CALL ERROR(112) SPCQ 227 55 CONTINUE SPCQ 228 C SPCQ 229 C SUPPRESS KNOWN CONSTITUENTS SPCQ 230 C REPLACE F WITH RESIDUAL TIME SERIES SPCQ 231 C COMPUTE QUADRATIC NORM OF F SPCQ 232 C CHECK IF RMS VALUE OF RESIDUAL F IS LESS SPCQ 233 C THAN EPS * FMAX * ROUND, WHERE SPCQ 234 C EPS = EPSARG = SMALLEST NUMBER SO 1 + EPS .GT. 1 SPCQ 235 C FMAX = MAXIMUM ABSOLUTE VALUE OF ORIGINAL F SPCQ 236 C ROUND ACCOUNTS FOR ACCUMULATED ROUNDOFF IN SPCQ 237 C COMPUTING RESIDUAL F SPCQ 238 IF(NK .GT. 0) CALL RESID(T, F, NF, SPCQ 239 $ NK, DAT, NDAT, LT, PER, NPER, NBASE, SPCQ 240 $ A, B, C, NKDIM) SPCQ 241 FNORM = 0.0D0 SPCQ 242 DO 60 I = 1,NF SPCQ 243 FNORM = FNORM + F(I) ** 2 SPCQ 244 60 CONTINUE SPCQ 245 C SPCQ 246 C CHECK IF RESIDUAL F CONSISTS OF ROUNDOFF SPCQ 247 ICRIT = 1 SPCQ 248 IF(DSQRT(FNORM/DFLOAT(NF)) .LT. SPCQ 249 $ EPSARG*FMAX*ROUND) ICRIT = 0 SPCQ 250 C SPCQ 251 C FOR EACH SPECTRAL PERIOD P(I),COMPUTE SPECTRAL VALUE S(I) SPCQ 252 C COMPUTE SCALAR PRODUCTS FCOS,FSIN,CC,CS,SS,U,V SPCQ 253 C COMPUTE BILINEAR FORMS UAU,UAV,VAV SPCQ 254 C COMPUTE PERCENTAGE VARIANCE S SPCQ 255 65 DO 130 I = 1,NW SPCQ 256 OMEGA = 2.0D0 * PI / P(I) SPCQ 257 FCOS = 0.0D0 SPCQ 258 FSIN = 0.0D0 SPCQ 259 CC = 0.0D0 SPCQ 260 CS = 0.0D0 SPCQ 261 SS = 0.0D0 SPCQ 262 IF(NK .EQ. 0) GO TO 75 SPCQ 263 DO 70 J = 1,NK SPCQ 264 U(J) = 0.0D0 SPCQ 265 V(J) = 0.0D0 SPCQ 266 70 CONTINUE SPCQ 267 75 DO 85 J = 1,NF SPCQ 268 WT = OMEGA * T(J) SPCQ 269 COSWT = DCOS(WT) SPCQ 270 SINWT = DSIN(WT) SPCQ 271 FCOS = FCOS + F(J) * COSWT SPCQ 272 FSIN = FSIN + F(J) * SINWT SPCQ 273 IF(NBASE .EQ. 0) GO TO 85 SPCQ 274 DO 80 L = 1, NBASE SPCQ 275 K = NDAT + LT + 2 * NPER + L SPCQ 276 FUNC = BASE(K, T(J), DAT, NDAT, LT, PER, NPER) SPCQ 277 U(K) = U(K) + FUNC * COSWT SPCQ 278 V(K) = V(K) + FUNC * SINWT SPCQ 279 80 CONTINUE SPCQ 280 85 CONTINUE SPCQ 281 Q = 0.5D0 * OMEGA * STEP SPCQ 282 SQ = DSIN(Q) SPCQ 283 CQ = DCOS(Q) SPCQ 284 IF(NPER .EQ. 0) GO TO 95 SPCQ 285 DO 90 J = 1,NPER SPCQ 286 SPMQ(J) = SP(J) * CQ - CP(J) * SQ SPCQ 287 IF(DABS(SPMQ(J)) .LT. EPSARG) SPMQ(J) = DSIGN(EPSARG,SPMQ(J)) SPCQ 288 SPPQ(J) = SP(J) * CQ + CP(J) * SQ SPCQ 289 IF(DABS(SPPQ(J)) .LT. EPSARG) SPPQ(J) = DSIGN(EPSARG,SPPQ(J)) SPCQ 290 90 CONTINUE SPCQ 291 95 DO 115 J = 1,NIVL SPCQ 292 XNQ = XN(J) * Q SPCQ 293 XLQ = XL(J) * Q SPCQ 294 SNQ = DSIN(XNQ) SPCQ 295 CNQ = DCOS(XNQ) SPCQ 296 SLQ = DSIN(XLQ) SPCQ 297 CLQ = DCOS(XLQ) SPCQ 298 CC = CC + SNQ * CNQ * CLQ * CLQ - SNQ * CNQ * SLQ * SLQ SPCQ 299 CS = CS + SNQ * CNQ * SLQ * CLQ SPCQ 300 IF(NK .EQ. 0) GO TO 115 SPCQ 301 IF(NDAT .EQ. 0) GO TO 100 SPCQ 302 K = IVL(J) SPCQ 303 U(K) = U(K) + SNQ * CLQ / SQ SPCQ 304 V(K) = V(K) + SNQ * SLQ / SQ SPCQ 305 100 IF(LT .EQ. 0) GO TO 105 SPCQ 306 K = NDAT + 1 SPCQ 307 SSCS = - SNQ * SLQ * CQ / SQ SPCQ 308 SCCS = SNQ * CLQ * CQ / SQ SPCQ 309 XNCS = XN(J) * CNQ * SLQ SPCQ 310 XNCC = - XN(J) * CNQ * CLQ SPCQ 311 XLSC = XL(J) * SNQ * CLQ SPCQ 312 XLSS = XL(J) * SNQ * SLQ SPCQ 313 STSQ = 0.5D0 * STEP / SQ SPCQ 314 U(K) = U(K) + STSQ * (SSCS + XNCS + XLSC) SPCQ 315 V(K) = V(K) + STSQ * (SCCS + XNCC + XLSS) SPCQ 316 105 IF(NPER .EQ. 0) GO TO 115 SPCQ 317 DO 110 L = 1,NPER SPCQ 318 K= NDAT + LT + 2 * L - 1 SPCQ 319 SNPMQ = (SNP(J,L) * CNQ - CNP(J,L) * SNQ) / SPMQ(L) SPCQ 320 SNPPQ = (SNP(J,L) * CNQ + CNP(J,L) * SNQ) / SPPQ(L) SPCQ 321 SINM = (SLP(J,L) * CLQ - CLP(J,L) * SLQ) * SNPMQ SPCQ 322 SINP = (SLP(J,L) * CLQ + CLP(J,L) * SLQ) * SNPPQ SPCQ 323 COSM = (CLP(J,L) * CLQ + SLP(J,L) * SLQ) * SNPMQ SPCQ 324 COSP = (CLP(J,L) * CLQ - SLP(J,L) * SLQ) * SNPPQ SPCQ 325 U(K) = U(K) + 0.5D0 * ( COSP + COSM) SPCQ 326 U(K+1) = U(K+1) + 0.5D0 * ( SINP + SINM) SPCQ 327 V(K) = V(K) + 0.5D0 * ( SINP - SINM) SPCQ 328 V(K+1) = V(K+1) + 0.5D0 * (-COSP + COSM) SPCQ 329 110 CONTINUE SPCQ 330 115 CONTINUE SPCQ 331 SQCQ = SQ * CQ SPCQ 332 CCSQCQ = CC / SQCQ SPCQ 333 SS = 0.5D0 * (DFLOAT(NF) - CCSQCQ) SPCQ 334 CC = 0.5D0 * (DFLOAT(NF) + CCSQCQ) SPCQ 335 CS = CS / SQCQ SPCQ 336 UAU = 0.0D0 SPCQ 337 UAV = 0.0D0 SPCQ 338 VAV = 0.0D0 SPCQ 339 IF(NK .EQ. 0) GO TO 125 SPCQ 340 DO 120 J = 1,NK SPCQ 341 DO 120 K = 1,NK SPCQ 342 UAU = UAU + U(J) * A(J,K) * U(K) SPCQ 343 UAV = UAV + U(J) * A(J,K) * V(K) SPCQ 344 VAV = VAV + V(J) * A(J,K) * V(K) SPCQ 345 120 CONTINUE SPCQ 346 125 S(I) = 0.0D0 SPCQ 347 DET = (CC-UAU) * (SS-VAV) - (CS-UAV) * (CS-UAV) SPCQ 348 IF(DABS(DET) .LT. EPSARG) GO TO 130 SPCQ 349 S(I) = 100.0D0 * ((SS - VAV) * FCOS * FCOS - SPCQ 350 $ 2.0D0 * (CS - UAV) * FCOS * FSIN + SPCQ 351 $ (CC - UAU) * FSIN * FSIN) / SPCQ 352 $ (DET * FNORM) SPCQ 353 130 CONTINUE SPCQ 354 RETURN SPCQ 355 END SPCQ 356 C*********************************************************************** SUBROUTINE SPECUN(T, F, NF, FNORM, SPUN 001 $ NK, DAT, NDAT, LT, PER, NPER, NBASE, A, C, SPUN 002 $ P, S, NW, IB, ICRIT) SPUN 003 IMPLICIT DOUBLE PRECISION (A-H,O-Z) SPUN 004 DIMENSION A(100,100), B(100), C(1), DAT(1), F(1), P(1), SPUN 005 $ PER(1), S(1), T(1), U(100), V(100) SPUN 006 DATA PI /3.141592653589793D0/, SPUN 007 $ ROUND /100000./, SPUN 008 $ NKDIM /100/ SPUN 009 C SPUN 010 C FUNCTION: SPECUN COMPUTES THE LEAST SQUARES SPECTRUM OF SPUN 011 C AN UNEQUALLY SPACED TIME SERIES SPUN 012 C AFTER SUPPRESSING KNOWN CONSTITUENTS SPUN 013 C SPUN 014 C CALLED FROM: DRIVER SPUN 015 C SPUN 016 C ARGUMENTS: SPUN 017 C SPECIFYING THE INPUT TIME SERIES SPUN 018 C T(NF) = INPUT TIME SERIES TIMES SPUN 019 C F(NF) = INPUT TIME SERIES VALUES SPUN 020 C = OUTPUT RESIDUAL TIME SERIES VALUES SPUN 021 C FNORM = OUTPUT QUADRATIC NORM OF RESIDUAL F SPUN 022 C SPUN 023 C SPECIFYING THE KNOWN CONSTITUENTS SPUN 024 C NK = INPUT TOTAL NUMBER OF KNOWN CONSTITUENTS SPUN 025 C DAT(NDAT) = INPUT TIMES NEW DATUM BEGINS SPUN 026 C LT = INPUT LINEAR TREND SWITCH (1 = USE TREND) SPUN 027 C (0 = DO NOT USE) SPUN 028 C PER(NPER) = INPUT FORCED PERIODS SPUN 029 C NBASE = INPUT NUMBER OF USER-DEFINED CONSTITUENTS SPUN 030 C A(NKDIM,NKDIM) = OUTPUT NORMAL EQUATION MATRIX RESULTING SPUN 031 C FROM SUPPRESSION OF KNOWN CONSTITUENTS SPUN 032 C C(NK) = OUTPUT PRELIMINARY AMPLITUDES OF KNOWN SPUN 033 C CONSTITUENTS SPUN 034 C SPUN 035 C SPECIFYING THE OUTPUT SPECTRUM SPUN 036 C P(NW) = INPUT SPECTRAL PERIODS SPUN 037 C DSIN,DSQRT SPUN 038 C S(NW) = OUTPUT SPECTRAL VALUES SPUN 039 C IB = INPUT SPECTRAL BAND LABEL SPUN 040 C ICRIT = OUTPUT ROUNDOFF FLAG SPUN 041 C (1 = OK. CONTINUE ANALYSIS) SPUN 042 C (0 = RESIDUAL TIME SERIES CONSISTS ONLY OF ROUNDOFF) SPUN 043 C SPUN 044 C EXTERNALS: DABS, DMAX1, BASE, DCOS, EPS, ERROR, DFLOAT, RESID, SPUN 045 C SPUN 046 C ERROR CONDITIONS: SPUN 047 C 1 = WARNING. ARGUMENT NDAT .LT. 0. (SET TO 0.) SPUN 048 C 2 = WARNING. ARGUMENT LT NOT 0 OR 1. (SET TO 0.) SPUN 049 C 3 = WARNING. ARGUMENT NPER .LT. 0. (SET TO 0.) SPUN 050 C 4 = WARNING. ARGUMENT NBASE .LT. 0. (SET TO 0.) SPUN 051 C 5 = WARNING. ARGUMENT NK .NE. NDAT+LT+2*NPER+NBASE. SPUN 052 C (SET TO NDAT + LT + 2 * NPER + NBASE.) SPUN 053 C 104 = FATAL. LESS THAN 3 TIME SERIES VALUES INPUT. SPUN 054 C 105 = FATAL. T ELEMENT VALUES NOT MONOTONIC INCREASING SPUN 055 C 106 = FATAL. NK TOO LARGE FOR DIMENSIONS OF A,B,U,V SPUN 056 C (LIMITATION NO. 2 BELOW) SPUN 057 C 107 = FATAL. DAT(1) .NE. T(1). (REQUIREMENT NO. 2 BELOW) SPUN 058 C 108 = FATAL. RESIDUAL TIME SERIES CONSISTS OF ROUNDOFF SPUN 059 C (NOW CALLED IN DRIVER) SPUN 060 C SPUN 061 C CALLING ROUTINE REQUIREMENTS: SPUN 062 C 1. WHEN NO KNOWN CONSTITUENTS ARE TO BE SUPPRESSED, THE SPUN 063 C CALLING ROUTINE MUST PASS ZERO VALUES FOR NK,NDAT, SPUN 064 C LT,NPER AND NBASE. SPUN 065 C 2. WHEN NDAT .GT. 0, THE CALLING ROUTINE MUST SET SPUN 066 C DAT(1) = T(1) SPUN 067 C 3. THE CALLING ROUTINE MUST SET SPUN 068 C NK = NDAT + LT + 2 * NPER + NBASE. SPUN 069 C 4. WHEN NBASE .GT. 0, THE USER MUST SUPPLY CODING IN SPUN 070 C FUNCTION BASE TO COMPUTE EACH USER-DEFINED SPUN 071 C CONSTITUENT. SPUN 072 C 5. ON INITIAL CALL, CALLING ROUTINE MUST SET IB = 1 TO SPUN 073 C COMPUTE RESIDUAL TIME SERIES. MANY SPECTRAL BANDS SPUN 074 C FOR THE SAME SPECTRUM CAN THEN BE COMPUTED BY SPUN 075 C SETTING IB .NE. 1, AND CALLING REPEATEDLY, SPUN 076 C CHANGING ONLY P(NW). SPUN 077 C 6. CALLING ROUTINE MUST DIMENSION ARGUMENT ARRAYS .GE. SPUN 078 C T(NF),F(NF),DAT(NDAT),PER(NPER),C(NK),P(NW),S(NW). SPUN 079 C 7. T ELEMENT VALUES ARE UNRESTRICTED AS TO SPACING, BUT SPUN 080 C MUST MONOTONICALLY INCREASE. P,DAT AND PER ELEMENT SPUN 081 C VALUES MUST BE IN THE SAME UNITS AS T. SPUN 082 C SPUN 083 C LIMITATIONS: SPUN 084 C 1. WHEN CALLED WITH IB = 1, AND NK .GT. 0, THE CONTENTS SPUN 085 C OF THE TIME SERIES F IS REPLACED BY THE RESIDUAL SPUN 086 C TIME SERIES VALUES. SPUN 087 C 2. WHEN NK .GT. NKDIM, A,B,U AND V MUST BE REDIMENSIONED SPUN 088 C .GE. NK, AND NKDIM CHANGED TO THE NEW DIMENSION. SPUN 089 C SPUN 090 IF(IB .NE. 1) GO TO 65 SPUN 091 C SPUN 092 C PROCESS INPUT ARGUMENTS SPUN 093 C CHECK NF .GE. 3 SPUN 094 C CHECK T INCREASES MONOTONICALLY SPUN 095 C COMPUTE FMAX = MAXIMUM ABSOLUTE VALUE IN F SPUN 096 C CHECK VALUES OF NDAT,LT,NPER,NBASE,AND NK SPUN 097 C CHECK DAT(1) .EQ. T(1) SPUN 098 IF(NF .LT. 3) CALL ERROR(104) SPUN 099 DO 5 I = 2, NF SPUN 100 IF(T(I) .LE. T(I-1)) CALL ERROR(105) SPUN 101 5 CONTINUE SPUN 102 FMAX = DABS(F(1)) SPUN 103 DO 10 I = 2, NF SPUN 104 10 FMAX = DMAX1(FMAX, DABS(F(I))) SPUN 105 IF(NDAT .GE. 0) GO TO 15 SPUN 106 CALL ERROR(1) SPUN 107 NDAT = 0 SPUN 108 15 IF(LT .EQ. 0 .OR. LT .EQ. 1) GO TO 20 SPUN 109 CALL ERROR(2) SPUN 110 LT = 0 SPUN 111 20 IF(NPER .GE. 0) GO TO 25 SPUN 112 CALL ERROR(3) SPUN 113 NPER = 0 SPUN 114 25 IF(NBASE .GE. 0) GO TO 30 SPUN 115 CALL ERROR(4) SPUN 116 NBASE = 0 SPUN 117 30 IF(NK .EQ. NDAT + LT + 2 * NPER + NBASE) GO TO 35 SPUN 118 CALL ERROR(5) SPUN 119 NK = NDAT + LT + 2 * NPER + NBASE SPUN 120 35 IF(NK .GT. NKDIM) CALL ERROR(106) SPUN 121 IF(NDAT .GE. 1 .AND. DAT(1) .NE. T(1)) CALL ERROR(107) SPUN 122 C SPUN 123 C SUPPRESS KNOWN CONSTITUENTS SPUN 124 C REPLACE F WITH RESIDUAL TIME SERIES SPUN 125 C COMPUTE QUADRATIC NORM OF F SPUN 126 C CHECK IF RMS VALUE OF RESIDUAL F IS LESS SPUN 127 C THAN EPS * FMAX * ROUND, WHERE SPUN 128 C EPS = EPSARG = SMALLEST NUMBER SO 1 + EPS .GT. 1 SPUN 129 C FMAX = MAXIMUM ABSOLUTE VALUE OF ORIGINAL F SPUN 130 C ROUND ACCOUNTS FOR ACCUMULATED ROUNDOFF IN SPUN 131 C COMPUTING RESIDUAL F SPUN 132 IF(NK .GT. 0) CALL RESID(T, F, NF, SPUN 133 $ NK, DAT, NDAT, LT, PER, NPER, NBASE, SPUN 134 $ A, B, C, NKDIM) SPUN 135 FNORM = 0.0D0 SPUN 136 EPSARG = EPS(ARG) SPUN 137 DO 60 I = 1, NF SPUN 138 60 FNORM = FNORM + F(I) * F(I) SPUN 139 C SPUN 140 C CHECK IF RESIDUAL F CONSISTS OF ROUNDOFF SPUN 141 ICRIT = 1 SPUN 142 IF(DSQRT(FNORM/DFLOAT(NF)) .LT. SPUN 143 $ EPSARG*FMAX*ROUND) ICRIT = 0 SPUN 144 C SPUN 145 C FOR EACH SPECTRAL PERIOD P(I),COMPUTE SPECTRAL VALUE S(I) SPUN 146 C COMPUTE SCALAR PRODUCTS FCOS,FSIN,CC,CS,SS,U,V SPUN 147 C COMPUTE BILINEAR FORMS UAU,UAV,VAV SPUN 148 C COMPUTE PERCENTAGE VARIANCE S SPUN 149 65 DO 130 I = 1, NW SPUN 150 OMEGA = 2.0D0 * PI / P(I) SPUN 151 FCOS = 0.0D0 SPUN 152 FSIN = 0.0D0 SPUN 153 CC = 0.0D0 SPUN 154 CS = 0.0D0 SPUN 155 SS = 0.0D0 SPUN 156 IF(NK .EQ. 0) GO TO 75 SPUN 157 DO 70 J = 1, NK SPUN 158 U(J) = 0.0D0 SPUN 159 70 V(J) = 0.0D0 SPUN 160 75 DO 85 J = 1, NF SPUN 161 WT = OMEGA * T(J) SPUN 162 COSWT = DCOS(WT) SPUN 163 SINWT = DSIN(WT) SPUN 164 FCOS = FCOS + F(J) * COSWT SPUN 165 FSIN = FSIN + F(J) * SINWT SPUN 166 CC = CC + COSWT * COSWT SPUN 167 CS = CS + COSWT * SINWT SPUN 168 SS = SS + SINWT * SINWT SPUN 169 IF(NK .EQ. 0) GO TO 85 SPUN 170 DO 80 K = 1, NK SPUN 171 FUNC = BASE(K, T(J), DAT, NDAT, LT, PER, NPER) SPUN 172 U(K) = U(K) + FUNC * COSWT SPUN 173 80 V(K) = V(K) + FUNC * SINWT SPUN 174 85 CONTINUE SPUN 175 UAU = 0.0D0 SPUN 176 UAV = 0.0D0 SPUN 177 VAV = 0.0D0 SPUN 178 IF(NK .EQ. 0) GO TO 125 SPUN 179 DO 120 J = 1, NK SPUN 180 DO 120 K = 1, NK SPUN 181 UAU = UAU + U(J) * A(J,K) * U(K) SPUN 182 UAV = UAV + U(J) * A(J,K) * V(K) SPUN 183 120 VAV = VAV + V(J) * A(J,K) * V(K) SPUN 184 125 S(I) = 0.0D0 SPUN 185 DET = (CC-UAU) * (SS-VAV) - (CS-UAV) * (CS-UAV) SPUN 186 IF(DABS(DET) .LT. EPSARG) GO TO 130 SPUN 187 S(I) = 100.0D0 * ((SS - VAV) * FCOS * FCOS - SPUN 188 $ 2.0D0 * (CS - UAV) * FCOS * FSIN + SPUN 189 $ (CC - UAU) * FSIN * FSIN) / SPUN 190 $ (DET * FNORM) SPUN 191 130 CONTINUE SPUN 192 RETURN SPUN 193 END SPUN 194