C C C MAIN PROGRAM SPANEQ C C IMPLICIT REAL*8(A-H,O-Z) C C THE FOLLOWING MUST BE DIMENSIONED AS INDICATED BY THE COMMENT C PRECEDING THE DIMENSION STATEMENT. NOTE THAT THE VARIABLES BEGINNING C WITH NDIM (ASSIGNED BELOW THE FOLLOWING DIMENSION STATEMENTS) MUST C AGREE WITH THE ACTUAL DIMENSIONED SIZES. C C NDIMTS BY NDIMNF : REAL*4 SF(12000) DIMENSION F(3,12000),T(3,12000),FF(3,12000) C NDIMTS : DIMENSION NF(3),FNFUNT(3),NDATMS(3) C NDIMTS BY NDIMDM : DIMENSION DATUMS(3,50) C NDIMDM BY 4 DIMENSION KIVL(50,4) C NDIMIV AND 2 BY NDIMIV : DIMENSION TA(55),TB(55),IVL(55),XN(55),XL(55),IDCODE(55) DIMENSION IVLIJ(2,55),DATS(2,55) C NDIMIV BY NDIMFF : DIMENSION SNP(55,70),CNP(55,70),SLP(55,70),CLP(55,70) C NDIMIV BY 2 DIMENSION IDAT(55,2) C NDIMFF : DIMENSION FOFREQ(70),FORPER(70),SP(70),CP(70),SPMQ(70),SPPQ(70) DIMENSION ACOEF(70),BCOEF(70) C NDIMFF BY 7 : DIMENSION FPFAZ(70,7) C NDIMSP : DIMENSION SPFREQ( 500),S( 500) C NDIMUK AND NDIMUK BY NDIMUK : DIMENSION RU(200),D(200),X(200),ITERM(200),RN(200,200) C NDIMBN : DIMENSION RPL(10),RPS(10),NSPVBN(10) NDIMTS=3 NDIMNF=12000 NDIMDM=50 NDIMIV=55 NDIMFF=70 NDIMSP=500 NDIMUK=200 NDIMBN=10 C READ INPUT DATA (INCLUDING TIME SERIES) C IJOB=0 1 CALL READ(T,F,NF,DATUMS,NDATMS,LNTRND,FOFREQ,FORPER,NFFREQ,SPFREQ, +NDIMTS,NDIMNF,NDIMDM,NDIMFF,NDIMSP,ITERM,NSYSNF,NSYSUD,NSYS,NDIMUK +,NPBAND,NDIMBN,RPL,RPS,NSPVBN,TUNITS,FUNITS,FPFAZ,NONSIM,NPRTS,NPR +RS,RHOMAX,NCOS,IDAT,ACOEF,BCOEF,KIVL,FNFUNT,DATS,IDCODE,NDIMIV, + IJOB,NPLTS,NPLRS,NPLSP,IHSCAL,MODE0,SPMAX) C C C ECHO INPUT DATA (INCLUDING TIME SERIES) C C CALL ECHO(T,F,NF,DATUMS,NDATMS,NDIMTS,NDIMNF,NDIMDM,LNTRND,NFFREQ, +FORPER,NSYSNF,NSYSUD,NSYS,NDIMFF,TUNITS,FUNITS,NPRTS,NCOS,FPFAZ,FN +FUNT) TEND=0.D0 NFN= 1 + NSYSNF C C SET UP CONTROL VECTORS FOR HANDLING DATUM BIASES C CALL DATIVL(T,DATUMS,1,NF,NDIMTS,NDIMNF,NDIMDM,TA,TB,IVL,NDIMIV,ND +ATMS,NIVL,IVLIJ,KIVL,XN,XL) IF(NPLTS.NE.0) CALL PLOTT(NFN,F, FUNITS,IHSCAL,1,S,TEND,0, @ IVLIJ,NIVL,T,0.D0,NDIMNF,NDIMTS,SF,NDIMIV,FF,MODE0,SPMAX) C C SET UP CONTROL VECTOR FOR HANDLING ALL SYSTEMATIC NOISE C CALL MITERM(ITERM,NDIMUK,NDATMS,NDIMTS,LNTRND,NFFREQ,NSYSNF,NSYSUD +,NCOS) C C COMPUTE TRIG FUNCTIONS FOR FORCED PERIODS C IF(NFFREQ.NE.0) CALL TRIGPK(NFFREQ,FOFREQ,NIVL,SNP,CNP,SLP,CLP,NDI +MFF,NDIMIV,SP,CP,XN,XL) C C COMPUTE MAGNITUDES OF SYSTEMATIC NOISE AND SOME STATISTICS AND PRINT C THEM C CALL SYSTAT(NF,NDIMTS,NSYS,RN,RU,ITERM,NDIMUK,NIVL,NDIMIV,TA,TB,FO +FREQ,NDIMFF,NDIMNF,T,F,NFFREQ,KIVL,NDIMDM,FPFAZ,D,X,SVAR,VFAC,NDAT +MS,LNTRND,FORPER,NSYSNF,NSYSUD,FUNITS,TUNITS,RHOMAX,NCOS,DATUMS,NP +RRS,RANSP,FNFUNT) IF(NPLRS.NE.0) CALL PLOTT(NFN,F, FUNITS,IHSCAL,1,S,TEND,0, @ IVLIJ,NIVL,T,0.D0,NDIMNF,NDIMTS,SF,NDIMIV,FF,MODE0,SPMAX) C C C COMPUTE SPECTRA C C IF(NPBAND.LE.0) GO TO 85 I1=NSYS+1 I2=NSYS+2 C DO 90 K=1,NPBAND NSPFRQ=NSPVBN(K) C CALL SPKVAL(RPL(K),RPS(K),NSPFRQ,SPFREQ,NDIMSP) C CALL SPCOMP(NSPFRQ,SPFREQ,NONSIM,NF,NSYS,RN,RU,I1,I2,T,F,NCOS,NDAT +MS,LNTRND,NFFREQ,FOFREQ,FPFAZ,KIVL,NDIMDM,NDIMTS,NDIMNF,SP,CP,XN,X +L,CNP,SNP,CLP,SLP,SPMQ,SPPQ,NIVL,IVL,SVAR,NDIMSP,NDIMUK,NDIMFF,NDI +MIV,S,NSYSNF,NSYSUD) C C PRINT SPECTRA C DO 50 I=1,NSPFRQ 50 SPFREQ(I)=3.141592653589793D0*2.D0/SPFREQ(I) C CALL SPLOT(SPFREQ,S,NSPFRQ,K,NDIMSP,TUNITS,RANSP) IF(NPLSP.EQ.0) GO TO 90 TEND=100.D0/RPS(K) TST=100.D0 / RPL(K) DO 60 I=1,NSPFRQ IF(S(I).LT.0.D0) S(I)=0.D0 60 CONTINUE CALL PLOTT(1,F, FUNITS,IHSCAL,2,S,TEND,NSPFRQ, @ IVLIJ,NIVL,T,TST,NDIMNF,NDIMTS,SF,NDIMIV,FF,MODE0,SPMAX) 90 CONTINUE C 85 CONTINUE GOTO1 END SUBROUTINE CCSUM(N1,N2,OMEGA1,OMEGA2,SUM) C C THIS ROUTINE COMPUTES THE SUM FOR K=N1 TO N2 OF C COS(OMEGA1*K)*COS(OMEGA2*K) FOR OMEGA1.NE.2*K*PI.AND.OMEGA2. C NE.2*K*PI C C ROUTINE COSSUM IS CALLED C IMPLICIT REAL*8 (A-H,O-Z) OMEGAM=OMEGA1-OMEGA2 OMEGAP=OMEGA1+OMEGA2 CALL COSSUM(N1,N2,OMEGAM,SUM1) CALL COSSUM(N1,N2,OMEGAP,SUM2) SUM=(SUM1+SUM2)/2.D0 RETURN END SUBROUTINE COSSUM(N1,N2,OMEGA,SUM) C C THIS SUBROUTINE COMPUTES THE SUM FOR K= N1 TO N2 OF COS(K*OMEGA) C FOR OMEGA.NE.2*K*PI C IMPLICIT REAL*8 (A-H,O-Z) N=N2-N1+1 R1=DFLOAT(N1) RN=DFLOAT(N) IF(OMEGA.NE.0.D0) GO TO 1 SUM=RN RETURN 1 ALPHA=(R1-1.D0)*OMEGA SNX2=DSIN(RN*OMEGA/2.D0) SX2=DSIN(OMEGA/2.D0) CAN1X2=DCOS(ALPHA+(RN+1.D0)*OMEGA/2.D0) SUM=SNX2*CAN1X2/SX2 RETURN END SUBROUTINE DATIVL(T,DATUMS,ISER,NF,NDIMTS,NDIMNF,NDIMDM,TA,TB,IVL, +NDIMIV,NDATMS,NIVL,IVLIJ,KIVL,XN,XL) C C C THIS ROUTINE SETS UP CONTROL VECTORS FOR HANDLING DATUM BIASES C C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION T(NDIMTS,NDIMNF),DATUMS(NDIMTS,NDIMDM),NDATMS(NDIMTS),NF +(NDIMTS),TA(NDIMIV),TB(NDIMIV),IVL(NDIMIV),IVLIJ(2,NDIMIV),KIVL(ND +IMDM,4),XN(NDIMIV),XL(NDIMIV) NDAT=NDATMS(ISER) IDAT=2 NIVL=0 TTA=T(ISER,1) IVLIJ(1,1)=1 N=NF(ISER) DO 45 I=2,N NEWIVL=0 IF(I.EQ.N.OR.(T(ISER,I)-T(ISER,I-1)).GT.1.5D0) NEWIVL=1 IF(NDAT.EQ.0.OR.IDAT.GT.NDAT) GO TO 40 IF(DATUMS(ISER,IDAT).GT.T(ISER,I)+0.5D0) GO TO 40 IDAT=IDAT+1 NEWIVL=2 40 IF(NEWIVL.EQ.0) GO TO 45 TTB=T(ISER,I-1) IF(I.EQ.N) TTB=T(ISER,N) IDIVL=IDAT-NEWIVL KIVL(IDIVL,2)=I-1 IF(I.EQ.N) KIVL(IDIVL,2)=N NIVL=NIVL+1 IVLIJ(2,NIVL)=I-1 IF(I.EQ.N) IVLIJ(2,NIVL)=N KIVL(IDIVL,4)=NIVL TA(NIVL)=TTA TB(NIVL)=TTB XN(NIVL)=TTB-TTA+1.D0 XL(NIVL)=TTB+TTA IVL(NIVL)=IDIVL TTA=T(ISER,I) IF(I.EQ.N) GO TO 45 IVLIJ(1,NIVL+1)=I 45 CONTINUE KIVL(1,1)=1 KIVL(1,3)=1 IF(NDAT.LT.2) RETURN DO 50 I=2,NDAT KIVL(I,1)=KIVL(I-1,2)+1 50 KIVL(I,3)=KIVL(I-1,4)+1 RETURN END FUNCTION DESIGN(I,J,IJ,KIVL,NDIMDM,T,F,NDIMTS,NDIMNF,OMEGA,ISER,IF +UNC,FAZ) C C THIS ROUTINE COMPUTES ELEMENTS OF THE DESIGN MATRIX FOR USE IN C ROUTINES NORMEQ, SYSTAT, RNADD AND SPCOMP. C IMPLICIT REAL*8(A-H,O-Z) DIMENSION KIVL(NDIMDM,4),T(NDIMTS,NDIMNF),F(NDIMTS,NDIMNF) GO TO(25,35,45,55,65,75,85),IJ C C DATUM BIASES C 25 DESIGN=1.D0 IF(I.LT.KIVL(J,1).OR.I.GT.KIVL(J,2)) DESIGN=0.D0 RETURN C C LINEAR TREND C 35 DESIGN=T(1,I) RETURN C C COSINE TERM C 45 DESIGN=DCOS(OMEGA*T(1,I)) RETURN C C SINE TERM C 55 DESIGN=DSIN(OMEGA*T(1,I)) RETURN C C COS TERM (FIXED PHASE) C 65 DESIGN=DCOS(OMEGA*T(1,I)-FAZ) RETURN C C NUMERICAL FUNCTION C 75 DESIGN=F(ISER,I) RETURN C C USER DEFINED FUNCTION C 85 DESIGN=FUNCTN(T(1,I),IFUNC) RETURN END SUBROUTINE ECHO(T,F,NF,DATUMS,NDATMS,NDIMTS,NDIMNF,NDIMDM,LNTRND,N +FFREQ,FORPER,NSYSNF,NSYSUD,NSYS,NDIMFF,TUNITS,FUNITS,NPRTS,NCOS,FP +FAZ,FNFUNT) C C THIS ROUTINE ECHOS THE INPUT DATA (INCLUDING ANY TIME SERIES BUT NOT C INCLUDING INPUT PARAMETERS OF GENERATED TIME SERIES) C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION T(NDIMTS,NDIMNF),F(NDIMTS,NDIMNF),NF(NDIMTS) DIMENSION DATUMS(NDIMTS,NDIMDM),NDATMS(NDIMTS),FORPER(NDIMFF) DIMENSION FPFAZ(NDIMFF,7),FNFUNT(NDIMTS) ISER=1 NDAT=NDATMS(ISER) PRINT 1001 IF(NSYS.EQ.0) GO TO 35 PRINT 1002 J=1 IF(NDAT.EQ.0) GO TO 5 IF(J.EQ.NDAT) PRINT 1003,J IF(J.NE.NDAT) PRINT 1004,J,NDAT 5 J=NDAT+1 IF(LNTRND.EQ.0) GO TO 10 PRINT 1005,J J=J+1 10 IF(NFFREQ.EQ.0) GO TO 20 DO 15 I=1,NFFREQ K=J+1 PRINT 1006,J,K,FORPER(I) 15 J=J+2 20 IF(NCOS.EQ.0) GO TO 25 DO 21 I=1,NCOS FAZ=FPFAZ(NFFREQ+I,2)*180.D0/3.14159265359D0 PRINT 1007,J,FORPER(NFFREQ+I),FAZ J=J+1 21 CONTINUE 25 IF(NSYSNF.EQ.0) GO TO 30 K=NSYSNF+J-1 IF(K.EQ.J) PRINT 1008,J IF(K.NE.J) PRINT 1009,J,K J=K+1 30 IF(NSYSUD.EQ.0) GO TO 35 K=NSYSUD+J-1 IF(K.EQ.J) PRINT 1010,J IF(K.NE.J) PRINT 1011,J,K J=K 35 PRINT 1012 IF(NPRTS.GT.0) GO TO 45 PRINT 1001 PRINT 1013,TUNITS,FUNITS CALL FPLOT(ISER,T,F,NF,DATUMS,NDATMS,NDIMTS,NDIMNF,NDIMDM) PRINT 1001 IF(NSYSNF.EQ.0) GO TO 45 DO 40 I=1,NSYSNF J=I+1 PRINT 1014,I,FNFUNT(J) CALL FPLOT(J,T,F,NF,DATUMS,NDATMS,NDIMTS,NDIMNF,NDIMDM) PRINT 1001 40 CONTINUE 1001 FORMAT(' ',110('-'),//) 1002 FORMAT(' ',10X,'FORCED CONSTITUENTS (SYSTEMATIC NOISE):',/) 1003 FORMAT(' ',10X,I2,6X,'DATUM BIAS') 1004 FORMAT(' ',10X,I2,' -',I3,1X,'DATUM BIASES') 1005 FORMAT(' ',10X,I2,6X,'LINEAR TREND') 1006 FORMAT(' ',10X,I2,' -',I3,1X,'COS AND SIN OF FORCED PERIOD=',F12.6 +) 1007 FORMAT(' ',10X,I2,6X,'COSINE TERM WITH PERIOD=',F12.6,' PHASE LAG= +',F8.2,' DEGREES') 1008 FORMAT(' ',10X,I2,6X,'NUMERICAL FUNCTION') 1009 FORMAT(' ',10X,I2,' -',I3,1X,'NUMERICAL FUNCTIONS') 1010 FORMAT(' ',10X,I2,6X,'USER DEFINED FUNCTION') 1011 FORMAT(' ',10X,I2,' -',I3,1X,'USER DEFINED FUNCTION') 1012 FORMAT(' ',//) 1013 FORMAT(' ',10X,'INPUT TIME SERIES: ','(TIME UNITS=',A8,',TIME SERI +ES UNITS=',A8,')',/) 1014 FORMAT(' ',//,' ',10X,'NUMERICAL FUNCTION NUMBER',I2,'(UNITS= ',A8 +,') :',/) 45 RETURN END SUBROUTINE FPLOT(ISER,T,F,NF,DATUMS,NDATMS,NDIMTS,NDIMNF,NDIMDM) C C THIS ROUTINE PLOTS A TIME SERIES ON THE LINE PRINTER. C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION T(NDIMTS,NDIMNF),F(NDIMTS,NDIMNF),NF(NDIMTS),DATUMS(NDIM +TS,NDIMDM),NDATMS(NDIMTS),IPLOT(80) DATA IBLNK/' '/,ISTAR/'*'/ PRINT 1001 DO 5 I=1,80 5 IPLOT(I)=IBLNK FMIN=F(ISER,1) FMAX=FMIN N=NF(ISER) DO 10 I=1,N IF(F(ISER,I).LT.FMIN) FMIN=F(ISER,I) IF(F(ISER,I).GT.FMAX) FMAX=F(ISER,I) 10 CONTINUE ND=NDATMS(ISER) IDAT=1 INTVAL=0 TIME=T(ISER,1) DO 20 I=1,N IF(ISER.GT.1) GO TO 11 IF(ND.EQ.0.OR.IDAT.GT.ND) GO TO 11 IF(DATUMS(ISER,IDAT).GT.T(ISER,I)) GO TO 11 PRINT 1002, IDAT IDAT=IDAT+1 INTVAL=INTVAL+1 GO TO 15 11 IF(T(ISER,I)-TIME.LT.1.5D0.AND.INTVAL.GT.0) GO TO 15 INTVAL=INTVAL+1 PRINT 1004,INTVAL 15 KF=1+IDINT(79.D0*(F(ISER,I)-FMIN)/(FMAX-FMIN)) IF(KF.LT.1) KF=1 IF(KF.GT.80) KF=80 IPLOT(KF)=ISTAR TIME=T(ISER,I) PRINT 1003,I,T(ISER,I),F(ISER,I),(IPLOT(J),J=1,80) 20 IPLOT(KF)=IBLNK 1001 FORMAT(' ',5X,'I',8X,'T(I)',8X,'F(I)',/) 1002 FORMAT(' ',30X,28('-'),'BEGINNING OF DATUM',I5,29('-')) 1003 FORMAT(' ',2X,I4,2D12.4,80A1) 1004 FORMAT(' ',30X,27('-'),'BEGINNING OF INTERVAL',I5,27('-')) RETURN END FUNCTION FUNCTN(TIME,KUSERD) C C THIS ROUTINE COMPUTES VALUES OF ANY USER DEFINED FUNCTIONS TO BE C FORCED AS SYSTEMATIC NOISE C IMPLICIT REAL*8 (A-H,O-Z) GO TO(1,2,3),KUSERD 1 FUNCTN=3*DEXP(-TIME/25.D0) RETURN 2 FUNCTN=DEXP(-TIME/12.D0) RETURN 3 FUNCTN=1.D0 RETURN END SUBROUTINE MAGFAZ(A,B,SAA,SAB,SBB,VFAC,AMP,PHAZ,SAMP,SPHAZ) C C THIS ROUTINE COMPUTES THE MAGNITUDES AND PHASE LAG (AND THEIR C STANDARD DEVIATIONS) OF A TRIGONOMETRIC TERM. C IMPLICIT REAL*8 (A-H,O-Z) PI=3.141592653589793D0 AA=A*A AB=A*B BB=B*B AMP=DSQRT(AA+BB) SAMP=DSQRT(VFAC*(AA*SAA+2.D0*AB*SAB+BB*SBB))/AMP PHAZ=2.D0*DATAN(B/(A+AMP)) IF(PHAZ.LT.0.D0) PHAZ=PHAZ+2.D0*PI PHAZ=PHAZ*180.D0/PI SPHAZ=DSQRT(VFAC*(BB*SAA-2.D0*AB*SAB+AA*SBB))/AMP/AMP*180.D0/PI RETURN END SUBROUTINE MITERM(ITERM,NDIMUK,NDATMS,NDIMTS,LNTRND,NFFREQ,NSYSNF, +NSYSUD,NCOS) C C THIS ROUTINE SETS UP A CONTROL VECTOR FOR HANDLING SYSTEMATIC NOISE C COMPUTATIONS C DIMENSION ITERM(NDIMUK),NDATMS(NDIMTS) NDAT=NDATMS(1) C C DATUM BIASES C IF(NDAT.EQ.0) GO TO 10 DO 5 J=1,NDAT 5 ITERM(J)=1 10 NCOLMN=NDAT+1 C C LINEAR TREND C IF(LNTRND.NE.1) GO TO 15 ITERM(NCOLMN)=2 NCOLMN=NCOLMN+1 C C FORCED FREQUENCIES C 15 IF(NFFREQ.EQ.0) GO TO 35 DO 20 J=1,NFFREQ ITERM(NCOLMN)=3 ITERM(NCOLMN+1)=4 20 NCOLMN=NCOLMN+2 C CFORCED FREQUENCIES WITH FORCED PHASE LAG C 35 IF(NCOS.EQ.0) GO TO 38 DO 36 J=1,NCOS ITERM(NCOLMN)=5 36 NCOLMN=NCOLMN+1 C C FORCED NUMERICAL FUNCTIONS 38 IF(NSYSNF.EQ.0) GO TO 45 DO 40 J=1,NSYSNF ITERM(NCOLMN)=6 40 NCOLMN=NCOLMN+1 C C FORCED USER DEFINED FUNCTIONS C 45 IF(NSYSUD.EQ.0) GO TO 55 DO 50 J=1,NSYSUD ITERM(NCOLMN)=7 50 NCOLMN=NCOLMN+1 55 RETURN END SUBROUTINE NORMEQ(I1,I2,RN,NDIMUK,RU,ITERM,NIVL,NDIMIV,TA,TB,FOFRE +Q,NDIMFF,NDIMNF,NF,T,NDIMTS,F,NFFREQ,KIVL,NDIMDM,FPFAZ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION RN(NDIMUK,NDIMUK),RU(NDIMUK),ICNT(2,7),ITERM(NDIMUK) DIMENSION TA(NDIMIV),TB(NDIMIV),FOFREQ(NDIMFF),NF(NDIMTS) DIMENSION T(NDIMTS,NDIMNF),F(NDIMTS,NDIMNF),KIVL(NDIMDM,4) DIMENSION FPFAZ(NDIMFF,7) N=NF(1) C C ZERO PORTIONS OF RN AND RU C DO 5 J=I1,I2 RU(J)=0.D0 DO 5 I=1,J 5 RN(I,J)=0.D0 C C INITIALIZE ICNT C DO 10 J=1,7 10 ICNT(2,J)=0 ICNT(2,6)=1 ICNT(2,5)=NFFREQ C C FILL RN C DO 20 J=1,I2 ICNT(2,ITERM(J))=ICNT(2,ITERM(J))+1 DO 15 K=1,7 15 ICNT(1,K)=0 ICNT(1,6)=1 ICNT(1,5)=NFFREQ DO 20 I=1,J ICNT(1,ITERM(I))=ICNT(1,ITERM(I))+1 CALL RNADD(I,J,RN,NDIMUK,ICNT,ITERM,NIVL,NDIMIV,TA,TB,FOFREQ,NDIMF +F,NDIMNF,N,KIVL,NDIMDM,T,F,NDIMTS,FPFAZ) 20 CONTINUE C C FILL RU C DO 25 K=1,7 25 ICNT(1,K)=0 ICNT(1,6)=1 ICNT(1,5)=NFFREQ DO 35 I=I1,I2 II=ITERM(I) ICNT(1,II)=ICNT(1,II)+1 IF(II.EQ.3.OR.II.EQ.4) OMEGA=FOFREQ(ICNT(1,3)) IF(II.NE.5) GO TO 30 FAZ=FPFAZ(ICNT(1,5),2) OMEGA=FOFREQ(ICNT(1,5)) 30 ISER=ICNT(1,6) IFUNC=ICNT(1,7) DO 35 J=1,N RU(I)=RU(I)+F(1,J)*DESIGN(J,I,II,KIVL,NDIMDM,T,F,NDIMTS,NDIMNF,OME +GA,ISER,IFUNC,FAZ) 35 CONTINUE RETURN END SUBROUTINE PAMPR(NDATMS,NDIMTS,LNTRND,NFFREQ,FORPER,NDIMFF,NSYSNF, +NSYSUD,X,NDIMUK,FPFAZ,FUNITS,VFAC,RN,TUNITS,RHOMAX,NCOS,NSYS,FNFUN +T) C C THIS ROUTINE PRINTS THE MAGNITUDES (AND THEIR STANDARD DEVIATIONS) C OF THE SYSTEMATIC NOISE COMPONENTS WHICH ARE COMPUTED IN ROUTINE C SYSTAT. C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION NDATMS(NDIMTS),FORPER(NDIMFF),X(NDIMUK),FPFAZ(NDIMFF,7) DIMENSION RN(NDIMUK,NDIMUK),FNFUNT(NDIMTS) PI=3.141592653589793D0 PRINT 1000 NDAT=NDATMS(1) C C DATUM BIASES C IF(NDAT.EQ.0) GO TO 10 DO 5 I=1,NDAT STD=DSQRT(RN(I,I)*VFAC) PRINT 1001,I,X(I),STD,FUNITS 5 CONTINUE 10 K=NDAT+1 C C LINEAR TREND C IF(LNTRND.EQ.0) GO TO 15 STD=DSQRT(RN(K,K)*VFAC) PRINT 1002,X(K),STD,FUNITS,TUNITS K=K+1 C C FORCED FREQUENCIES C 15 IF(NFFREQ.EQ.0) GO TO 25 DO 20 I=1,NFFREQ CALL MAGFAZ(X(K),X(K+1),RN(K,K),RN(K,K+1),RN(K+1,K+1),VFAC,AMP,PHA +Z,SAMP,SPHAZ) IF(FPFAZ(I,1).LT.0.5D0) GO TO 16 RAMP=AMP/FPFAZ(I,2) DFAZ=PHAZ-FPFAZ(I,4) IF(FPFAZ(I,6).LT.0.5D0) GO TO 18 DFAZ=DFAZ-180.D0 18 IF(DFAZ.LT.0.D0) DFAZ=DFAZ+360.D0 SRAMP=DSQRT(AMP*AMP/FPFAZ(I,2)**4*FPFAZ(I,3)**2+SAMP*SAMP/FPFAZ(I, +2)**2) SDFAZ=DSQRT(SPHAZ*SPHAZ+FPFAZ(I,5)**2) 16 PRINT 1003,FORPER(I),TUNITS STD=DSQRT(VFAC*RN(K,K)) PRINT 1004,X(K),STD,FUNITS STD=DSQRT(VFAC*RN(K+1,K+1)) PRINT 1005,X(K+1),STD,FUNITS PRINT 1006,AMP,SAMP,FUNITS PRINT 1007,PHAZ,SPHAZ IF(FPFAZ(I,1).LT.0.5D0) GO TO 19 PRINT 1008,FPFAZ(I,2),FPFAZ(I,3),FPFAZ(I,7) PRINT 1009,FPFAZ(I,4),FPFAZ(I,5) PRINT 1010,RAMP,SRAMP,FUNITS,FPFAZ(I,7) PRINT 1011,DFAZ,SDFAZ 19 PRINT 1012 20 K=K+2 C C FORCED FREQUENCIES WITH FORCED PHASE C 25 IF(NCOS.EQ.0) GO TO 28 DO 27 I=1,NCOS FPER=FORPER(NFFREQ+I) COEF=X(K) FAZ=FPFAZ(NFFREQ+I,2)*180.D0/PI STD=DSQRT(RN(K,K)*VFAC) PRINT 1017,FPER,TUNITS,FAZ PRINT 1016,COEF,STD,FUNITS 27 K=K+1 28 CONTINUE C C FORCED NUMERICAL FUNCTIONS C IF(NSYSNF.EQ.0) GO TO 35 DO 40 I=1,NSYSNF STD=DSQRT(RN(K,K)*VFAC) PRINT 1018,I,X(K),STD,FUNITS,FNFUNT(I+1) 40 K=K+1 35 CONTINUE C C FORCED USER DEFINED FUNCTIONS C IF(NSYSUD.EQ.0) GO TO 45 DO 50 I=1,NSYSUD STD=DSQRT(RN(K,K)*VFAC) PRINT 1019,I,X(K),STD,FUNITS 50 K=K+1 45 CONTINUE PRINT 1013 IF(NSYS.LT.2) GO TO 30 IF(RHOMAX.GT.0.999D0) GO TO 30 PRINT 1015,RHOMAX IEND=NSYS-1 DO 26 I=1,IEND JSTART=I+1 DO 26 J=JSTART,NSYS RHO=RN(I,J)/DSQRT(RN(I,I)*RN(J,J)) IF(DABS(RHO).LT.RHOMAX) GO TO 26 PRINT 1014,RHO,I,J 26 CONTINUE PRINT 1012 PRINT 1013 1000 FORMAT(' ',//,' ',110('-'),/,' ',10X,'PRELIMINARY MAGNITUDES OF SY +STEMATIC NOISE COMPONENTS:',//) 1001 FORMAT(' ',10X,'DATUM BIAS #',I5,' =',D14.6,6X,'+/-',D13.6,1X,A8,/ +) 1002 FORMAT(' ',10X,'LINEAR TREND',6X,'=',D14.6,6X,'+/-',D13.6,1X,A8,'/ +',A8,/) 1003 FORMAT(' ',10X,'FORCED PERIOD',5X,'=',D14.6,1X,A8,':') 1004 FORMAT(' ',16X,'COSINE COEFF= ',D13.6,6X,'+/-',D13.6,1X,A8) 1005 FORMAT(' ',16X,'SINE COEFF = ',D13.6,6X,'+/-',D13.6,1X,A8) 1006 FORMAT(' ',16X,'AMPLITUDE = ',D13.6,6X,'+/-',D13.6,1X,A8) 1007 FORMAT(' ',16X,'PHASE',7X,'= ',D13.6,6X,'+/-',D13.6,' DEGREES') 1008 FORMAT(' ',16X,'FORCING AMPL= ',D13.6,6X,'+/-',D13.6,1X,A8) 1009 FORMAT(' ',16X,'FORCNG PHASE= ',D13.6,6X,'+/-',D13.6,' DEGREES') 1010 FORMAT(' ',16X,'AMPL RATIO= ',D13.6,6X,'+/-',D13.6,1X,A8,'/',A8) 1011 FORMAT(' ',16X,'PHASE DIFF = ',D13.6,6X,'+/-',D13.6,' DEGREES') 1012 FORMAT(' ',/) 1013 FORMAT(' ',110('-'),//) 1014 FORMAT(' ',16X,'CORRELATION COEFFICIENT=',F8.4,' FOR PARAMETERS',I +4,' AND',I4) 1015 FORMAT(' ',16X,'CORRELATION COEFFICIENTS GREATER IN MAGNITUDE THAN +',F7.4,' :',///) 1016 FORMAT(' ',16X,'COSINE COEFF= ',D13.6,6X,'+/-',D13.6,1X,A8) 1017 FORMAT(' ',10X,'FORCED PERIOD',5X,'=',D14.6,1X,A8,5X,'FORCED PHASE + LAG',2X,'=',D14.6,1X,'DEGREES :') 1018 FORMAT(' ',10X,'FORCED NUMERICAL FUNCTION #',I3,':',/,' ',16X,'COE +FFICIENT = ',D13.6,6X,'+/-',D13.6,1X,A8,'/',A8,/) 1019 FORMAT(' ',10X,'FORCED USER DEFINED FUNCTION #',I3,':',/,' ',16X,' +COEFFICIENT = ',D13.6,6X,'+/-',D13.6,1X,A8,/) 30 RETURN END SUBROUTINE PARGEN(NGEN,LT,COEFLT,NBASE,IDAT,DATS,FOFREQ,ACOEF,BCOE +F,NPER,NDIMDM,NDIMTS,NDIMFF,PI2,RNOISE,ISER,IDCODE,NDIMIV) C C THIS ROUTINE PRINTS INPUT PARAMETERS OF ANY GENERATED TIME SERIES. C IMPLICIT REAL*8(A-H,O-Z) DIMENSION IDAT(NDIMIV,2),DATS(NDIMTS,NDIMIV),FOFREQ(NDIMFF) DIMENSION ACOEF(NDIMFF),BCOEF(NDIMFF),IDCODE(NDIMIV) PRINT 1007 PRINT 1001,ISER KDAT=0 DO 1 I=1,NGEN IF(ISER.NE.1) GO TO 5 IF(IDCODE(I).NE.0) GO TO 5 KDAT=KDAT+1 PRINT 1009,KDAT 5 PRINT 1002,I,(IDAT(I,J),J=1,2),DATS(ISER,I) 1 CONTINUE IF(LT.EQ.0) GO TO 2 PRINT 1003,COEFLT 2 IF(NPER.EQ.0) GO TO 4 DO 3 I=1,NPER PER=PI2/FOFREQ(I) PRINT 1004,I,PER,ACOEF(I),BCOEF(I) 3 CONTINUE 4 PRINT 1008,RNOISE IF(NBASE.EQ.0) PRINT 1005 IF(NBASE.GT.0) PRINT 1006,NBASE PRINT 1007 1001 FORMAT(' ',//,' ',10X,'PARAMETERS OF GENERATED TIME SERIES # ',I2, +' :',//) 1002 FORMAT(' ',10X,'INTERVAL BIAS ',I3,' FROM T= ',I6,' TO ',I6,' = ', +D13.6,/) 1003 FORMAT(' ',10X,'LINEAR TREND= ',D13.6,'*T',/) 1004 FORMAT(' ',10X,'PERIOD ',I3,' = ',D13.6,' (COSINE COEF. = ',D13.6, +' SINE COEF. = ',D13.6,' )',/) 1005 FORMAT(' ',10X,'NO OTHER CONSTITUENTS',/) 1006 FORMAT(' ',10X,'THERE ARE ',I3,' OTHER CONSTITUENTS (SEE FUNCTION +SUBPROGRAM FUNCTN)',/) 1007 FORMAT(' ',110('-')) 1008 FORMAT(' ',10X,'RANDOM NOISE WITH LARGEST ABSOLUTE VALUE = ',D13.6 +,/) 1009 FORMAT(' ',10X,'DATUM #',I3,' :',/) RETURN END SUBROUTINE PLOTT(NFN,G, TITLE2,HSCALE,MODE,S,T,NS,IVLIJ,NIVL,TT, & T1,NDIMNF,NDIMTS,SF,NDIMIV,F,MODE0,SPMAX) INTEGER HSCALE REAL*8 TITLE2,G,S,T,F,TT,T1 REAL*8 UNIT1(6),UNIT2(6) DIMENSION G(NDIMTS,1),TITLE2(NDIMTS),SF(NDIMNF),S(1),NF(3), @ IVLIJ(2,NDIMIV),F(NDIMTS,NDIMNF),TT(NDIMTS,1) DATA UNIT1/8H MSEC , 8H MBAR , & 8H DEG C , 8H % VAR , & 8H CM , 8H INCHES / DATA UNIT2/8H MSEC \, 8H MBAR \, & 8H DEG C \, 8H % VAR \, & 8H CM \, 8H INCHES\/ IF(MODE.EQ.2) GOTO 35 DO 31 I=1,NFN DO 31 J=1,NDIMNF 31 F(I,J)=0.D0 DO 33 I=1,NIVL N1=IVLIJ(1,I) N2=IVLIJ(2,I) DO 32 J=N1,N2 JJ=IDINT(TT(1,J)) 32 F(1,JJ)=G(1,J) 33 CONTINUE DO 34 I=1,NFN DO 34 J=1,NDIMNF IF(F(I,J).NE.0.D0) NF(I)=J 34 CONTINUE 35 NPLOT=0 CALL AREA(6.,8.) CALL SETPLT(0.,0.,6.,8.) CALL SHIFT(0.,8.) SAXIS=0. IF(MODE.EQ.1) GO TO 21 DO 20 I=1,NS SF(I)=S(I) 20 CONTINUE FMAX=SPMAX FMIN=0. STIME=1. ETIME=NS NN=1 NFN=1 HMON=(T-T1)/6. SAXIS=T1 IREP=0 GO TO 22 21 NMAX=1 DO 10 I=1,NFN IF(NF(I).GT.NMAX) NMAX=NF(I) 10 CONTINUE LENGTH=NMAX NN=1 N=1 N1=0 STIME=1. ETIME=HSCALE*6. 1 IF(FLOAT(LENGTH).GT.ETIME) IREP=1 IF(IREP.NE.1) GO TO 8 C FIND MAX AND MIN 8 IF(N.EQ.N1) GO TO 9 FMAX=-999999. FMIN=999999. DO 3 J=1,LENGTH SF(J)=F(N,J) IF(SF(J).LT.FMIN) FMIN=SF(J) IF(SF(J).GT.FMAX) FMAX=SF(J) 3 CONTINUE MAX=(FMAX+10.)/10. FMAX=10*MAX IF(FMAX.LT.10.) FMAX=10. MIN=FMIN/10. FMIN=10*MIN FMIN =FMIN - 10 IF(MODE0.NE.0) FMIN=0. FMAX=FMAX + 10 C PLOT AXES 9 HMON=HSCALE/720. 22 HINC=(FMAX-FMIN)/2. DO 25 I=1,6 IF(TITLE2(N).EQ.UNIT1(I)) KK=I 25 CONTINUE IF(MODE.EQ.2) KK=4 CALL AREA(5.,1.5) CALL SETPLT(0.,0.,5.,1.5) CALL SHIFT(.5,.3) 4 CALL SETPLT(STIME,FMIN,ETIME,FMAX) IF(MODE.EQ.1) CALL AXS2(STIME,FMIN,ETIME,FMIN,'X 30 DAYS\',SAXIS, @ HMON,-6,.12,.1,(.05,.05)) IF(MODE.EQ.2) CALL AXS2(STIME,FMIN,ETIME,FMIN,'CPH X 100\', & SAXIS,HMON,-6,.12,.1,(.05,.05)) CALL AXS2(STIME,FMIN,STIME,FMAX,UNIT2(KK),FMIN,HINC,2,.12,.1, & (.05,.05)) C PLOT FUNCTION IC=0 ISTIME=STIME NGAP=ETIME IF(MODE.EQ.2) GOTO 36 IF(NGAP.GT.NF(N)) NGAP=NF(N) 36 CONTINUE DO 5 K=ISTIME,NGAP X=K IF(MODE.EQ.2) GO TO 26 IF(MODE0.NE.0) GO TO 26 IF(SF(K).EQ.0.) IC=0 IF(SF(K).EQ.0.) GO TO 5 26 CALL NOWPLT(IC,X,SF(K)) IC=1 5 CONTINUE NPLOT=NPLOT+1 IF(NPLOT.LT.3) GO TO 14 CALL ENDPLT NPLOT=0 CALL AREA(6.,8.) CALL SETPLT(0.,0.,6.,8.) CALL SHIFT(0.,8.) 14 IF(IREP.NE.1) N=N+1 IF(IREP.NE.1) GO TO 6 STIME=ETIME+1 ETIME=STIME+HSCALE*6-1 IREP=0 IF(FLOAT(LENGTH).GT.ETIME) IREP=1 NN=NN+1 N1=N SAXIS=SAXIS+HMON*6. IF(NPLOT.EQ.0) GO TO 15 CALL SETPLT(0.,0.,5.,1.5) CALL SHIFT(-.5,-2.5) 15 IF(N.EQ.1) GO TO 1 GO TO 8 6 IF(N.GT.NFN) GO TO 7 IF(NPLOT.EQ.0) GO TO 16 CALL SETPLT(0.,0.,5.,1.5) CALL SHIFT(-.5,-2.5) 16 SAXIS=0. NN=1 N1=N-1 STIME=1. ETIME=HSCALE*6 IF(FLOAT(LENGTH).GT.ETIME) IREP=1 GO TO 8 7 CALL ENDPLT RETURN END SUBROUTINE READ(T,F,NF,DATUMS,NDATMS,LNTRND,FOFREQ,FORPER,NFFREQ,S +PFREQ,NDIMTS,NDIMNF,NDIMDM,NDIMFF,NDIMSP,ITERM,NSYSNF,NSYSUD,NSYS, +NDIMUK,NPBAND,NDIMBN,RPL,RPS,NSPVBN,TUNITS,FUNITS,FPFAZ,NONSIM,NPR +TS,NPRRS,RHOMAX,NCOS,IDAT,ACOEF,BCOEF,KIVL,FNFUNT,DATS,IDCODE,NDIM +IV,IJOB,NPLTS,NPLRS,NPLSP,IHSCAL,MODE0,SPMAX) C C THIS ROUTINE READS THE INPUT DATA. C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION T(NDIMTS,NDIMNF),F(NDIMTS,NDIMNF),NF(NDIMTS) DIMENSION DATUMS(NDIMTS,NDIMDM),NDATMS(NDIMTS),FOFREQ(NDIMFF) DIMENSION FORPER(NDIMFF),SPFREQ(NDIMSP),ITERM(NDIMUK) DIMENSION RPL(NDIMBN),RPS(NDIMBN),IDAT(NDIMIV,2),ACOEF(NDIMFF) DIMENSION NSPVBN(NDIMBN),FPFAZ(NDIMFF,7),BCOEF(NDIMFF),KIVL(NDIMDM +,4) DIMENSION FNFUNT(NDIMTS),TITLE(10),DATS(NDIMTS,NDIMIV) DIMENSION IDCODE(NDIMIV),PEREM(20),AMREM(20),PHREM(20) READ(5,1003,END=99) (TITLE(I), I=1,10) PRINT 1009, (TITLE(I), I=1,10) PI2=3.141592653589793*2.D0 IR=5 IJOB=IJOB+1 IF(IJOB.GT.1) GOTO 45 C C READ CONTROL VARIABLES C READ(IR,1001) LNTRND,NPBAND,NSYSNF,NSYSUD,NONSIM,NPRTS,NPRRS, + MAXRHO,NGENTS,NGENSY,NREMOV,NPLTS,NPLRS,NPLSP,IHSCAL + ,MODE0 RHOMAX=DFLOAT(MAXRHO)/10.D0 IF(MAXRHO.EQ.0) RHOMAX=1.D0 NFFREQ=0 NCOS=0 C C READ PERIODS, AMPLITUDES AND PHASE LAGS OF CONSTITUENTS TO BE REMOVED C FROM MAIN TIME SERIES BEFORE ANALYSIS (IF ANY) C IF(NREMOV.LE.0) GO TO 2 PRINT 1011 DO 1 I=1,NREMOV READ(IR,1010) PEREM(I),AMREM(I),PHREM(I) PRINT 1012, PEREM(I),AMREM(I),PHREM(I) PEREM(I)=PI2/PEREM(I) PHREM(I)=PHREM(I)*PI2/360.D0 1 CONTINUE 2 CONTINUE C C READ FORCED PERIODS,ETC C DO 5 I=1,NDIMFF READ(IR,1002) FORPER(I),(FPFAZ(I,J),J=1,7) IF(FORPER(I).EQ.0.D0) GO TO 10 IF(FPFAZ(I,1).LT.1.5D0) NFFREQ=NFFREQ+1 IF(FPFAZ(I,1).LT.1.5D0) GO TO 5 NCOS=NCOS+1 FPFAZ(I,2)=FPFAZ(I,2)*PI2/360.D0 5 CONTINUE C C READ SPECTRA PARAMETERS (IF ANY) C 10 IF(NPBAND.EQ.0) GO TO 45 DO 75 I=1,NPBAND READ(IR,1005) RPL(I),RPS(I),NSPVBN(I),SPMAX 75 CONTINUE IF(SPMAX.LT.0.01D0) SPMAX=100.D0 C C C 45 ISER=1 C C READ TIME SERIES C IF(IJOB.GT.1) GOTO 46 READ(IR,1004) TUNITS,FUNITS IF(NGENTS.GT.0) GO TO 70 46 CALL READTS(ISER,T,F,NF,FAC,CNST,DATUMS,NDATMS,NDIMTS,NDIMNF,NDIMD +M,IDCODE,NDIMIV,IJOB) IF(IJOB.GT.1) RETURN GO TO 110 C C READ DATA FOR GENERATING TIME SERIES AND GENERATE TIME SERIES C 70 CALL TSGEN(NGENTS,ISER,IDAT,NDIMDM,DATUMS,NDIMTS,NDIMNF,T,F,NDIMFF +,FOFREQ,ACOEF,BCOEF,PI2,DATS,NF,NDATMS,IDCODE,NDIMIV) C C READ TIME SERIES FOR FORCED NUMERICAL FUNCTION C 110 IF(NSYSNF.EQ.0) GO TO 65 DO 50 I=1,NSYSNF J=I+1 READ(IR,1004) FNFUNT(J) IF(NGENSY.GT.0) GO TO 80 CALL READTS(J,T,F,NF,FAC,CNST,DATUMS,NDATMS,NDIMTS,NDIMNF,NDIMDM,I +DCODE,NDIMIV) GO TO 50 80 CALL TSGEN(NGENSY,J,IDAT,NDIMDM,DATUMS,NDIMTS,NDIMNF,T,F,NDIMFF,FO +FREQ,ACOEF,BCOEF,PI2,DATS,NF,NDATMS,IDCODE,NDIMIV) 50 CONTINUE 65 NSYS=LNTRND+NDATMS(1)+2*NFFREQ+NSYSNF+NSYSUD+NCOS IF(NFFREQ.EQ.0) GO TO 115 DO 120 I=1,NFFREQ 120 FOFREQ(I)=PI2/FORPER(I) 115 IF(NCOS.EQ.0) GO TO 125 NCS=NCOS+NFFREQ NST=NFFREQ+1 DO 130 I=NST,NCS 130 FOFREQ(I)=PI2/FORPER(I) 125 CONTINUE IF(NREMOV.LE.0) RETURN N=NF(1) DO 135 I=1,N TI=T(1,I) DO 135 J=1,NREMOV WT=PEREM(J)*TI-PHREM(J) F(1,I)=F(1,I)-AMREM(J)*DCOS(WT) 135 CONTINUE 1001 FORMAT(16I5) 1002 FORMAT(7F10.6,A8) 1003 FORMAT(10A8) 1004 FORMAT(2A8) 1005 FORMAT(2F10.5,I10,F10.5) 1009 FORMAT(' ',10A8) 1010 FORMAT(3F10.6) 1011 FORMAT(' ',10X,'CONSTITUENTS TO BE REMOVED BEFORE ANALYSIS:',/,' ' +,10X,1X,'PERIOD',5X,'AMPLITUDE',3X,'PHASE LAG',/) 1012 FORMAT(' ',7X,F12.6,F10.3,F12.2,/) RETURN 99 STOP END SUBROUTINE READTS(ISER,T,F,NF,FACTOR,CONST,DATUMS,NDATMS,NDIMTS,ND +IMNF,NDIMDM,IDCODE,NDIMIV,IJOB) C C THIS ROUTINE READS THE INPUT TIME SERIES. C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION T(NDIMTS,NDIMNF),F(NDIMTS,NDIMNF),DATUMS(NDIMTS,NDIMDM) DIMENSION NDATMS(NDIMTS),IV(12),NF(NDIMTS),IDCODE(NDIMIV) IF(IJOB.GT.1) GOTO 6 READ(5,1001) FACTOR, CONST IF(FACTOR.EQ.0.D0) FACTOR=1.D0 PRINT 1004,FACTOR,CONST,ISER IF(ISER.EQ.1) READ(5,1003) (IDCODE(I),I=1,NDIMIV) 6 NF(ISER)=0 IPREV=0 KK=0 J=0 INTNO=0 1 READ(5,1002) IE,(IV(K),K=1,12) IF(IE.EQ.999) GO TO 5 DO 4 K=1,12 KK=KK+1 IF(IV(K).NE.0) GO TO 2 IPREV=0 GO TO 4 2 IF(IPREV.EQ.1) GO TO 3 IPREV=1 INTNO=INTNO+1 IF(IDCODE(INTNO).NE.0) GO TO 3 J=J+1 DATUMS(ISER,J)=DFLOAT(KK) 3 NF(ISER)=NF(ISER)+1 F(ISER,NF(ISER))=IV(K)*FACTOR+CONST T(ISER,NF(ISER))=DFLOAT(KK) 4 CONTINUE GO TO 1 5 NDATMS(ISER)=J 1001 FORMAT(2F10.3) 1002 FORMAT(I3,5X,12I6) 1003 FORMAT(80I1) 1004 FORMAT(' ',10X,'FACTOR= ',D13.6,' CONSTANT= ',D13.6,' FOR TIME S +ERIES # ',I3,//) RETURN END SUBROUTINE RNADD(I,J,RN,NDIMUK,ICNT,ITERM,NIVL,NDIMIV,TA,TB,FOFREQ +,NDIMFF,NDIMNF,N,KIVL,NDIMDM,T,F,NDIMTS,FPFAZ) C C THIS ROUTINE COMPUTES CONTRIBUTIONS TO THE NORMAL EQUATIONS C COEFFICIENT MATRIX FOR SYSTEMATIC NOISE. C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION RN(NDIMUK,NDIMUK),ICNT(2,7),ITERM(NDIMUK),TA(NDIMIV) DIMENSION TB(NDIMIV),FOFREQ(NDIMFF),KIVL(NDIMDM,4) DIMENSION T(NDIMTS,NDIMNF),F(NDIMTS,NDIMNF),FPFAZ(NDIMFF,7) II=ITERM(I) IJ=ITERM(J) IF(II.EQ.3.OR.II.EQ.4) OMEGA1=FOFREQ(ICNT(1,3)) IF(IJ.EQ.3.OR.IJ.EQ.4) OMEGA2=FOFREQ(ICNT(2,3)) IF(II.NE.5) GO TO 2 FAZ1=FPFAZ(ICNT(1,5),2) OMEGA1=FOFREQ(ICNT(1,5)) 2 IF(IJ.NE.5) GO TO 3 FAZ2=FPFAZ(ICNT(2,5),2) OMEGA2=FOFREQ(ICNT(2,5)) 3 GO TO(5,60,75,110,50,50,50),II 5 GO TO(10,20,30,40,50,50,50),IJ C C CONTRIBUTION TO RN FOR DATUM BIASES C 10 IF(I.NE.J) RETURN K1=KIVL(I,3) K2=KIVL(I,4) DO 15 K=K1,K2 15 RN(I,J)=RN(I,J)+TB(K)-TA(K)+1.D0 RETURN C C CONTRIBUTION TO RN FOR DATUM BIAS AND LINEAR TREND 20 K1=KIVL(I,3) K2=KIVL(I,4) DO 25 K=K1,K2 25 RN(I,J)=RN(I,J)+(TB(K)*(TB(K)+1.D0)-TA(K)*(TA(K)-1.D0))/2.D0 RETURN C C CONTRIBUTION TO RN FOR DATUM BIAS AND COS (TRIG TERM) C 30 K1=KIVL(I,3) K2=KIVL(I,4) OMEGA=FOFREQ(ICNT(2,3)) DO 35 K=K1,K2 N1=IDINT(TA(K)) N2=IDINT(TB(K)) CALL COSSUM(N1,N2,OMEGA,SUM) 35 RN(I,J)=RN(I,J)+SUM RETURN C C CONTRIBUTION TO RN FOR DATUM BIAS AND SIN (TRIG TERM) C 40 K1=KIVL(I,3) K2=KIVL(I,4) OMEGA=FOFREQ(ICNT(2,4)) DO 45 K=K1,K2 N1=IDINT(TA(K)) N2=IDINT(TB(K)) CALL SINSUM(N1,N2,OMEGA,SUM) 45 RN(I,J)=RN(I,J)+SUM RETURN C C CONTRIBUTION TO RN FOR ALL TERMS WITH NO EFFICIENT ALGORITHM C 50 DO 55 K=1,N DESGN1=DESIGN(K,I,II,KIVL,NDIMDM,T,F,NDIMTS,NDIMNF,OMEGA1,ICNT(1,6 +),ICNT(1,7),FAZ1) DESGN2=DESIGN(K,J,IJ,KIVL,NDIMDM,T,F,NDIMTS,NDIMNF,OMEGA2,ICNT(2,6 +),ICNT(2,7),FAZ2) 55 RN(I,J)=RN(I,J)+DESGN1*DESGN2 RETURN 60 GO TO(190,65,50,50,50,50,50),IJ C C CONTRIBUTION TO RN FOR LINEAR TREND C 65 DO 70 K=1,NIVL 70 RN(I,J)=RN(I,J)+(TB(K)*(TB(K)+1.D0)*(2.D0*TB(K)+1.D0)-(TA(K)-1.D0) +*TA(K)*(2.D0*TA(K)-1.D0))/6.D0 RETURN 75 GO TO(190,190,80,95,50,50,50),IJ C C CONTRIBUTION TO RN FOR COS, COS OF TRIG TERMS C 80 DO 85 K=1,NIVL N1=IDINT(TA(K)) N2=IDINT(TB(K)) CALL CCSUM(N1,N2,OMEGA1,OMEGA2,SUM) 85 RN(I,J)=RN(I,J)+SUM RETURN C C CONTRIBUTION TO RN FOR COS , SIN OF TRIG TERMS C 95 DO 100 K=1,NIVL N1=IDINT(TA(K)) N2=IDINT(TB(K)) CALL SCSUM(N1,N2,OMEGA2,OMEGA1,SUM) 100 RN(I,J)=RN(I,J)+SUM RETURN C C CONTRIBUTION TO RN FOR SIN, COS OF TRIG TERM C 105 TEMP=OMEGA2 OMEGA2=OMEGA1 OMEGA1=TEMP GO TO 95 110 GO TO(190,190,105,115,50,50,50),IJ C C CONTRIBUTION TO RN FOR SIN, SIN OF TRIG TERM C 115 DO 120 K=1,NIVL N1=IDINT(TA(K)) N2=IDINT(TB(K)) CALL SSSUM(N1,N2,OMEGA1,OMEGA2,SUM) 120 RN(I,J)=RN(I,J)+SUM 190 RETURN END SUBROUTINE SCSUM(N1,N2,OMEGA1,OMEGA2,SUM) C C THIS ROUTINE COMPUTES THE SUM FOR K = N1 TO N2 OF C SIN(OMEGA1*K)*COS(OMEGA2*K) FOR OMEGA1.NE.2*K*PI.AND OMEGA2. C NE.2*K*PI C C ROUTINE SINSUM IS CALLED C IMPLICIT REAL*8 (A-H,O-Z) OMEGAM=OMEGA1-OMEGA2 OMEGAP=OMEGA1+OMEGA2 CALL SINSUM(N1,N2,OMEGAM,SUM1) CALL SINSUM(N1,N2,OMEGAP,SUM2) SUM=(SUM1+SUM2)/2.D0 RETURN END SUBROUTINE SINSUM(N1,N2,OMEGA,SUM) C C THIS SUBROUTINE COMPUTES THE SUM FOR K=N1 TO N2 OF SIN(K*OMEGA) C FOR OMEGA.NE.2*K*PI C IMPLICIT REAL*8 (A-H,O-Z) N=N2-N1+1 R1=DFLOAT(N1) RN=DFLOAT(N) ALPHA=(R1-1.D0)*OMEGA SNX2=DSIN(RN*OMEGA/2.D0) SX2=DSIN(OMEGA/2.D0) IF(SX2.EQ.0.D0) GO TO 1 SAN1X2=DSIN(ALPHA+(RN+1.D0)*OMEGA/2.D0) SUM=SNX2*SAN1X2/SX2 RETURN 1 SUM=0.D0 RETURN END SUBROUTINE SPCOMP(NSPFRQ,SPFREQ,NONSIM,NF,NSYS,RN,RU,I1,I2,T,F,NCO +S,NDATMS,LNTRND,NFFREQ,FOFREQ,FPFAZ,KIVL,NDIMDM,NDIMTS,NDIMNF,SP,C +P,XN,XL,CNP,SNP,CLP,SLP,SPMQ,SPPQ,NIVL,IVL,SVAR,NDIMSP,NDIMUK,NDIM +FF,NDIMIV,S,NSYSNF,NSYSUD) C C THIS ROUTINE COMPUTES THE SPECTRUM. C IMPLICIT REAL*8(A-H,O-Z) DIMENSION SPFREQ(NDIMSP),NF(NDIMTS),RN(NDIMUK,NDIMUK),RU(NDIMUK) DIMENSION T(NDIMTS,NDIMNF),F(NDIMTS,NDIMNF),NDATMS(NDIMTS) DIMENSION FOFREQ(NDIMFF),FPFAZ(NDIMFF,7),KIVL(NDIMDM,4),SP(NDIMFF) DIMENSION CP(NDIMFF),XN(NDIMIV),XL(NDIMIV),CNP(NDIMIV,NDIMFF) DIMENSION SNP(NDIMIV,NDIMFF),CLP(NDIMIV,NDIMFF),SLP(NDIMIV,NDIMFF) DIMENSION SPMQ(NDIMFF),SPPQ(NDIMFF),IVL(NDIMIV),S(NDIMSP) N=NF(1) C C LOOP THROUGH FOR EACH SPECTRAL VALUE C DO 130 I=1,NSPFRQ OMEGA=SPFREQ(I) IF(NONSIM.EQ.1.OR.NSYS.EQ.0) GO TO 81 C C ZERO CORRESPONDING PARTS OF NORMAL EQUATIONS. C DO 80 J=1,NSYS RN(J,I1)=0.D0 80 RN(J,I2)=0.D0 81 RU(I1)=0.D0 RU(I2)=0.D0 CC=0.D0 CS=0.D0 SS=0.D0 C C COMPUTE CORRESPONDING PART OF CONSTANT VECTOR C DO 85 J=1,N WT=OMEGA*T(1,J) COSWT=DCOS(WT) SINWT=DSIN(WT) RU(I1)=RU(I1)+F(1,J)*COSWT RU(I2)=RU(I2)+F(1,J)*SINWT C C COMPUTE CONTRIBUTION TO PART OF COEFFICIENT MATRIX FROM FORCED C FREQUENCIES WITH FORCED PHASE LAG. C IF(NCOS.EQ.0) GO TO 86 DO 82 L=1,NCOS K=NDATMS(1)+LNTRND+2*NFFREQ+L OMEG=FOFREQ(NFFREQ+L) FAZ=FPFAZ(NFFREQ+L,2) DESIG =DESIGN(J,K,5,KIVL,NDIMDM,T,F,NDIMTS,NDIMNF,OMEG,1,1,FAZ) RN(K,I1)=RN(K,I1)+DESIG *COSWT 82 RN(K,I2)=RN(K,I2)+DESIG *SINWT 86 CONTINUE C C . . . FOR FORCED NUMERICAL FUNCTIONS . . . C IF(NSYSNF.EQ.0) GO TO 84 DO 83 L=1,NSYSNF K=NDATMS(1)+LNTRND+2*NFFREQ+NCOS+L DESIG=F(L+1,J) RN(K,I1)=RN(K,I1)+DESIG*COSWT 83 RN(K,I2)=RN(K,I2)+DESIG*SINWT 84 CONTINUE C C . . . FOR FORCED USER DEFINED FUNCTIONS . . . C IF(NSYSUD.EQ.0) GO TO 87 DO 88 L=1,NSYSUD K=NDATMS(1)+LNTRND+2*NFFREQ+NCOS+NSYSNF+L DESIG=FUNCTN(T(1,J),L) RN(K,I1)=RN(K,I1)+DESIG*COSWT 88 RN(K,I2)=RN(K,I2)+DESIG*SINWT 87 CONTINUE 85 CONTINUE Q=OMEGA/2.D0 SQ=DSIN(Q) CQ=DCOS(Q) IF(NFFREQ.EQ.0.OR.NONSIM.EQ.1) GO TO 95 DO 90 J=1,NFFREQ SPCQ=SP(J)*CQ CPSQ=CP(J)*SQ SPMQ(J)=SPCQ-CPSQ SPPQ(J)=SPCQ+CPSQ IF(DABS(SPMQ(J)).LT.1D-30) SPMQ(J)=DSIGN(1D-30,SPMQ(J)) 90 CONTINUE 95 DO 115 J=1,NIVL XNQ=XN(J)*Q XLQ=XL(J)*Q SNQ=DSIN(XNQ) CNQ=DCOS(XNQ) SLQ=DSIN(XLQ) CLQ=DCOS(XLQ) SNQCNQ=SNQ*CNQ SNQCLQ=SNQ*CLQ SNQSLQ=SNQ*SLQ CC=CC+SNQCNQ*(CLQ*CLQ-SLQ*SLQ) CS=CS+SNQCNQ*SLQ*CLQ IF(NONSIM.EQ.1.OR.NSYS.EQ.0) GO TO 115 IF(NDATMS(1).EQ.0) GO TO 100 K=IVL(J) RN(K,I1)=RN(K,I1)+SNQCLQ/SQ RN(K,I2)=RN(K,I2)+SNQSLQ/SQ C C . . . FOR FORCED LINEAR TREND . . . C 100 IF(LNTRND.EQ.0) GO TO 105 K=NDATMS(1)+1 RN(K,I1)=RN(K,I1)+(-SNQSLQ*CQ/SQ+XN(J)*CNQ*SLQ+XL(J)*SNQCLQ)/2.D0/ +SQ RN(K,I2)=RN(K,I2)+(SNQCLQ*CQ/SQ-XN(J)*CNQ*CLQ+XL(J)*SNQSLQ)/2.D0/S +Q C C. . . AND FINISH FOR FORCED FREQUENCIES . . . C 105 IF(NFFREQ.EQ.0) GO TO 115 DO 110 L=1,NFFREQ K=NDATMS(1)+LNTRND+2*L-1 P1=SNP(J,L)*CNQ P2=CNP(J,L)*SNQ SNPMQ=(P1-P2)/SPMQ(L) SNPPQ=(P1+P2)/SPPQ(L) P1=SLP(J,L)*CLQ P2=CLP(J,L)*SLQ SINM=(P1-P2)*SNPMQ SINP=(P1+P2)*SNPPQ P1=CLP(J,L)*CLQ P2=SLP(J,L)*SLQ COSM=(P1+P2)*SNPMQ COSP=(P1-P2)*SNPPQ RN(K,I1)=RN(K,I1)+(COSP+COSM)/2.D0 RN(K+1,I1)=RN(K+1,I1)+(SINP+SINM)/2.D0 RN(K,I2)=RN(K,I2)+(SINP-SINM)/2.D0 110 RN(K+1,I2)=RN(K+1,I2)+(-COSP+COSM)/2.D0 115 CONTINUE S(I)=0.D0 P1=DFLOAT(N)/2.D0 P2=CC/(2.D0*SQ*CQ) SS=P1-P2 CC=P1+P2 CS=CS/(SQ*CQ) IF(NONSIM.EQ.1.OR.NSYS.EQ.0) GO TO 125 UNU=0.D0 UNV=0.D0 VNV=0.D0 C C COMPUTE BILINEAR FORM FOR SPECTRAL VALUE COMPUTATION. C DO 120 J=1,NSYS DO 120 K=1,NSYS P1=RN(J,I1)*RN(J,K) UNU=UNU+P1*RN(K,I1) UNV=UNV+P1*RN(K,I2) 120 VNV=VNV+RN(J,I2)*RN(J,K)*RN(K,I2) CCUNU=CC-UNU CSUNV=CS-UNV SSVNV=SS-VNV DET=CCUNU*SSVNV-CSUNV*CSUNV IF(DABS(DET).LT.1.D-30) GO TO 130 C C COMPUTE SPECTRAL VALUE C S(I)=100.D0*(SSVNV*RU(I1)**2-2.D0*CSUNV*RU(I1)*RU(I2)+CCUNU*RU(I2) +**2)/DET/SVAR GO TO 130 125 DET=CC*SS-CS*CS IF(DABS(DET).LT.1.D-30) GO TO 130 S(I)=100.D0*(SS*RU(I1)**2-2.D0*CS*RU(I1)*RU(I2)+CC*RU(I2)**2)/DET/ +SVAR 130 CONTINUE RETURN END SUBROUTINE SPKVAL(RPL,RPS,NSPFRQ,SPFREQ,NDIMSP) C C THIS ROUTINE COMPUTES SPECTRAL FREQUENCIES AND STORES THEM IN SPFREQ. C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION SPFREQ(NDIMSP) PI2=3.141592653589793D0*2.D0 RWL=PI2/RPL RN1=NSPFRQ-1 IF(NSPFRQ.EQ.1) RN1=1.D0 DW=(PI2/RPS-RWL)/RN1 WT=RWL DO 1 I=1,NSPFRQ SPFREQ(I)=WT 1 WT=WT+DW RETURN END SUBROUTINE SPLOT(P,S,NW,IB,NDIMSP,TUNITS,RANSP) C C THIS ROUTINE PLOTS SPECTRUM ON THE LINE PRINTER. C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(NDIMSP),S(NDIMSP),IPLOT(88) DATA IBLNK/' '/,ISTAR/'*'/,IDOT/':'/ PRINT 1001,IB,NW,P(1),P(NW),TUNITS,TUNITS SMAX=S(1) DO 1 I=2,NW IF(S(I).GT.SMAX) SMAX=S(I) 1 CONTINUE KRS=IDINT(RANSP*88.D0/SMAX) IF(KRS.LT.1) KRS=1 IF(KRS.GT.88) KRS=88 DO 5 I=1,88 5 IPLOT(I)=IBLNK DO 10 I=1,NW KS=IDINT(S(I)*88.D0/SMAX) IF(KS.LT.1)KS=1 IF(KS.GT.88)KS=88 IPLOT(KS)=ISTAR PRINT 1002,P(I),S(I),(IPLOT(J),J=1,88) IPLOT(KS)=IBLNK IPLOT(KRS)=IDOT PRINT 1003,(IPLOT(J),J=1,88) 10 IPLOT(KRS)=IBLNK 1001 FORMAT(' ',///,' ',12X,'SPECTRAL BAND ',I5,/,' ',10X,I5,' SPECTRAL + VALUES BETWEEN',F14.6,' AND',F14.6,' ',A8,//,' ',5X,A8,3X,'% VAR' +,/) 1002 FORMAT(' ',F13.5,1X,F8.2,88(A1)) 1003 FORMAT('+',22X,88(A1)) RETURN END SUBROUTINE SSSUM(N1,N2,OMEGA1,OMEGA2,SUM) C C THIS ROUTINE COMPUTES THE SUM FOR K=N1 TO N2 OF C SIN(OMEGA1*K)*SIN(OMEGA2*K) FOR OMEGA1.NE.2*K*PI. AND OMEGA2. C NE.2*K*PI C C ROUTINE COSSUM IS CALLED C IMPLICIT REAL*8 (A-H,O-Z) OMEGAM=OMEGA1-OMEGA2 OMEGAP=OMEGA1+OMEGA2 CALL COSSUM(N1,N2,OMEGAM,SUM1) CALL COSSUM(N1,N2,OMEGAP,SUM2) SUM=(SUM1-SUM2)/2.D0 RETURN END SUBROUTINE SYSTAT(NF,NDIMTS,NSYS,RN,RU,ITERM,NDIMUK,NIVL,NDIMIV,TA +,TB,FOFREQ,NDIMFF,NDIMNF,T,F,NFFREQ,KIVL,NDIMDM,FPFAZ,D,X,SVAR,VFA +C,NDATMS,LNTRND,FORPER,NSYSNF,NSYSUD,FUNITS,TUNITS,RHOMAX,NCOS,DAT +UMS,NPRRS,FCRIT,FNFUNT) C C THIS ROUTINE COMPUTES SYSTEMATIC NOISE AND SOME STATISTICS. C IMPLICIT REAL*8(A-H,O-Z) DIMENSION NF(NDIMTS),RN(NDIMUK,NDIMUK),RU(NDIMUK),ITERM(NDIMUK) DIMENSION D(NDIMUK),X(NDIMUK),TA(NDIMIV),TB(NDIMIV),FOFREQ(NDIMFF) DIMENSION FORPER(NDIMFF),T(NDIMTS,NDIMNF),F(NDIMTS,NDIMNF) DIMENSION KIVL(NDIMDM,4),FPFAZ(NDIMFF,7),DATUMS(NDIMTS,NDIMDM) DIMENSION ICNT(2,7),NDATMS(NDIMTS),FNFUNT(NDIMTS) N=NF(1) IDF=N-NSYS PRINT 1001,IDF IF(IDF.LT.0) STOP RANSP=200.D0/DFLOAT(IDF) PRINT 1004,RANSP DF=IDF-2 FCRIT=100.D0/(1.D0+1.D0/(0.05D0**(-2.D0/DF)-1.D0)) PRINT 1006,FCRIT I1=1 I2=NSYS IF(I2.EQ.0) GO TO 40 CALL NORMEQ(I1,I2,RN,NDIMUK,RU,ITERM,NIVL,NDIMIV,TA,TB,FOFREQ,NDIM +FF,NDIMNF,NF,T,NDIMTS,F,NFFREQ,KIVL,NDIMDM,FPFAZ) CALL XSOL(RN,NSYS,RU,D,X,NDIMUK,1,0) DO 5 K=1,7 5 ICNT(1,K)=0 ICNT(1,5)=NFFREQ ICNT(1,6)=1 DO 15 J=1,NSYS IJ=ITERM(J) ICNT(1,IJ)=ICNT(1,IJ)+1 IF(IJ.EQ.3.OR.IJ.EQ.4) OMEGA=FOFREQ(ICNT(1,3)) IF(IJ.NE.5) GO TO 10 OMEGA=FOFREQ(ICNT(1,5)) FAZ=FPFAZ(ICNT(1,5),2) 10 ISER=ICNT(1,6) IFUNC=ICNT(1,7) DO 15 I=1,N F(1,I)=F(1,I)-X(J)*DESIGN(I,J,IJ,KIVL,NDIMDM,T,F,NDIMTS,NDIMNF,OME +GA,ISER,IFUNC,FAZ) 15 CONTINUE 40 SVAR=0.D0 DO 20 I=1,N 20 SVAR=SVAR+F(1,I)*F(1,I) IF(IDF.EQ.0) GO TO 25 VFAC=SVAR/DFLOAT(IDF) SVF=DSQRT(VFAC) PRINT 1002,SVF GO TO 30 25 VFAC=1.D0 30 IF(NSYS.GT.0) CALL PAMPR(NDATMS,NDIMTS,LNTRND,NFFREQ,FORPER,NDIMFF +,NSYSNF,NSYSUD,X,NDIMUK,FPFAZ,FUNITS,VFAC,RN,TUNITS,RHOMAX,NCOS,NS +YS,FNFUNT) IF(NPRRS.GT.0) GO TO 35 PRINT 1003,TUNITS,FUNITS CALL FPLOT(1,T,F,NF,DATUMS,NDATMS,NDIMTS,NDIMNF,NDIMDM) PRINT 1005 1001 FORMAT(' ',///,' ','DEGREES OF FREEDOM = ',I6,/) 1002 FORMAT(' ','SQUARE ROOT OF VARIANCE FACTOR = ',D12.4,///) 1003 FORMAT(' ','RESIDUAL SERIES: (TIME UNITS = ',A8,' ,TIME SERIES UNI +TS = ',A8,' )',/) 1004 FORMAT(' ','LEVEL OF RANDOM NOISE SPECTRUM = ',F7.3,/) 1005 FORMAT(' ',110('-'),/) 1006 FORMAT(' ','95% CRITICAL VALUE FOR SPECTRUM = ',F7.3,/) 35 RETURN END SUBROUTINE TRIGPK(NFFREQ,FOFREQ,NIVL,SNP,CNP,SLP,CLP,NDIMFF,NDIMIV +,SP,CP,XN,XL) C C THIS ROUTINE COMPUTES AND STORES TRIG FUNCTIONS OF FORCED FREQUENCIES C FOR SPECTRAL COMPUTATIONS. C IMPLICIT REAL*8(A-H,O-Z) DIMENSION FOFREQ(NDIMFF),SNP(NDIMIV,NDIMFF),CNP(NDIMIV,NDIMFF) DIMENSION SLP(NDIMIV,NDIMFF),CLP(NDIMIV,NDIMFF),SP(NDIMFF) DIMENSION CP(NDIMFF),XN(NDIMIV),XL(NDIMIV) DO 1 I=1,NFFREQ PK=FOFREQ(I)/2.D0 SP(I)=DSIN(PK) CP(I)=DCOS(PK) DO 1 J=1,NIVL XNK=XN(J)*PK XLK=XL(J)*PK SNP(J,I)=DSIN(XNK) CNP(J,I)=DCOS(XNK) SLP(J,I)=DSIN(XLK) 1 CLP(J,I)=DCOS(XLK) RETURN END SUBROUTINE TSGEN(NGEN,ISER,IDAT,NDIMDM,DATUMS,NDIMTS,NDIMNF,T,F,ND +IMFF,FOFREQ,ACOEF,BCOEF,PI2,DATS,NF,NDATMS,IDCODE,NDIMIV) C C THIS ROUTINE READS THE INPUT PARAMETERS FOR GENERATION OF TIME SERIES C AND GENERATES TIME SERIES. C IMPLICIT REAL*8(A-H,O-Z) DIMENSION IDAT(NDIMIV,2),DATUMS(NDIMTS,NDIMDM),T(NDIMTS,NDIMNF) DIMENSION F(NDIMTS,NDIMNF),FOFREQ(NDIMFF),ACOEF(NDIMFF) DIMENSION BCOEF(NDIMFF),DATS(NDIMTS,NDIMIV),NF(NDIMTS) DIMENSION NDATMS(NDIMTS),IDCODE(NDIMIV) READ(5,1001) LT,COEFLT,NBASE,RNOISE DO 1 I=1,NGEN READ(5,1002) (IDAT(I,J),J=1,2),DATS(ISER,I) 1 CONTINUE IF(ISER.EQ.1) READ(5,1004) (IDCODE(I),I=1,NGEN) NPER=0 DO 2 I=1,NDIMFF READ(5,1003) FOFREQ(I),ACOEF(I),BCOEF(I) IF(FOFREQ(I).EQ.0.D0) GO TO 3 FOFREQ(I)=PI2/FOFREQ(I) NPER=NPER+1 2 CONTINUE 3 NF(ISER)=0 KINT=0 DO 7 I=1,NGEN ISTART=IDAT(I,1) IF(ISER.NE.1) GO TO 11 IF(IDCODE(I).NE.0) GO TO 12 KINT=KINT+1 DATUMS(1,KINT)=DFLOAT(ISTART) GO TO 12 11 DATUMS(ISER,I)=DFLOAT(ISTART) 12 IEND=IDAT(I,2) DO 7 J=ISTART,IEND NF(ISER)=NF(ISER)+1 T(ISER,NF(ISER))=DFLOAT(J) F(ISER,NF(ISER))=DATS(ISER,I)+COEFLT*DFLOAT(J)+RNOISE*RND(5) IF(NPER.EQ.0) GO TO 5 RJ=DFLOAT(J) DO 4 K=1,NPER WT=RJ*FOFREQ(K) F(ISER,NF(ISER))=F(ISER,NF(ISER))+ACOEF(K)*DCOS(WT)+BCOEF(K)*DSIN( +WT) 4 CONTINUE 5 CONTINUE IF(NBASE.EQ.0) GO TO 7 DO 6 K=1,NBASE F(ISER,NF(ISER))=F(ISER,NF(ISER))+FUNCTN(T(ISER,NF(ISER)),K) 6 CONTINUE 7 CONTINUE IF(ISER.NE.1) GO TO 9 ISUM=0 DO 8 I=1,NGEN 8 ISUM=ISUM+IDCODE(I) NDATMS(ISER)=NGEN-ISUM GO TO 10 9 NDATMS(ISER)=NGEN 10 CONTINUE CALL PARGEN(NGEN,LT,COEFLT,NBASE,IDAT,DATS,FOFREQ,ACOEF,BCOEF,NPER +,NDIMDM,NDIMTS,NDIMFF,PI2,RNOISE,ISER,IDCODE,NDIMIV) 1001 FORMAT(I5,F10.6,I5,F10.6) 1002 FORMAT(2I5,F10.6) 1003 FORMAT(3F10.6) 1004 FORMAT(80I1) RETURN END SUBROUTINE SPKVAL(RPL,RPS,NSPFRQ,SPFREQ,NDIMSP) C C THIS ROUTINE COMPUTES SPECTRAL FREQUENCIES AND STORES THEM IN SPFREQ. C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION SPFREQ(NDIMSP) PI2=3.141592653589793D0*2.D0 RWL=PI2/RPL RN1=NSPFRQ-1 IF(NSPFRQ.EQ.1) RN1=1.D0 DW=(PI2/RPS-RWL)/RN1 WT=RWL DO 1 I=1,NSPFRQ SPFREQ(I)=WT 1 WT=WT+DW RETURN END SUBROUTINE XSOL(T,N,B,D,X,NR,JSTART,JINV) C C THIS ROUTINE SOLVES THE SYSTEM N*X=U FOR POSITIVE DEFINITE SYMMETRIC C N (STORED IN T). USED FOR SYSTEMATIC NOISE COMPUTATION. C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION T(NR,NR),D(NR),B(NR),X(NR) C C COMPUTE CHOLESKI SQUARE ROOT IN PLACE C DO 4 J=JSTART,N DO 4 I=JSTART,J IF(I.EQ.JSTART) GO TO 2 M=I-1 SUM=0.D0 DO 1 K=JSTART,M 1 SUM=SUM+T(K,I)*T(K,J) T(I,J)=T(I,J)-SUM 2 IF(I.EQ.J) GO TO 3 T(I,J)=T(I,J)/T(I,I) GO TO 4 3 T(I,I)=DSQRT(T(I,I)) 4 CONTINUE C C COMPUTE SOLUTION IF JINV = 0 C IF(JINV.GT.0) GO TO 9 D(1)=B(1)/T(1,1) IF(N.EQ.1) GO TO 15 DO 6 I=2,N SUM=0.D0 K=I-1 DO 5 J=1,K 5 SUM=SUM+T(J,I)*D(J) 6 D(I)=(B(I)-SUM)/T(I,I) 15 X(N)=D(N)/T(N,N) IF(N.EQ.1) GO TO 9 M=N-1 DO 8 I=1,M SUM=0.D0 J=N-I+1 L=N-I DO 7 K=J,N 7 SUM=SUM+T(L,K)*X(K) 8 X(L)=(D(L)-SUM)/T(L,L) C C COMPUTE INVERSE C 9 DO 12 J=JSTART,N DO 12 I=JSTART,J IF(I.LT.J) GO TO 10 T(J,J)=1.D0/T(J,J) GO TO 12 10 SUM=0.D0 M=J-1 DO 11 K=I,M 11 SUM=SUM-T(I,K)*T(K,J) T(I,J)=SUM/T(J,J) 12 CONTINUE C DO 14 J=JSTART,N DO 14 I=JSTART,J SUM=0.D0 DO 13 K=J,N 13 SUM=SUM+T(I,K)*T(J,K) T(I,J)=SUM 14 T(J,I)=SUM RETURN END