      REAL FAP(21),P(20),XT(50),YT(50,5),TX(250),CY(250)
      CHARACTER XS*10,YS*10,NA*10
      NS=50
      CALL FILTA(FAP)
      WRITE(*,*) '############################'
      WRITE(*,*) '# SIMULATION BY FILT(V2.5) #'
      WRITE(*,*) '############################'
      WRITE(*,*)
      CALL FMA(NP,LN,NL,TX,CY,NA)
      WRITE(*,*) 'NUMBER OF EQUATIONS ?'
      READ(*,*) NE
      WRITE(*,*) 'NUMBER OF PARAMETERS ?'
      READ(*,*) M
      WRITE(*,*) 'X0,X9 ?'
      READ(*,*) T0,T9
      WRITE(*,*) 'Y0,Y9?'
      READ(*,*) Y0,Y9
      WRITE(*,*) 'NAME OF X-AXIS?'
      READ(*,'(A)') XS
      WRITE(*,*) 'NAME OF Y-AXIS?'
      READ(*,'(A)') YS
      DT=(T9-T0)/(NS-1)
    5 DO 10 I=1,M
      WRITE(*,'('' P('',I2,'') ?'')') I
   10 READ(*,*) P(I)
      DO 20 I=1,NS
      T=T0+DT*(I-1)
      XT(I)=T
      DO 20 J=1,NE
      CALL FILTB(J,T,P,CP,FAP)
      YT(I,J)=CP
   20 CONTINUE
      DO 30 I=1,NS
      WRITE(*,'('' T='',F10.5,'',  CP='',5F10.5)')
     &XT(I),(YT(I,J),J=1,NE)
   30 CONTINUE
      DO 180 I=1,NS
      DO 180 J=1,NE
      IF(Y9.LT.YT(I,J))Y9=YT(I,J)
  180 CONTINUE
      CALL PLOT2(NS,NE,XT,YT,T0,T9,Y0,Y9,XS,YS,NP,TX,CY)
      WRITE(*,*)'SAVE(Y/N) ?'
      READ(*,'(A)') YN
      IF(YN.EQ.'Y') CALL SAVE2(NS,NE,XT,YT,NP,TX,CY,M,P)
      GOTO 5
      END
      SUBROUTINE FILTA(FAP)
      DIMENSION FAP(21),FCP(21)
      NFP=6
      FCP(1)=1.
      DO 10 I=1,NFP
   10 FCP(I+1)=FCP(I)*FLOAT(NFP+1-I)/FLOAT(I)
      FAP(NFP+2)=0.
      DO 20 I=NFP,0,-1
   20 FAP(I+1)=FAP(I+2)+FCP(I+1)
      END
      SUBROUTINE FILTB(J,T,P,CP,FAP)
      REAL P(20),FAP(21)
      COMPLEX CCP,S
      PI=3.141593
      NFP=6
      NFA=6
      NFK=50
      CP=0.
      IF(T.EQ.0.) T=1.E-5
      FAX=EXP(FLOAT(NFA))/T
      DO 10 I=1,NFK
      S=CMPLX(FLOAT(NFA)/T,(FLOAT(I)-.5)*PI/T)
      CALL FFUNC(J,P,S,CCP)
      FF=AIMAG(CCP)
      IF(MOD(I,2).NE.0) FF=-FF
      CP=CP+FF
   10 CONTINUE
      CP=CP*FAP(1)
      DO 20 I=1,NFP
      S=CMPLX(FLOAT(NFA)/T,(FLOAT(NFK+I)-.5)*PI/T)
      CALL FFUNC(J,P,S,CCP)
      FF=AIMAG(CCP)
      IF(MOD(NFK+I,2).NE.0) FF=-FF
      CP=CP+FAP(I+1)*FF
   20 CONTINUE
      CP=FAX*CP/FAP(1)
      END
      SUBROUTINE PLOT2(NS,NE,XT,YT,X0,X9,Y0,Y9,XS,YS,NP,TX,CY)
      DIMENSION XT(50),YT(50,5),XG(6),TX(250),CY(250)
      CHARACTER XS*10,YS*10,IP(51,17)
      XM=50.0
      YM=16.0
      DO 100 I=1,51
      DO 100 J=1,17
      IP(I,J)=' '
  100 CONTINUE
      DO 110 I=1,NS
      IX=INT(XM*(XT(I)-X0)/(X9-X0)+0.5)+1
      DO 110 J=1,NE
      IY=INT(YM*(Y9-YT(I,J))/(Y9-Y0)+0.5)+1
      IP(IX,IY)=CHAR(J+48)
  110 CONTINUE
      DO 107 I=1,NP
      IX=INT(XM*(TX(I)-X0)/(X9-X0)+0.5)+1
      IY=INT(YM*(Y9-CY(I))/(Y9-Y0)+0.5)+1
      IP(IX,IY)=CHAR(79)
  107 CONTINUE
      WRITE(*,'(5X,A10)') YS
      DO 130 K=1,17
      JJ=K-4*INT(K/4)
      IF(JJ.EQ.1) THEN
      YY=(Y9-Y0)*(16-K+1)/16.0+Y0
      IF((YY.GE.10000.0).OR.((YY.LT.0.1).AND.(YY.NE.0.0))) THEN
      WRITE(*,'(1X,E12.4,''+'',51A1)') YY,(IP(I,K),I=1,51)
      ELSE
      WRITE(*,'(5X,F8.3,''+'',51A1)') YY,(IP(I,K),I=1,51)
      ENDIF
      ELSE
      WRITE(*,'(13X,''I'',51A1)') (IP(I,K),I=1,51)
      ENDIF
  130 CONTINUE
      WRITE(*,'(13X,''++'',5A10)') ('---------+',I=1,5)
      DO 140 I=1,6
      XG(I)=(X9-X0)*(I-1)/5.0+X0
  140 CONTINUE
      WRITE(*,'(7X,6F10.2)') (XG(I),I=1,6)
      WRITE(*,'(35X,A10)') XS
      END
      SUBROUTINE TMP
      CHARACTER TP
      WRITE(*,*) 'DO YOU CONTINUE(Y/N)?'
      READ(*,'(A1)') TP
      IF(TP.EQ.'N') STOP
      END
C   ====== FILE I/O ======= 
      SUBROUTINE FMA(N,LN,NL,X,Y,FN)
      DIMENSION NL(5),X(50),Y(50)
      CHARACTER FN*10,IC*1
  100 WRITE(*,*) 'E)XIT, I)NPUT, R)EAD, W)RITE, S)TOP ?'
      READ(*,'(A)') IC
      IF(IC.EQ.'S') STOP
      IF(IC.EQ.'E') RETURN
      IF(IC.EQ.'I') GOTO 1100
      IF(IC.EQ.'R') GOTO 1200
      IF(IC.EQ.'W') GOTO 1300
      GOTO 100
 1100 WRITE(*,*) 'NUMBER OF LINES (1-5) ?'
      READ(*,*) LN
      N=0
      IB=0
      DO 1115 J=1,LN
      WRITE(*,'('' NUMBER OF POINTS ('',I2,'') ?'')') J
      READ(*,*) NL(J)
      N=N+NL(J)
      DO 1127 I=1,NL(J)
      WRITE(*,1125) J,I,J,I
 1125 FORMAT(' T',I2,'(',I2,')  , CP',I2,'(',I2,') ?')
      READ(*,*) X(IB+I),Y(IB+I)
 1127 CONTINUE
      IB=IB+NL(J)
 1115 CONTINUE 
      GOTO 100
 1200 WRITE(*,*) 'FILE NAME TO READ ?'
      READ(*,'(A)') FN
      OPEN(UNIT=5,FILE=FN)
      READ(5,*) LN
      N=0
      IB=0
      DO 1220 J=1,LN
      READ(5,*) NL(J)
      N=N+NL(J)
      WRITE(*,*)
      DO 1210 I=1,NL(J)
      READ(5,*) X(IB+I),Y(IB+I)
      WRITE(*,1225) J,I,X(IB+I),J,I,Y(IB+I)
 1225 FORMAT(' T',I2,'(',I2,')=',F10.5,'  , CP',I2,'(',I2,')=',F10.5)
 1210 CONTINUE
      IB=IB+NL(J)
 1220 CONTINUE
      CLOSE(UNIT=5)
      GOTO 100
 1300 WRITE(*,*) 'FILE NAME TO WRITE ?'
      READ(*,'(A)') FN
      OPEN(UNIT=6,FILE=FN,STATUS='NEW')
      WRITE(6,*) LN
      N=0
      IB=0
      DO 1315 J=1,LN
      WRITE(6,*) NL(J)
      N=N+NL(J)
      DO 1327 I=1,NL(J)
      WRITE(6,'(F12.5,'','',F12.5)') X(IB+I),Y(IB+I)
 1327 CONTINUE
      IB=IB+NL(J)
 1315 CONTINUE 
      CLOSE(UNIT=6)
      GOTO 100
      END
C  ======  GRAPHIC DATA SAVE =======
      SUBROUTINE SAVE2(NS,NE,XT,YT,NP,TX,CY,M,P)
      DIMENSION XT(NS),YT(NS,NE),TX(NP),CY(NP),P(20)
      CHARACTER NA*12
      WRITE(*,*)'FILE NAME ?'
      READ(*,'(A)') NA
      OPEN(UNIT=6,FILE=NA,STATUS='NEW')
      WRITE(6,'(I10,'','',I10)') NS,NE
      DO 10 I=1,NS
      WRITE(6,'(6(F12.5,'',''))') XT(I),(YT(I,J),J=1,NE)
   10 CONTINUE
      WRITE(6,*) NP
      DO 20 I=1,NP
      WRITE(6,'(F12.5,'','',F12.5)') TX(I),CY(I)
   20 CONTINUE
      DO 30 I=1,M
      WRITE(6,'(''P('',I2,'')=,'',F12.5)')I,P(I)
   30 CONTINUE
      CLOSE(UNIT=6)
      END
C-----------------------------------------------------------------------
      SUBROUTINE FFUNC(J,P,S,CCP)
      DIMENSION P(20)
      COMPLEX S,CCP
      GOTO(100,200,300,400,500)J
  100 CCP=P(1)/(S+P(2))
      RETURN
  200 CCP=P(1)*P(3)*P(4)/(S+P(2))/(S+P(3))
      RETURN
  300 CONTINUE
      RETURN
  400 CONTINUE
      RETURN
  500 CONTINUE
      RETURN
      END
