C PROGRAM 2.2 PTSCAT C C ... This program draws histograms and scatter plot of C multivariate time series ... C C The following inputs are required in the subroutine READMD. C TITLE: title of data C N: data length C L: number of variates C IFM: written format C = 1 ((Y(I,J),J=1,L),I=1,N) C = 2 ((Y(I,J),I=1,N),J=1,L) C Y(I,J): multivariate time series C Parameters: C MJ: adjustable dimension of Y (MJ.GE.N) C MJ1: adjustable dimension of Y (MJ1.GE.L) C MAR.24,1989: modified 12/20/90 C PARAMETER( MJ=1000,MJ1=7) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(MJ,MJ1), OUTMIN(10), OUTMAX(10) DATA OUTMIN /10*-1.0D30/, OUTMAX /10*1.0D30/ C C ... read in multivariate time series ... C CALL READMD( 1,MJ,Y,N,L ) C C ... draw histograms and scatter-plot ... C CALL PTSCAT( Y,L,MJ,N,OUTMIN,OUTMAX ) C STOP E N D SUBROUTINE HISTGR( Y,N,OUTMIN,OUTMAX,HIST,K,HMAX ) C C ... This subroutine counts number of observations for each bin ... C C Inputs: C Y(I): data C N: number of observations C OUTMIN: bound for outliers in low side C OUTMAX: bound for outliers in high side C Outputs: C HIST(I): number of observations in I-th bin C K: number of bins C HMAX: the maximum of HIST(I) C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(N), HIST(40) C CALL RANGE( Y,N,OUTMIN,OUTMAX,YMIN,YMAX ) K = 20 IF( N.LT.500 ) K = 15 IF( N.LT.200 ) K = 10 40 CONTINUE DO 10 I=1,K 10 HIST(I) = 0 DO 20 I=1,N J = (Y(I)-YMIN)/( YMAX-YMIN )*K + 1 20 IF( J.GE.1 .AND. J.LE.K ) HIST(J) = HIST(J) + 1 HMAX = -1.0D30 DO 50 J=1,K 50 IF( HIST(J).GT.HMAX ) HMAX = HIST(J) RETURN E N D SUBROUTINE PLOT4( Y,N,YMIN,YMAX,WX,WY ) C C ... This subroutine draws histogram ... C C Inputs: C Y(I): number of counts at I-th bin C N: number of bins C YMIN: minimum in Y axis C YMAX: maximum in Y axis C WX: width of the figure C WY: hight of the figure C REAL*8 Y, YMIN, YMAX, WX, WY DIMENSION Y(N) C DX = WX/N DY = WY/(YMAX-YMIN) C CALL PLOT( 0.0,0.0,-3 ) DO 100 I=1,N XX = DX*(I-1) XX1 = DX*I YY = (Y(I) - YMIN)*DY CALL PLOT( XX,YY,2 ) CALL PLOT( XX1,YY,2 ) 100 CONTINUE C RETURN E N D SUBROUTINE PLOTSC( X,Y,N,XMIN,XMAX,YMIN,YMAX,WX,WY,SIZE,CH ) C C ... This subroutine draws scatter plot ... C C Inputs: C X(I): data for x-axis C Y(I): data for y-axis C N: data length C XMIN,XMAX: edges of x-axis C YMIN,YMAX: edges of y-axis C WX,WY: width andhight of the figure C SIZE: size of the symbol C CH: character code C REAL*8 X, Y, XMIN, XMAX, YMIN, YMAX, WX, WY, SIZE CHARACTER CH*1 DIMENSION X(N), Y(N) C DX = WX/(XMAX-XMIN) DY = WY/(YMAX-YMIN) C DO 10 I=1,N XX = (X(I)-XMIN)*DX-SIZE/2 YY = (Y(I)-YMIN)*DY-SIZE/2 IF( XX.GE.0.0.AND.XX.LE.WX ) THEN IF( YY.GE.0.0.AND.YY.LE.WY ) THEN CALL SYMBOL( XX,YY,SIZE,CH,0.0,1 ) END IF END IF 10 CONTINUE C RETURN E N D SUBROUTINE PLTBOX ( WX,WY,XMIN,XMAX,XSTEP,YMIN,YMAX,YSTEP,IPEN ) C C ... This subroutine draws box window ... C C Inputs: C WX,WY: width and height of the box window C XMIN,XMAX: C XSTEP: C YMIN,YMAX: C YSTEP: C IPEN: thickness or colour of the window C REAL*8 WX, WY, XMIN, XMAX, XSTEP, YMIN, YMAX, YSTEP C CALL NEWPEN( IPEN ) CALL PLOT( 0.0,0.0,3 ) CALL PLOT( SNGL(WX),0.0,2 ) CALL PLOT( SNGL(WX),SNGL(WY),2 ) CALL PLOT( 0.0,SNGL(WY),2 ) CALL PLOT( 0.0,0.0,2 ) C CALL NEWPEN(1) DX = XSTEP/(XMAX-XMIN) NX = 1.0D0/DX DX = DX*WX NX = ( XMAX-XMIN )/XSTEP DO 10 I=1,NX CALL PLOT( DX*I,0.0,3 ) 10 CALL PLOT( DX*I,0.1,2 ) DY = YSTEP/( YMAX-YMIN ) NY = 1.0D0/DY DY = DY*WY DO 20 I=1,NY CALL PLOT( 0.0,DY*I,3 ) 20 CALL PLOT( 0.1,DY*I,2 ) RETURN E N D SUBROUTINE PTSCAT( Y,ID,MJ,N,OUTMIN,OUTMAX ) C C ... This subroutin plots cross-correlation function ... C C Inputs: C Y(I,J): multivariate time series C ID: number of variates C MJ: adjustable dimension of Y (MJ.GE.N) C N: data length C OUTMIN: C OUTMAX: C IMPLICIT REAL*8(A-H,O-Z) CHARACTER*8 VNAME REAL*4 SWX DIMENSION Y(MJ,ID), HIST(20), OUTMIN(*), OUTMAX(*) DIMENSION VNAME(3), VALUE(3) C WX = 16.0/ID -0.5 SWX = WX VNAME(1) = 'N =' VNAME(2) = 'ID =' VALUE(1) = N VALUE(2) = ID C CALL PLOTS C call plots( 1,0,0,1,0 ) C call form( 1 ) C call factor( 10.0 ) CALL HEADER( 'PROGRAM 2.2: HISTOGRAM AND SCATTER PLOT', * 40,2,VNAME,VALUE ) CALL PLOT( 2.0,16.5-SWX,-3 ) C DO 10 I=1,ID IF( I.GE.2 ) CALL PLOT( -(SWX+0.5)*(ID-I+1),-SWX-0.5,-3 ) CALL RANGE ( Y(1,I),N,OUTMIN(I),OUTMAX(I),YMIN,YMAX ) DO 10 J=I,ID IF( J.GE.2 ) CALL PLOT( SWX+0.5,0.0,-3 ) IPEN = 1 IF( I.EQ.J ) THEN CALL HISTGR( Y(1,I),N,OUTMIN(I),OUTMAX(I),HIST,NH,HMAX ) CALL PLTBOX( WX,WX,0.0D0,DBLE(NH),10.0D0,-1.0D0,1.0D0,1.0D0,2 ) CALL PLOT4( HIST,NH,0.0D0,HMAX,WX,WX ) ELSE CALL PLTBOX( WX,WX,0.0D0,DBLE(NH),10.D0,-1.D0,1.D0,1.D0,IPEN ) CALL RANGE ( Y(1,J),N,OUTMIN(J),OUTMAX(J),XMIN,XMAX ) CALL PLOTSC( Y(1,J),Y(1,I),N,XMIN,XMAX,YMIN,YMAX,WX,WX,0.1D0, * '+' ) END IF 10 CONTINUE CALL PLOTE C call plot( 0.0,0.0,999 ) C RETURN E N D SUBROUTINE RANGE ( Y,N,OUTMIN,OUTMAX,YMIN,YMAX ) C C ... This subroutine finds the minimum and the maximuM of data .. C C Inputs: C Y(I): data C N: data length C OUTMIN: bound for outliers in low side C OUTMAX: bound for outliers in high side C Outputs: C YMIN: minimum of Y(I), I=1,...,N C YMAX: maximum of Y(I), I=1,...,N C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(N) C YMIN = 1.0D30 YMAX =-1.0D30 C DO 10 I=1,N IF( Y(I).LT.YMIN .AND. Y(I).GT.OUTMIN ) YMIN = Y(I) 10 IF( Y(I).GT.YMAX .AND. Y(I).LT.OUTMAX ) YMAX = Y(I) C RETURN E N D SUBROUTINE READMD( IDEV,NMAX,Y,N,L ) C C ... THIS SUBROUTINE READS IN MULTI-VARIATE TIME SERIES ... C C INPUTS: C IDEV: INPUT FILE SPECIFICATION C NMAX: ADJUSTABLE DIMENSION OF Y(I,J) C OUTPUTS: C Y(I,J): TIME SERIES (I=1,N;J=1,L) C N: DATA LENGTH C L: DIMENSION OF THE TIME SERIES C CHARACTER*72 TITLE REAL*8 Y(NMAX,*) COMMON /CMDATA/ TITLE C C OPEN( IDEV,FILE='temp2.dat' ) READ( IDEV,1 ) TITLE READ( IDEV,* ) N, L, IFM IF( IFM.EQ.0 ) READ( IDEV,* ) ((Y(I,J),J=1,L),I=1,N) IF( IFM.GE.1 ) READ( IDEV,* ) ((Y(I,J),I=1,N),J=1,L) C CLOSE( IDEV ) C RETURN 1 FORMAT( A72 ) 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 = 4.0 ELSE X1 = X1 + 5.5 X2 = X2 + 5.5 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