C PROGRAM 14.1 NGSMTH C C ... NON-GAUSSIAN SMOOTHING ... C C INPUTS: C NOISEV: TYPE OF SYSTEM NOISE DENSITY (0,1,2,3) C NOISEW: TYPE OF OBSER. NOISE DENSITY (0,1,2,3,4) C TAU2: VARIANCE OR DISPERSION OF SYSTEM NOISE C SIG2: VARIANCE OR DISPERSION OF OBSERVATION NOISE C BV: SHAPE PARAMETER OF SYSTEM NOISE (FOR NOISEV = 2) C BW: SHAPE PARAMETER OF OBSER. NOISE (FOR NOISEW = 2) C PARAMETERS: C MJ: ADJUSTABLE DIMENSION OF Y, PS, FS, SS (MJ >= N) C K: NUMBER OF INTERVALS + 1 C @NFILTER.F: 5/14/85, 6/11/87, 5/19/88, 2/14/91 C MODIFIED 2/16/93 C PARAMETER( MJ=500, K=201 ) IMPLICIT REAL*8(A-H,O-Z) CHARACTER TITLE*72 DIMENSION Y(MJ), P(K), F(K), S(K), T(K), Q(-K:K) DIMENSION YM(10), ST(7), TREND(MJ,7), LOC(MJ) REAL*4 PS(K,MJ), FS(K,MJ), SS(K,MJ) COMMON /COMAAA/ Y COMMON /C91214/ XMIN, XMAX, FIGMIN, FIGMAX, YM COMMON /C91215/ NOISEV, NOISEW, INITD, ITF, ITH COMMON /C91216/ B, OUTMIN, OUTMAX COMMON /C91218/ SIG2, TAU2, BW, BV COMMON /C91219/ N, KK, NS, NFE, NPE COMMON /CMDATA/ TITLE COMMON / DDD / FF , AIC , SD C NAMELIST /PARAM/ XMIN, XMAX, TAU20, SIG20, BC, IOPT, B, C * OUTMIN, OUTMAX, NS, NFE, NPE, FIGMIN, FIGMAX, C * NOISEV, NOISEW, ITF, ITH C EXTERNAL GAUSS EXTERNAL PEARSN EXTERNAL TWOEXP EXTERNAL USERV C C ... READ DATA ... C CALL READTS( 1,Y,N ) CALL MOMENT( Y,N,YM(1),YM(2) ) C LOC(1) = 0 KK = K CALL DEFALT READ( 5,* ) NOISEV, TAU2, BV READ( 5,* ) NOISEW, SIG2, BW C DX = (XMAX-XMIN)/(K-1) CALL IDIST( F,K,YM(1),YM(2),XMIN,DX ) CALL NORMLZ( F,K,DX,SSUM ) IF( NOISEV.EQ.0 ) CALL TRANS( USERV ,K,DX,TAU2,BV,Q ) IF( NOISEV.EQ.1 ) CALL TRANS( GAUSS ,K,DX,TAU2,BV,Q ) IF( NOISEV.EQ.2 ) CALL TRANS( PEARSN,K,DX,TAU2,BV,Q ) IF( NOISEV.EQ.3 ) CALL TRANS( TWOEXP,K,DX,TAU2,BV,Q ) CALL NGSMTH( Y,P,F,S,T,N,K,DX,XMIN,Q,FF,PS,FS,SS,LOC ) C DO 30 I=1,NPE DO 10 J=1,K 10 S(J) = SS(J,I) CALL PINTVL( S,K,XMIN,DX,ST ) DO 20 J=1,7 20 TREND(I,J) = ST(J) + DX*LOC(I) 30 CONTINUE C CALL PRNGSM( NOISEV,NOISEW,TAU2,SIG2,FF,TREND,MJ,N ) C CALL PTNGSM( N,NPE,NOISEV,NOISEW,TAU2,SIG2,FF,B,Y,TREND,MJ, * FIGMIN,FIGMAX,SS,K,LOC ) STOP E N D SUBROUTINE DEFALT C C ... THIS SUBROUTINE SETS DEFAULT VALUES OF PARAMETERS ... C IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(500) COMMON /COMAAA/ Y COMMON /C91214/ XMIN, XMAX, FIGMIN, FIGMAX, YM(10) COMMON /C91215/ NOISEV, NOISEW, INITD, ITF, ITH COMMON /C91216/ B, OUTMIN, OUTMAX COMMON /C91219/ N, KK, NS, NFE, NPE NOISEV = 2 NOISEW = 1 INITD = 1 ITF = 1 ITH = 1 B = 1.00D0 OUTMIN = -1.0D30 OUTMAX = 1.0D30 CALL MAXMIN( Y,N,XMIN,XMAX,DY ) NS = 1 NFE = N NPE = N FIGMIN = XMIN FIGMAX = XMAX C RETURN E N D SUBROUTINE NGSMTH( Y,P,F,S,T,N,K,DX,XMIN,Q,FF,PS,FS,SS,LOC ) C C ... NON-GAUSSIAN SMOOTHER ... C C INPUTS: C Y(I): TIME SERIES C P(I): INITIAL DENSITY C N: DATA LENGTH C K: NUMBER OF INTERVALS IN STEP FUNCTION APPROXIMATION C DX: WIDTH OF INTERVAL C XMIN: MINIMUM OF THE INTERVAL C Q: SYSTEM NOISE DENSITY C OUTPUTS: C FF: LOG-LIKELIHOOD C LOC(I): LOCATION OF THE CENTER OF THE INTERVAL AT STEP I C SS: SMOOTHED DENSITY C IMPLICIT REAL*8(A-H,O-Z) DIMENSION P(K), F(K), S(K), T(K), Y(N), Q(-K:K), LOC(N) REAL*4 FS(K,N), PS(K,N), SS(K,N) COMMON /C91215/ NOISEV, NOISEW, INITD, ITF, ITH COMMON /C91216/ B, OUTMIN, OUTMAX COMMON /C91219/ NN, KK, NS, NFE, NPE C FF = 0.0D0 C DO 200 II=NS,NPE C C ... CONVOLUTION (SYSTEM NOISE) ... C CALL CONVOL( Q,F,K,P ) CALL NORMLZ( P,K,DX,PSUM ) C C ... BAYES FORMULA ... C IF( Y(II).LE.OUTMIN .OR. Y(II).GE.OUTMAX .OR. II.GT.NFE ) THEN DO 110 I=1,K 110 F(I) = P(I) ELSE CALL BAYES( P,K,XMIN,DX,Y(II),F,LOC(II) ) CALL NORMLZ( F,K,DX,FINT ) C C ... LIKELIHOOD COMPUTATION ... C FF = FF + DLOG( FINT ) C IF( MOD(II,10).EQ.0) WRITE(6,*) II,FF END IF C C ... SAVE FOR SMOOTHING ... C DO 130 I=1,K PS(I,II) = P(I) 130 FS(I,II) = F(I) C C ... SHIFT ORIGIN ... C CALL SHIFT( F,K,T,II,N,LOC ) C 200 CONTINUE C C ... SMOOTHING ... C DO 190 J=NFE,NPE DO 190 I=1,K 190 SS(I,J) = FS(I,J) DO 195 I=1,K 195 S(I) = FS(I,NFE) C DO 300 II=NFE-1,NS,-1 C IF(MOD(II,10).EQ.0) WRITE(6,*) II DO 210 I=1,K T(I) = 0.0D0 P(I) = 0.0D0 210 F(I) = FS(I,II) DO 220 I=1,K J = I - (LOC(II+1)-LOC(II)) IF( J.GE.1.AND.J.LE.K ) P(I) = PS(J,II+1) 220 IF( J.GE.1.AND.J.LE.K ) T(I) = S(J) DO 230 I=1,K 230 S(I) = T(I) C CALL SCONVL( Q,S,P,F,K,T ) CALL NORMLZ( T,K,DX,TSUM ) C DO 240 I=1,K S(I) = T(I) 240 SS(I,II) = S(I) 300 CONTINUE C RETURN E N D SUBROUTINE PINTVL( P,K,XMIN,DX,Y ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION P(K), Y(7), PROB(7), P1(401) DATA PROB /0.0013D0, 0.0227D0, 0.1587D0, 0.5000D0, 0.8413D0, * 0.9773D0,0.9987D0/ C P1(1) = 0.0 DO 10 I=2,K 10 P1(I) = P1(I-1) + (P(I-1) + P(I))*DX/2.0 C DO 30 J=1,7 PP = PROB(J) DO 20 I=2,K IF(P1(I-1).LE.PP .AND. P1(I).GT.PP) GO TO 30 20 CONTINUE 30 Y(J) = XMIN + (I-2)*DX + DX*(PP - P1(I-1))/(P1(I) - P1(I-1)) C RETURN E N D SUBROUTINE CONVOL( Q,S,K,P ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION S(K), P(K), Q(-K:K) C DO 20 I=1,K J1 = 1-I J2 = K-I SUM = 0.0 DO 10 J=J1,J2 10 SUM = SUM + S(I+J)*Q(J) 20 P(I) = SUM C RETURN E N D SUBROUTINE SCONVL( Q,P,R,S,K,T ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION S(K), P(K), R(K), T(K), Q(-K:K) C DO 20 I=1,K J1 = 1-I J2 = K-I SUM = 0.0D0 DO 10 J=J1,J2 10 IF(P(I+J).GT.0.0D0) SUM = SUM + P(I+J)/R(I+J)*Q(J) 20 T(I) = S(I)*SUM RETURN E N D SUBROUTINE TRANS( FUNCT,K,DX,TAU2,BV,Q ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Q(-K:K), PARAM(3) EXTERNAL FUNCT C PARAM(1) = 0.0D0 PARAM(2) = TAU2 PARAM(3) = BV C DO 20 I=1-K,K-1 X0 = -DX*I - DX/2 SUM = (FUNCT(X0,PARAM) + FUNCT(X0+DX,PARAM))/2 DO 10 J=1,49 X = X0 + (DX*J)/50 10 SUM = SUM + FUNCT( X,PARAM ) 20 Q(I) = SUM*DX/50 C RETURN E N D SUBROUTINE NORMLZ( P,K,DX,SUM ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION P(K) C SUM = 0.0D0 DO 10 I=1,K 10 SUM = SUM + P(I) SUM = SUM*DX DO 30 I=1,K 30 P(I) = P(I)/SUM C RETURN E N D SUBROUTINE IDIST( P,K,P1,P2,XMIN,DX ) IMPLICIT REAL*8( A-H,O-Z ) DIMENSION P(K), PARAM(3) COMMON /C91215/ NOISEV, NOISEW, INITD, ITF, ITH C PARAM(1) = P1 PARAM(2) = P2 C DO 10 I=1,K X = XMIN + DX*(I-1) IF( INITD.EQ.0 ) P(I) = USERI( X,PARAM ) IF( INITD.EQ.1 ) P(I) = GAUSS( X,PARAM ) IF( INITD.EQ.2 ) P(I) = UNIF ( X,PARAM ) 10 CONTINUE RETURN END SUBROUTINE BAYES( P,K,XMIN,DX,Y,F,LSHIFT ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION P(K), F(K), PARAM(3) COMMON /C91215/ NOISEV, NOISEW, INITD, ITF, ITH COMMON /C91218/ SIG2, TAU2, BW, BV EXTERNAL USERW EXTERNAL GAUSS EXTERNAL PEARSN EXTERNAL TWOEXP EXTERNAL DBLEXP C PARAM(2) = SIG2 PARAM(3) = BW C DO 10 I=1,K PARAM(1) = XMIN + DX*(I-1+LSHIFT) IF( NOISEW.EQ.0 ) F(I) = P(I)*USERW ( Y,PARAM ) IF( NOISEW.EQ.1 ) F(I) = P(I)*GAUSS ( Y,PARAM ) IF( NOISEW.EQ.2 ) F(I) = P(I)*PEARSN( Y,PARAM ) IF( NOISEW.EQ.3 ) F(I) = P(I)*TWOEXP( Y,PARAM ) IF( NOISEW.EQ.4 ) F(I) = P(I)*DBLEXP( Y,PARAM ) 10 CONTINUE RETURN E N D DOUBLE PRECISION FUNCTION USERW( Y,PARAM ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION PARAM(3) DATA C1 /2.506628275D0/ C YMEAN = 0.0D0 VAR = PARAM(1) USERW = DEXP( -(Y-YMEAN)**2/(2*VAR) )/(C1*DSQRT( VAR )) RETURN E N D DOUBLE PRECISION FUNCTION USERV( X,PARAM ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION PARAM(3) DATA C1 /2.506628275D0/ C USERV = DEXP( -(X-PARAM(1))**2/(2*PARAM(2)) ) * /(C1*DSQRT( PARAM(2) )) RETURN E N D DOUBLE PRECISION FUNCTION TWOEXP( X,PARAM ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION PARAM(2) C TWOEXP = DEXP( -DABS(X-PARAM(1))*PARAM(2) )*PARAM(2)/2.0D0 RETURN C E N D SUBROUTINE SHIFT( F,K,T,II,N,LOC ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION F(K), T(K), LOC(*) C C ... FIND THE POSTERIOR MODE AND SHIFT ORIGIN ... C PMAX = 0.0D0 DO 10 I=1,K IF( F(I).LE.PMAX ) GO TO 10 PMAX = F(I) IMAX = I 10 CONTINUE IF(II.LT.N) LOC(II+1) = LOC(II) + IMAX - (K+1)/2 DO 20 I=1,K J = I + IMAX - (K+1)/2 T(I) = 0.0D0 20 IF(J.GE.1.AND.J.LE.K ) T(I) = F(J) DO 30 I=1,K 30 F(I) = T(I) C RETURN E N D SUBROUTINE POST3D( F,LOC,K,N ) IMPLICIT REAL*8(A-H,O-Z) REAL*4 F(K,N) DIMENSION YH(2000), FF(801), LOC(N) COMMON /C91214/ X0, X1, F0, F1, YM(10) ANGLE = 70.0 C WX = 8.0 WY = 2.0 WZ =12.0 CALL DELX( X0,X1,DX ) Z0 = 0.0 Z1 = N KK = 1 NDIF = 1 IF( N.GE.100 ) NDIF = 2 IF( N.GE.200 ) NDIF = 4 IF( N.GE.300 ) NDIF = 6 IF( N.GE.500 ) NDIF = N/50 DZ = NDIF*10 NN = N/NDIF N0 = NDIF/2 + 1 DO 10 I=1,K 10 FF(I) = F(I,1) CALL MAXMIN( FF,K,Y0,Y1,DY ) C CALL PLOT3A( K,NN,WX,WY,WZ,ANGLE,X0,X1,DX,Y0,Y1,DY,Z0,Z1,DZ, * KN,KK,YH ) DO 100 J=N0,N,NDIF DO 80 I=1,K 80 FF(I) = 0.0D0 II = LOC(J) I1 = MAX0( 1,II ) I2 = MIN0( K,K+II ) DO 90 I=I1,I2 90 FF(I+II) = F(I,J) CALL PLOT3B( FF,K,NN,WX,WY,WZ,ANGLE,Y0,Y1,KK,YH,ZS,ZC ) 100 CONTINUE CALL PLOT( -SNGL(ZC),-SNGL(ZS),-3 ) C RETURN E N D SUBROUTINE PTNGSM( N,NPE,NOISEV,NOISEW,TAU2,SIG2,FF,B,Y,TREND, * MJ,FIGMIN,FIGMAX,SS,K,LOC ) IMPLICIT REAL*8(A-H,O-Z) CHARACTER VNAME*8 REAL*4 SS(K,MJ) DIMENSION Y(N), TREND(MJ,7), VNAME(10), VALUE(10) DIMENSION LOC(N) C CALL PLOTS C call plots( 1,0,0,1,0 ) C call form( 1 ) C call factor( 10.0 ) WX = 10.0 WY = 8.0 VNAME(1) = 'N ' VNAME(2) = 'NOISEV ' VNAME(3) = 'NOISEW ' VNAME(4) = 'TAU2 ' VNAME(5) = 'SIG2 ' VNAME(6) = 'LOG-LK ' VNAME(7) = 'B ' VALUE(1) = N VALUE(2) = NOISEV VALUE(3) = NOISEW VALUE(4) = TAU2 VALUE(5) = SIG2 VALUE(6) = FF VALUE(7) = B CALL HEADER( 'PROGRAM 14.1: NON-GAUSSIAN SMOOTHING' * ,37,7,VNAME,VALUE ) DO 40 I=N+1,NPE 40 Y(I) = Y(N) C CALL DELX( 0.0D0,DBLE(NPE),DDX ) CALL DELX( FIGMIN,FIGMAX,DY ) CALL AXISXY( 2.0D0,16.0D0-WY,WX,WY,0.D0,DFLOAT(NPE),FIGMIN, * FIGMAX,DDX,DY,0.2D0,1,10,2) CALL NEWPEN( 1 ) CALL PLOTY( Y,NPE,FIGMIN,FIGMAX,WX,WY,1,1 ) C CALL AXISXY( WX+2.0,0.0D0,WX,WY,0.0D0,DFLOAT(NPE),FIGMIN, * FIGMAX,DDX,DY,0.2D0,1,10,2) DO 30 J=1,7 CALL NEWPEN( 1 ) IF( J.EQ.4 ) CALL NEWPEN(2) 30 CALL PLOTY( TREND(1,J),NPE,FIGMIN,FIGMAX,WX,WY,1,1 ) C CALL PLOTI C call plot( 0.0,0.0,777 ) CALL HEADER( 'PROGRAM 14.1: ',13,7,VNAME,VALUE ) CALL PLOT( 5.0,2.0,-3 ) CALL POST3D( SS,LOC,K,NPE ) CALL PLOTE C call plot( 0.0,0.0,999 ) C RETURN E N D SUBROUTINE PRNGSM( NOISEV,NOISEW,TAU2,SIG2,FF,TREND,MJ,N ) IMPLICIT REAL*8(A-H,O-Z) CHARACTER TITLE*72 DIMENSION TREND(MJ,7) COMMON /CMDATA/ TITLE C WRITE(6,600) WRITE(6,610) TITLE WRITE(6,620) NOISEV, NOISEW WRITE(6,630) TAU2, SIG2, FF WRITE(6,640) WRITE(6,650) (TREND(I,4),I=1,N) RETURN 600 FORMAT( 1H1,'PROGRAM 14.1: NON-GAUSSIAN SMOOTHING' ) 610 FORMAT( 1H ,A72 ) 620 FORMAT( 1H ,'NOISEV =',I2,3X,'NOISEW =',I2 ) 630 FORMAT( 1H ,'TAU1 =',D12.5,3X,'SIG2 =',D12.4,3X,'LOG-LIKELIHOOD =' * ,F15.5 ) 640 FORMAT( 1H ,'*** ESTIMATED TREND ***' ) 650 FORMAT( 1H ,5D15.6 ) 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 DELX( XMIN,XMAX,DX ) C C ... This subroutine determines step width in drawing axis ... C C Inputs: C XMIN: the minimum value in the axis C XMAX: the maximum value in the axis C Output: C DX: step size C IMPLICIT REAL*8(A-H,O-Z) DDX = XMAX - XMIN DX = INT( DLOG10( DDX ) ) IF( DDX.LT.1.0D0 ) DX = DX - 1 DX = 10.0D0**DX IF( DDX/DX.LT.3.0D0 ) DX = DX/2.0 IF( DDX/DX.LT.3.0D0 ) DX = DX/2.0 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 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 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 DOUBLE PRECISION FUNCTION GAUSS( X,PARAM ) C C ... Gaussian (normal) distribution ... C C Inputs: C X: C PARAM(1): mean C PARAM(2): variance C Output: C GAUSS: density at X C IMPLICIT REAL*8(A-H,O-Z) DIMENSION PARAM(2) DATA C1 /2.506628275D0/ C GAUSS = DEXP( -(X-PARAM(1))**2/(2*PARAM(2)) )/(C1*DSQRT(PARAM(2))) RETURN E N D DOUBLE PRECISION FUNCTION PEARSN( X,PARAM ) C C ... Pearson family of distributions ... C C Inputs: C X: C PARAM(1): location parameter, mu C PARAM(2): dispersion parameter, tau square C PARAM(3): shape parameter C Output: C PEARSN: density at X C IMPLICIT REAL*8(A-H,O-Z) DIMENSION PARAM(3) DATA PI/3.1415926535D0/ C PEARSN = DGAMMA(PARAM(3))/DGAMMA(PARAM(3)-0.5D0) * /DSQRT(PI)*PARAM(2)**(PARAM(3)-0.5D0) * /((X-PARAM(1))**2 + PARAM(2))**PARAM(3) RETURN C END DOUBLE PRECISION FUNCTION DGAMMA( X ) C C ... Gamma function ... C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(0:10) DATA A /0.999999999999269D0, 0.42278433696202D0, * 0.41184025179616D0, 0.08157821878492D0, 0.0742379070629D0, * -0.00021090746731D0, 0.01097369584174D0,-0.00246674798054D0, * 0.00153976810472D0,-0.00034423420456D0, 0.00006771057117D0/ C DGAM = 1.0D0 Y = X IF( X.GT.3.0D0 ) THEN 10 Y = Y-1 DGAM = DGAM*Y IF( Y.GT.3.0D0 ) GO TO 10 END IF IF( X.LE.2.0D0 ) THEN 20 DGAM = DGAM/Y Y = Y+1 IF( Y.LE.2.0D0 ) GO TO 20 END IF C Z = 1.0D0 SUM = 0.0D0 DO 30 I=0,10 SUM = SUM + A(I)*Z 30 Z = Z*(Y-2) DGAMMA = DGAM*SUM RETURN E N D DOUBLE PRECISION FUNCTION DBLEXP( X,PARAM ) C C ... double exponential distribution f(x) = exp(x - exp(x)) ... C C Inputs: C X: C Output: C DBLEXP: density at X C IMPLICIT REAL*8(A-H,O-Z) C DBLEXP = DEXP( X-DEXP(X) ) RETURN C E N D DOUBLE PRECISION FUNCTION USERI( X,PARAM ) C C ... User supplied density function ... C (The following is an example of two-sided exponential dist.) C C Inputs: C X: C PARAM(1): mean C PARAM(2): lambda C Output: C USERI: density at X C IMPLICIT REAL*8(A-H,O-Z) DIMENSION PARAM(2) C SIGMA = DSQRT( PARAM(2) ) USERI = SIGMA*DEXP( -SIGMA*DABS(X-PARAM(1)) )/2 RETURN E N D DOUBLE PRECISION FUNCTION UNIF( X,PARAM ) C C ... uniform distribution f(x) = 1 ... C C Output: C UNIF: density at X (=1) C IMPLICIT REAL*8(A-H,O-Z) C UNIF = 1.0D0 RETURN E N D SUBROUTINE PLOT3A( N,M,WX,WY,WZ,ANGLE,X0,X1,DX,Y0,Y1,DY,Z0,Z1, * DZ,KN,K,YH ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION YH(2000) C PI = 3.1415926535D0 THETA = ANGLE * ( PI/180.0 ) C CALL AXISXY( 0.0D0,0.0D0,WX,WY,X0,X1,Y0,Y1,DX,DY,0.2D0,0,5,5 ) CALL AXISZ( WX,WZ,THETA,Z0,Z1,DZ,0.2D0,1 ) ZS = WZ*DSIN(THETA) KN = K*(N-1) + 1 DO 30 I=1,2000 30 YH(I) = Y0 + 2*(ZS/N)/(WY/(Y1-Y0)) RETURN E N D SUBROUTINE PLOT3B( Y,N,M,WX,WY,WZ,ANGLE,Y0,Y1,K,YH,ZS,ZC ) IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(N+1), YH(2000) C PI = 3.1415926535D0 THETA = ANGLE * ( PI/180.0 ) C ZS = WZ*DSIN(THETA) ZC = WZ*DCOS(THETA) KN = K*(N-1) + 1 SHIFT = (ZC/M) / (WX/(KN-1)) IS = SHIFT DS = SHIFT-IS KN1 = KN - IS DO 40 I = 1,KN1 L = I+IS 40 YH(I) = YH(L)+ ((YH(L+1)-YH(L))*DS) - (ZS/M)/(WY/(Y1-Y0)) DO 50 I=KN1+1,KN 50 YH(I) = Y0 DO 60 I=1,KN L = (I-1)/K + 1 IF( I.LT.KN ) W = Y(L) +(( Y(L+1) - Y(L) )*MOD(I-1,K))/K IF( I.EQ.KN ) W = Y(L) IF(L.EQ.N) W = Y(L) 60 YH(I) = DMAX1( YH(I),W) ZC1 = ZC/M ZS1 = ZS/M CALL PLOT( SNGL(ZC1),SNGL(ZS1),-3 ) CALL PLOTY( YH,KN,Y0,Y1,WX,WY,0,1 ) C RETURN E N D SUBROUTINE AXISZ( WX,WZ,THETA,Z0,Z1,DZ,DWC,IZ ) IMPLICIT REAL*8(A-H,O-Z) REAL*4 WC C CALL NEWPEN(2) CALL PLOT( SNGL(WX),0.0,3 ) CALL PLOT( SNGL(WX+WZ*DCOS(THETA)),SNGL(WZ*DSIN(THETA)),2 ) C WC = DWC NZ = (Z1-Z0)/DZ + 0.001 MZ = -DLOG10(DZ) IF ( MZ.LE.0) MZ = -1 D = WZ*DZ /( Z1-Z0) DO 10 I = 1,NZ+1 ZN = Z0 + DZ*(I-1) ZZ = D*(I-1) CALL PLOT(SNGL(WX+ZZ*DCOS(THETA)),SNGL(ZZ*DSIN(THETA)),3 ) CALL PLOT(SNGL(WX+ZZ*DCOS(THETA)+(WC/2)),SNGL(ZZ*DSIN(THETA)),2) IF( MOD(I,IZ).EQ.1 .OR. IZ.EQ.1 ) * CALL NUMBER( SNGL(WX+ZZ*DCOS(THETA)+2*WC),SNGL(ZZ*DSIN(THETA) * -(WC/2)),WC,ZN,0.0,MZ ) 10 CONTINUE CALL NEWPEN( 1 ) RETURN E N D