C PROGRAM 9.1 SMOOTH C C ... Prediction and Interpolation of Time Series ... C C Inputs: C NFE: End point of filtering C NPE: End point of prediction C M: Order of the AR model C TAU2: Innovation variance of the AR model C AR(I): AR coefficients (II=1,M) C OUTMIN,OUTMAX: Lower and upper limits of observations C NMISS: Number of missed intervals C N0(I): Start position of missed intervals (I=1,NMISS) C NN(I): Number of missed observations (I=1,NMISS) C Inputs required in subroutine READTS C TITLE: Title of the data set C N: Data length C Y(I): Tiem series (I=1,N) C Parameters: C IDEV: Input device specification C NMAX: Adjustable dimension of Y, YMISS (NMAX >= N) C MJ: Highest dimension of the state (MJ >= M) C K: Dimension of the system noise C ISW: =1 (R is specified) C @TEST.FILTER2: SEP.08,1990, SEP.02,1992 C MODIFIED 2/15/93 C PARAMETER( NMAX=1000,MJ=40,K=1,ISW=1,IDEV=1 ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(NMAX), YMISS(NMAX) DIMENSION F(MJ,MJ), G(MJ), H(MJ), Q(K,K) DIMENSION XPS(MJ,NMAX), XFS(MJ,NMAX), XSS(MJ,NMAX) DIMENSION VPS(MJ,MJ,NMAX), VFS(MJ,MJ,NMAX), VSS(MJ,MJ,NMAX) DIMENSION XF(MJ), VF(MJ,MJ) DIMENSION N0(10), NN(10), AR(MJ) C READ( 5,* ) NFE, NPE READ( 5,* ) M READ( 5,* ) TAU2 READ( 5,* ) (AR(I),I=1,M) READ( 5,* ) OUTMIN, OUTMAX READ( 5,* ) NMISS DO 10 I=1,NMISS 10 READ( 5,* ) N0(I), NN(I) NS = 1 C C ... Read Time Series ... C CALL READTS( IDEV,Y,N ) CALL MOMENT( Y,N,YMEAN,YVAR ) C C ... SET MISSING OBSERVATIONS ... CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 20 I=1,N 20 YMISS(I) = Y(I) - YMEAN DO 30 J=1,NMISS DO 30 I=1,NN(J) 30 YMISS(N0(J)+I-1) = OUTMIN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ... Set State Space Model ... C CALL SETFGH( M,MJ,K,AR,TAU2,F,G,H,Q,R ) C C ... Set Initial State ... C CALL ISTATE( M,MJ,0.0D0,YVAR,XF,VF ) C C ... Kalman Filtering ... C CALL FILTER( YMISS,XF,VF,F,G,H,Q,R,M,K,ISW,NS,NFE,NPE,MJ,NMAX, * NMAX,OUTMIN,OUTMAX,VFS,VPS,XFS,XPS,FF,SIG2 ) C C ... Fixed Interval Smoothing ... C CALL SMOOTH( F,M,MJ,NDIM,NS,NFE,NPE,VFS,VPS,XFS,XPS, * VSS,XSS ) FLK = -FF AIC = -2*FF + 2*(M+1) C C ... Plot and Print ... C CALL PTSTAT( Y,M,TAU2,FLK,AIC,XSS,VSS,N,NFE,NPE,MJ, * YMEAN,NMISS,N0,NN ) CALL PRSTAT( M,TAU2,FLK,AIC,XSS,VSS,YMEAN,MJ,NFE,NPE, * NMISS,N0,NN ) C STOP E N D SUBROUTINE FILTER( Y,XF,VF,F,G,H,Q,R,M,K,ISW,NS,NFE,NPE,MJ,NMAX, * NDIM,OUTMIN,OUTMAX,VFS,VPS,XFS,XPS,FF,SIG2 ) C C ... Kalman Filter (General Form, L=1) ... C C Inputs: C Y: time series C NS: Start position of filtering C NFE: End position of filtering C NPE: End position of prediction C XF: Initial state vector C VF: Initial covariance matrix C M: Dimension of the state vector C K: Dimension of the system noise C F: M*M matrix C G: M*K matrix C H: M vector C Q: K*K matrix, system noise covariance C R: observation variance C ISW: = 0; R will be estimated C = 1; R is given C MJ: Adjustable dimension of XF, VF C NDIM: Adjustable dimension of XFS, XPS, VFS, VPS C = 0 XF, XP, VF, VP are not stored C > 0 They are stored for smoothing C NMAX Adjustable dimension of Y C OUTMIN: Lower limit for detecting outliers C OUTMAX: Upper limit for detecting outliers C Outputs: C VFS: Covariance matrices of the filter C VPS: Covariance matrices of the predictor C XFS: Mean vectors of the filter C XPS: Mean vectors of the predictor C FF: Log likelihood C SIG2: Estimated observational noise variance C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(NMAX) DIMENSION F(MJ,MJ), G(MJ,K), H(MJ), Q(K,K) DIMENSION XF(MJ), VF(MJ,MJ), XP(40), VP(40,40) DIMENSION XFS(MJ,NDIM), XPS(MJ,NDIM) DIMENSION VFS(MJ,MJ,NDIM), VPS(MJ,MJ,NDIM) DIMENSION WRK(40,40), VH(40), GAIN(40) DATA PI /3.1415926535D0/ C SIG2 = 0.0D0 SDET = 0.0D0 NSUM = 0 C DO 500 II=NS,NPE C C ... ONE STEP AHEAD PREDICTION ... C DO 20 I=1,M SUM = 0.0D0 DO 10 J=1,M 10 SUM = SUM + F(I,J)*XF(J) 20 XP(I) = SUM C DO 40 I=1,M DO 40 J=1,M SUM = 0.0D0 DO 30 JJ=1,M 30 SUM = SUM + F(I,JJ)*VF(JJ,J) 40 WRK(I,J) = SUM C DO 60 I=1,M DO 60 J=1,M SUM = 0.0D0 DO 50 JJ=1,M 50 SUM = SUM + WRK(I,JJ)*F(J,JJ) 60 VP(I,J) = SUM C DO 80 I=1,M DO 80 J=1,K SUM = 0.0D0 DO 70 JJ=1,K 70 SUM = SUM + G(I,JJ)*Q(JJ,J) 80 WRK(I,J) = SUM C DO 100 I=1,M DO 100 J=1,M SUM = VP(I,J) DO 90 JJ=1,K 90 SUM = SUM + WRK(I,JJ)*G(J,JJ) 100 VP(I,J) = SUM C C ... FILTERING ... C IF( Y(II).GT.OUTMIN.AND.Y(II).LT.OUTMAX.AND. II.LE.NFE ) THEN C DO 210 I=1,M SUM = 0.0D0 DO 200 J=1,M 200 SUM = SUM + VP(I,J)*H(J) 210 VH(I) = SUM C PERR = Y(II) PVAR = R DO 220 I=1,M PERR = PERR - H(I)*XP(I) 220 PVAR = PVAR + H(I)*VH(I) C DO 250 I=1,M 250 GAIN(I) = VH(I)/PVAR C DO 290 I=1,M 290 XF(I) = XP(I) + GAIN(I)*PERR C DO 310 I=1,M DO 310 J=1,M 310 VF(I,J) = VP(I,J) - GAIN(I)*VH(J) C SIG2 = SIG2 + PERR**2/PVAR SDET = SDET + DLOG(PVAR) NSUM = NSUM + 1 C C ... MISSING OBSERVATION ... C ELSE DO 350 I=1,M XF(I) = XP(I) DO 350 J=1,M 350 VF(I,J) = VP(I,J) END IF C C ... SAVE MEAN AND COVARIANCE ... C IF( NDIM.GT.1 ) THEN DO 360 I=1,M XPS(I,II) = XP(I) XFS(I,II) = XF(I) DO 360 J=1,M VPS(I,J,II) = VP(I,J) 360 VFS(I,J,II) = VF(I,J) END IF C 500 CONTINUE SIG2 = SIG2/NSUM IF(ISW.EQ.0) FF = -0.5D0*(NSUM*(DLOG(PI*2*SIG2) + 1) + SDET) IF(ISW.EQ.1) FF = -0.5D0*(NSUM*(DLOG(PI*2) + SIG2) + SDET) C RETURN E N D SUBROUTINE SMOOTH( F,M,MJ,NDIM,NS,NFE,NPE,VFS,VPS,XFS,XPS, * VSS,XSS ) C C ... Fixed Interval Smoother (General Form) ... C C Inputs: C NS: Start position of filtering C NFE: End position of filtering C NPE: End position of prediction C M: Dimension of the state vector C F: M*M matrix C MJ: Adjustable dimension of XF, VF C NDIM: Adjustable dimension of XFS, XPS, VFS, VPS C NMAX Adjustable dimension of Y C VFS: Covariance matrices of the filter C VPS: Covariance matrices of the predictor C XFS: Mean vectors of the filter C XPS: Mean vectors of the predictor C Outputs: C VSS: Covariance matrices of the smoother C XSS: Mean vectors of the smoother C IMPLICIT REAL*8(A-H,O-Z) DIMENSION F(MJ,MJ) DIMENSION XS(40), VS(40,40), VP(40,40) DIMENSION XFS(MJ,NDIM), XPS(MJ,NDIM), XSS(MJ,NDIM) DIMENSION VFS(MJ,MJ,NDIM), VPS(MJ,MJ,NDIM), VSS(MJ,MJ,NDIM) DIMENSION WRK(40,40), SGAIN(40,40) C C C ... SMOOTHING ... C DO 10 II=NFE,NPE DO 10 I=1,M XSS(I,II) = XFS(I,II) DO 10 J=1,M 10 VSS(I,J,II) = VFS(I,J,II) DO 20 I=1,M XS(I) = XFS(I,NFE) DO 20 J=1,M 20 VS(I,J) = VFS(I,J,NFE) C DO 500 II=NFE-1,NS,-1 C NZERO = 0 DO 100 I=1,M 100 IF( VFS(I,I,II).GT.1.0D-12 ) NZERO = NZERO + 1 C IF( NZERO.EQ.0 ) THEN DO 110 I=1,M XS(I) = XFS(I,II) XSS(I,II) = XFS(I,II) DO 110 J=1,M VS(I,J) = VFS(I,J,II) 110 VSS(I,J,II) = VFS(I,J,II) C ELSE DO 410 I=1,M DO 410 J=1,M 410 VP(I,J) = VPS(I,J,II+1) C CALL GINVRS( VP,VDET,M,40 ) C DO 425 I=1,M DO 425 J=1,M SUM = 0.0D0 DO 420 IJ=1,M 420 SUM = SUM + VFS(I,IJ,II)*F(J,IJ) 425 WRK(I,J) = SUM C DO 440 I=1,M DO 440 J=1,M SUM = 0.0D0 DO 430 IJ=1,M 430 SUM = SUM + WRK(I,IJ)*VP(IJ,J) 440 SGAIN(I,J) = SUM C DO 450 I=1,M XS(I) = XFS(I,II) DO 450 J=1,M WRK(I,J) = 0.0D0 450 VS (I,J) = VFS(I,J,II) C DO 460 J=1,M DO 460 I=1,M 460 XS(I) = XS(I) + SGAIN(I,J)*(XSS(J,II+1) - XPS(J,II+1)) C DO 470 J=1,M DO 470 IJ=1,M DO 470 I=1,M 470 WRK(I,J)=WRK(I,J) + SGAIN(I,IJ)*(VSS(IJ,J,II+1)-VPS(IJ,J,II+1)) C DO 480 J=1,M DO 480 IJ=1,M DO 480 I=1,M 480 VS(I,J) = VS(I,J) + WRK(I,IJ)*SGAIN(J,IJ) DO 485 I=1,M 485 IF( VS(I,I).LT.0.0D0 ) VS(I,I) = 0.0D0 C DO 490 I=1,M XSS(I,II) = XS(I) DO 490 J=1,M 490 VSS(I,J,II) = VS(I,J) END IF C 500 CONTINUE C RETURN E N D SUBROUTINE GINVRS( A,DET,M,MJ ) C C ... Generalized inverse of a square matrix A ... C C Inputs: C A: M*M matrix C M: Dimension of A C MJ: Adjustable dimension of A C Outputs: C A: Generalize inverse of A C DET: Determinant of A C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(MJ,MJ), IND(50) EPS = 1.0D-10 C DO 10 I=1,M 10 IND(I) = I DO 60 L=1,M AMAX = 0.0D0 DO 20 I=L,M IF( A(IND(I),IND(I)).GT.AMAX ) THEN AMAX = A(IND(I),IND(I)) I0 = I END IF 20 CONTINUE IF( AMAX.GT.EPS*A(IND(1),IND(1)) ) THEN IMAX = IND(I0) DO 30 I=I0,L+1,-1 30 IND(I) = IND(I-1) IND(L) = IMAX LMAX = L DO 40 I=L+1,M A(IND(I),IMAX) = -A(IND(I),IMAX)/A(IMAX,IMAX) DO 40 J=L+1,M 40 A(IND(I),IND(J)) = A(IND(I),IND(J)) * + A(IND(I),IMAX)*A(IMAX,IND(J)) ELSE DO 50 I=L,M DO 50 J=L,M 50 A(IND(I),IND(J)) = 0.0D0 GO TO 70 END IF 60 CONTINUE C 70 DET = 1.0D0 DO 80 I=1,M C DET = DET*A(IND(I),IND(I)) C IF( A(IND(I),IND(I)).GT.EPS*AMAX ) THEN IF( A(IND(I),IND(I)).GT.0.0D0 ) THEN A(IND(I),IND(I)) = 1.0D0/A(IND(I),IND(I)) ELSE A(IND(I),IND(I)) = 0.0D0 END IF 80 CONTINUE C MS = MIN0( M-1,LMAX ) DO 200 L=MS,1,-1 DO 100 J=L+1,M SUM = 0.0D0 DO 90 I=L+1,M 90 SUM = SUM + A(IND(I),IND(L))*A(IND(I),IND(J)) 100 A(IND(L),IND(J)) = SUM SUM = A(IND(L),IND(L)) DO 110 I=L+1,M 110 SUM = SUM + A(IND(L),IND(I))*A(IND(I),IND(L)) A(IND(L),IND(L)) = SUM DO 120 I=L+1,M 120 A(IND(I),IND(L)) = A(IND(L),IND(I)) IMAX = IND(L) DO 130 I=L+1,M IF( IMAX.GT.IND(I) ) THEN IND(I-1) = IND(I) IND(I) = IMAX END IF 130 CONTINUE 200 CONTINUE C RETURN E N D SUBROUTINE ISTATE( M,MJ,XMEAN,XVAR,XF,VF ) C C ... Initial state ... C C Inputs: C M: Dimension of the state vector C MJ: Adjustable dimension of F C XMEAN: Mean value C XVAR: Varaince C Outputs: C XF: State vector, X(0|0) C VF: State covarance matrix, V(0|0) C IMPLICIT REAL*8(A-H,O-Z) DIMENSION XF(MJ), VF(MJ,MJ) C DO 10 I=1,MJ DO 10 J=1,MJ 10 VF(I,J) = 0.0D0 C DO 20 I=1,M 20 XF(I) = XMEAN C DO 30 I=1,M 30 VF(I,I) = XVAR C RETURN E N D SUBROUTINE PRSTAT( M,TAU2,FF,AIC,XSS,VSS,YMEAN,MJ,NFE,NPE, * NMISS,N0,NN ) C C ... Print out Estimates ... C C Inputs: C M: Model order C TAU2: Variance of the model C FF: Log-Likelihood of the mode C AIC: AIC of the model C XSS: Smoothed state C VSS: Covariance matrices C YMEAN: Mean of the original time series C MJ: Adjustable dimension of F and G C NFE: End position of filtering C NPE: End position of prediction C NMISS: Number of missed intervals C N0(I): Start position of I-th missed interval C NN(I): Numbers of missed observations in I-th interval C IMPLICIT REAL*8(A-H,O-Z) CHARACTER TITLE*72 DIMENSION XSS(MJ,*), VSS(MJ,MJ,*), N0(10), NN(10) COMMON /CMDATA/ TITLE C WRITE(6,600) WRITE(6,610) TITLE WRITE(6,620) M, TAU2, FF, AIC WRITE(6,630) DO 10 J=1,NMISS DO 10 I=1,NN(J) II = N0(J)+I-1 XM = XSS(1,II) + YMEAN 10 WRITE(6,640) II, XM, VSS(1,1,II) DO 20 I=NFE+1,NPE XM = XSS(1,I) + YMEAN 20 WRITE(6,640) I, XM, VSS(1,1,I) RETURN C 600 FORMAT( 1H ,'PROGRAM 9.1: PREDICTION AND INTERPOLATION' ) 610 FORMAT( 1H ,'DATA: ',A72 ) 620 FORMAT( 1H ,'M =',I3,3X,'TAU2 =',D12.5,3X,'LOG-L =',F12.4, * 3X,'AIC =',F12.4 ) 630 FORMAT( 1H ,' I',10X,'MEAN',12X,'VARANCE' ) 640 FORMAT( 1H ,I4,3X,F15.6,5X,D13.6 ) E N D SUBROUTINE PTSTAT( Y,M,TAU2,FF,AIC,XSS,VSS,N,NFE,NPE, * MJ,YMEAN,NMISS,N0,NN ) C C ... Plot Original time series, estimated ans errors ... C C Inputs: C Y(I): Time Series C M: Model order C TAU2: Variance of the model C FF: Log-Likelihood of the mode C AIC: AIC of the model C XSS: Smoothed state C VSS: Covariance matrices C N: Data length C NFE: End position of filtering C NPE: End position of prediction C YMEAN: Mean of the original time series C MJ: Adjustable dimension of F and G C NMISS: Number of missed intervals C N0(I): Start position of I-th missed interval C NN(I): Numbers of missed observations in I-th interval C MODIFIED 2/15/93 C IMPLICIT REAL*8(A-H,O-Z) CHARACTER VNAME*8 REAL*4 XX, YY DIMENSION Y(N) DIMENSION XSS(MJ,N), VSS(MJ,MJ,N) DIMENSION DATA(2000), VNAME(6), VALUE(6), N0(10), NN(10) C VNAME(1) = 'M = ' VNAME(2) = 'TAU2 = ' VNAME(3) = 'FF = ' VNAME(4) = 'AIC = ' VALUE(1) = M VALUE(2) = TAU2 VALUE(3) = FF VALUE(4) = AIC CALL PLOTS CALL HEADER( '@TEST.FILTER2: KALMAN FILTER/SMOOTHER ',38,4, * VNAME,VALUE ) WY = 6.0 WX = 20.0 DX = 12.0 FN = NPE IY = 10 CALL MAXMIN( Y,N,YMIN1,YMAX1,DY ) CALL NEWPEN( 2 ) CALL PLOT( 3.0,16.5-SNGL(WY),-3 ) CALL AXISXY( 0.0D0,0.0D0,WX,WY,0.0D0,FN,YMIN1,YMAX1,DX,DY, * 0.2D0,1,IY,2 ) CALL SYMBOL(0.25,SNGL(WY)+0.25,0.25,'TIME SERIES AND ESTIMATED', * 0.0,25 ) C DO 120 J=3,5,1 DO 100 I=1,NPE 100 DATA(I) = XSS(1,I) + DSQRT( VSS(1,1,I) )*(J-4) + YMEAN C CALL NEWPEN( 1 ) IF(J.EQ.4) CALL NEWPEN( 2 ) CALL PLOTY( DATA,NPE,YMIN1,YMAX1,WX,WY,1,1 ) 120 CONTINUE DO 130 J=1,NMISS DO 130 I=1,NN(J) II = N0(J)+I-1 XX = (WX*II)/NPE YY = WY*(Y(II)-YMIN1)/(YMAX1-YMIN1) 130 CALL SYMBOL( XX,YY,0.1,1,0.0,-1 ) DO 140 II=NFE+1,NPE XX = (WX*II)/NPE YY = WY*(Y(II)-YMIN1)/(YMAX1-YMIN1) 140 CALL SYMBOL( XX,YY,0.1,1,0.0,-1 ) C C DO 160 I=1,NPE 160 DATA(I) = Y(I) - XSS(1,I) - YMEAN CALL MAXMIN( DATA,NPE,DMIN,DMAX,DDY ) YMIN3 = DY*INT( DMIN/DY ) YMAX3 = DY*INT( DMAX/DY ) IF( DMIN.LT.YMIN3 ) YMIN3 = YMIN3 - DY IF( DMAX.GT.YMAX3 ) YMAX3 = YMAX3 + DY WY3 = WY*(YMAX3-YMIN3)/(YMAX1-YMIN1) C CALL PLOT( 0.0,-SNGL(WY3+2.0),-3 ) CALL NEWPEN( 2 ) CALL AXISXY( 0.0D0,0.0D0,WX,WY3,0.0D0,FN,YMIN3,YMAX3,DX,DY, * 0.2D0,1,IY,2) CALL SYMBOL( 0.25,WY3+0.25,0.25,'ESTIMATION ERROR',0.0,16 ) C CALL NEWPEN( 1 ) IF(J.EQ.4) CALL NEWPEN( 2 ) CALL PLOTY2( DATA,NPE,YMIN3,YMAX3,WX,WY3,1,1,0.0D0 ) C CALL PLOTE RETURN E N D SUBROUTINE SETFGH( M,MJ,K,AR,TAU2,F,G,H,Q,R ) C C ... Specify the state space representation of AR model ... C C Inputs: C M: Order of AR model C MJ: Adjustable dimension of F C AR(I): AR coefficents C TAU2: Innovation variance C Outputs: C F: State transition matrix C G,H: M-vector C Q: System noise varaince C R: Observational noise varaince (=0) C MODIFIED 2/15/93 C IMPLICIT REAL*8(A-H,O-Z) DIMENSION F(MJ,MJ), G(MJ), H(MJ), AR(MJ), Q(K,K) C DO 10 J=1,MJ G(J) = 0.0D0 H(J) = 0.0D0 DO 10 I=1,MJ 10 F(I,J) = 0.0D0 C DO 20 I=1,M 20 F(1,I) = AR(I) DO 30 I=2,M 30 F(I,I-1) = 1.0D0 G(1) = 1.0D0 H(1) = 1.0D0 Q(1,1)=TAU2 R = 0.0D0 C RETURN E N D SUBROUTINE MOMENT( Y,N,YMEAN,VAR ) C C ... Mean and variance of the data ... C C Inputs: C Y(I): data C N: data length C Outputs: C YMEAN: mean C VAR: variance C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(N) C SUM = 0.0D0 DO 10 I=1,N 10 SUM = SUM + Y(I) YMEAN = SUM/N SUM = 0.0D0 DO 20 I=1,N 20 SUM = SUM + (Y(I)-YMEAN)**2 VAR = SUM/N C RETURN E N D SUBROUTINE AXISXY(X,Y,WX,WY,X0,X1,Y0,Y1,DX,DY,DWC,IBOX,IX,IY) C C ... This subroutine draws X and Y axes. C C Inputs: C X,Y: location of the left bottom of the figure C WX,WY: width and height of the figure C X0,X1: lower and upper bounds of the X axis C Y0,Y1: lower and upper bounds of the Y axis C DX,DY: step width in X and Y axes C DWC: size of characters C IBOX: = 0 draw X and Y axes only C = 1 draw window doundary C IX,IY: number of additional click in each step C Modified: 8/31/90, 11/21/90 C REAL*8 X,Y,WX,WY,X0,X1,Y0,Y1,DX,DY,DWC C FX = X FY = Y FWX = WX FWY = WY FWC = DWC CALL NEWPEN( 2 ) CALL PLOT( FX,FY,-3 ) CALL PLOT( 0.0,FWY,3 ) CALL PLOT( 0.0,0.0,2 ) CALL PLOT( FWX,0.0,2 ) IF ( IBOX.EQ.1 ) THEN CALL PLOT( FWX,FWY,2 ) CALL PLOT( 0.0,FWY,2 ) END IF C C ... draw X axis ... C WC = DWC NX = (X1-X0)/DX + 0.001 NY = (Y1-Y0)/DY + 0.001 MX = -DLOG10( DX ) + 0.001 IF(DX.LT.1.0D0) MX = -DLOG10( DX ) + 0.999 IF( MX.LE.0 ) MX = -1 D = WX*DX/(X1-X0) DO 10 I=1,NX+1 XN = X0 + DX*(I-1) XX = D*(I-1) CALL PLOT( XX,-FWC/2,3 ) CALL PLOT( XX,0.0,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( XX,FWY-FWC/2,3 ) CALL PLOT( XX,FWY,2 ) END IF LX = 0 IF(XN.GT.0.0D0) LX = ALOG10( XN ) + 0.001 IF( MX.GT.0 ) XXX = WC*(LX+MX)*0.5 + WC*4/7.0 IF( MX.LT.0 ) XXX = WC*LX*0.5 + WC*4/7.0 CALL NUMBER( XX-XXX,-(1.5*FWC+0.1),FWC,XN,0.0,MX ) 10 CONTINUE IF( IX.GT.1 ) THEN CALL NEWPEN(1) DO 20 I=1,NX+1 DO 20 J=1,IX XX = D*(I-1) + (D*J)/IX IF( XX.LE.WX ) THEN CALL PLOT( XX,FWC/2,3 ) CALL PLOT( XX,0.0,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( XX,FWY-FWC/2,3 ) CALL PLOT( XX,FWY,2 ) END IF END IF 20 CONTINUE CALL NEWPEN( 2 ) END IF C C ... draw Y axis ... C D = WY*DY/(Y1-Y0) MY = -DLOG10( DY ) + 0.99999 IF( MY.LE.0 ) MY = -1 DO 30 I=1,NY+1 YN = Y0 + DY*(I-1) YY = D*(I-1) CALL PLOT( -FWC/2,YY,3 ) CALL PLOT( 0.0,YY,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( FWX-FWC/2,YY,3 ) CALL PLOT( FWX,YY,2 ) END IF LY = 1 IF(YN.LT. 0.0) LY = 2 IF(YN.GE. 1.0) LY = ALOG10( YN )+1.00001 IF(YN.LE.-1.0) LY = ALOG10(-YN )+2.00001 IF( MY.GT.0 ) YYY = WC*(LY+MY+1) + 0.20 IF( MY.LT.0 ) YYY = WC*LY + 0.20 CALL NUMBER( -YYY,YY-FWC/2,FWC,YN,0.0,MY ) 30 CONTINUE IF( IX.GT.1 ) THEN CALL NEWPEN(1) DO 40 I=1,NY+1 DO 40 J=1,IY YY = D*(I-1) + (D*J)/IY IF( YY.LE.WY ) THEN CALL PLOT( FWC/2,YY,3 ) CALL PLOT( 0.0,YY,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( FWX-FWC/2,YY,3 ) CALL PLOT( FWX,YY,2 ) END IF END IF 40 CONTINUE CALL NEWPEN( 2 ) END IF C RETURN E N D SUBROUTINE HEADER( PNAME,NP,M,VNAME,VALUE ) C C ... This subroutine draw title of data, program name, names C and values of parameters ... C C Inputs: C PNAME: program name C NP: number of characters of program name C M: number of parameters C VNAME(I): name of the I-th parameter C VALUE(I): value of the I-th parameter C Input via common: C TITLE: title of the data set C CHARACTER PNAME*80, TITLE*72 CHARACTER*8 VNAME(M) REAL*8 VALUE DIMENSION VALUE(M), DAY(2), TIME(3) COMMON /CMDATA/ TITLE C CALL NEWPEN( 2 ) CALL DATE( DAY ) CALL CLOCK( TIME,1 ) C CALL SYMBOL( 2.0,18.7,0.25,TITLE,0.0,72 ) CALL SYMBOL( 2.0,18.3,0.20,PNAME,0.0,NP ) CALL NEWPEN( 1 ) CALL SYMBOL( 2.0,17.9,0.20,DAY,0.0,8 ) CALL SYMBOL( 4.0,17.9,0.20,TIME,0.0,12 ) YY = 17.9 DO 10 I=1,M IF( MOD(I,4).EQ.1 ) THEN YY = YY - 0.40 X1 = 2.0 X2 = 3.5 ELSE X1 = X1 + 5.0 X2 = X2 + 5.0 END IF CALL SYMBOL( X1,YY,0.20,VNAME(I),0.0,8 ) 10 CALL NUMBER( X2,YY,0.20,SNGL(VALUE(I)),0.0,6 ) C RETURN E N D SUBROUTINE MAXMIN( X,N,XMIN0,XMAX0,DXL ) C C ... This subroutine determines the minimum, the maximum and C the step width ... C C Inputs: C X(I): data C N: data length C Outputs: C XMIN0: the minimum bound for the figure C XMAX0: the maximum bound for the figure C DXL: step width C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(N) C XMIN = 1.0D30 XMAX =-1.0D30 C DO 10 I=1,N IF( X(I) .LT. XMIN ) XMIN = X(I) 10 IF( X(I) .GT. XMAX ) XMAX = X(I) DX = XMAX-XMIN IF( DLOG10(DX) .GE. 0.0D0 ) DXL = INT( DLOG10(DX) ) IF( DLOG10(DX) .LT. 0.0D0 ) DXL = INT( DLOG10(DX) )-1.0 DXL = 10.0**DXL IF( DX/DXL.GT.6.0D0 ) DXL = DXL*2.0 DIF = INT( DX/DXL ) XMIN0 = INT( XMIN/DXL )*DXL XMAX0 = XMIN0 + DIF*DXL IF( XMIN0 .GT. XMIN ) XMIN0 = XMIN0 - DXL 30 IF( XMAX0 .GE. XMAX ) GO TO 40 XMAX0 = XMAX0 + DXL GO TO 30 40 CONTINUE C RETURN E N D SUBROUTINE PLOTY( Y,N,YMIN,YMAX,WX,WY,IPOS,ISW ) C C ... This subroutine draws time vs data plot of the data ... C C Inputs: C Y(I): data C N: data length C YMIN: minimum bound for the Y axis C YMAX: maximum bound for the Y axis C WX,WY: width and height of the figure C IPOS: starting position of the data C ISW: = 1 connect data by straight line C = 2 connect data by dashed line C = 3 express data by character (character code=ISW) C Modified: 8/31/90 C REAL*8 Y(N), YMIN, YMAX, WX, WY C DX = WX/(N-1+IPOS) DY = WY/(YMAX - YMIN) C CALL NEWPEN( 1 ) DO 100 I=1,N XX = DX*(I-1+IPOS) YY = (Y(I) - YMIN)*DY IF( I.EQ.1 .AND. ISW.LT.3 ) THEN CALL PLOT( XX,YY,3 ) ELSE IF(ISW.EQ.1) CALL PLOT( XX,YY,2 ) IF(ISW.EQ.2) CALL DASHP( XX,YY,0.2 ) IF(ISW.GE.3) CALL SYMBOL( XX,YY,0.2,ISW,0.0,-1 ) END IF 100 CONTINUE C RETURN E N D SUBROUTINE PLOTY2( Y,N,YMIN,YMAX,WX,WY,IOFF,ISW,YLEVEL ) C C ... This subroutine draws box car graph of the data ... C C Inputs: C Y(I): data C N: data length C YMIN: minimum bound for the Y axis C YMAX: maximum bound for the Y axis C WX,WY: width and height of the figure C IPOS: starting position of the data C ISW: = 1 draw box with thick line C = 2 draw box with thin line C YLEVEL: ground level for the box car C REAL*8 Y(N), YMIN, YMAX, WX, WY, YLEVEL C DX = WX/(N-1+IOFF) DY = WY/(YMAX - YMIN) C CALL NEWPEN( 3 ) IF( N.GT.75 ) CALL NEWPEN( 2 ) IF( N.GT.150) CALL NEWPEN( 1 ) IF( ISW.EQ.2 ) CALL NEWPEN( 1 ) DO 100 I=1,N XX = DX*(I-1+IOFF) YY = (Y(I) - YMIN)*DY Y0 = (YLEVEL-YMIN)*DY CALL PLOT( XX,Y0,3 ) CALL PLOT( XX,YY,2 ) 100 CONTINUE CALL NEWPEN( 2 ) C RETURN E N D SUBROUTINE READTS( IDEV,Y,N ) REAL*8 Y(N) CHARACTER*72 TITLE COMMON /CMDATA/ TITLE C C OPEN( IDEV,FILE='temp.dat' ) READ( IDEV,1 ) TITLE READ( IDEV,* ) N READ( IDEV,* ) (Y(I),I=1,N) C CLOSE( IDEV ) RETURN 1 FORMAT( A72 ) E N D