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