Help!Help!!一个关于求解安托尼方程参数的程序
小白是化学系研究生，现跟着老板做计算，老板让学fortran,然后就甩了这么个程序给我。小白折腾了很久这个程序也只看懂了皮毛，有好多问题希望各路大神能帮忙解答，小白在此先行谢过！！这个程序是用fortran77编的，语法比较古老，还望海涵！
Q1：求出AntA,AntB,AntC的算法是什么？
Q2：那个求逆的算法是什么？小白感觉和一般的高斯消元不一样
Q3：下面那个计算最优参数的那个子程序，我就觉得跟天书一样，taila
C This programme is used to correlate the Antoine constants C
C for vapor pressure of pure components C
C C
C The Antoine equation of vapor pressure is C
C C
C lg(p)=AntA-AntB/(t+AntC) C
C where C
C p in kPa and t in C
PROGRAM LIU53
CHARACTER DATFILE*12,RESFILE*12,SYSNAME*30,SYSFORM*10
DIMENSION A(7,7),B(7),X(7),PC(60),DP(60)
COMMON /TP/ TE(60),PE(60)
WRITE(*,'(1X////1X)')
WRITE(*,*) 'Input the file name of data file'
READ(*,'(A12)') DATFILE
OPEN(4,FILE=DatFile,STATUS='OLD')
READ(4,*) SYSNAME,SYSFORM
READ(4,*) ND,IPUNIT,ITUNIT
C
C ITUNITindicating the unit of temperature
C =0 for K =1 for C
C
C IPUNITindicating the unit of pressure
C =0 forkPa
C =1 formmHg
C =2 foratm
C =3 forbar
C =4 forMPa
C
READ(4,*) (TE(I),PE(I),I=1,ND)
CLOSE(4)
TMIN=TE(1)
TMAX=TE(ND)
DO 10 I=1,ND
IF (ITUNIT.EQ.0) TE(I)=TE(I)-273.15
IF (IPUNIT.EQ.1) PE(I)=PE(I)/760.*101.325
IF (IPUNIT.EQ.2) PE(I)=PE(I)*101.325
IF (IPUNIT.EQ.3) PE(I)=PE(I)*100.0
IF (IPUNIT.EQ.4) PE(I)=PE(I)*1000.0
10 CONTINUE
DO 180 I=1,7
DO 190 J=1,7
A(I,J)=0.0
190 CONTINUE
B(I)=0.0
180 CONTINUE
DO 200 I=1,ND
PLG=ALOG(PE(I))/2.302585
T=TE(I)
B(1)=B(1)+T*T*PLG
B(2)=B(2)+T*PLG*PLG
B(3)=B(3)+T*PLG
A(1,1)=A(1,1)+T*T
A(1,3)=A(1,3)+T
A(2,2)=A(2,2)+PLG*PLG
A(2,3)=A(2,3)+PLG
200 CONTINUE
A(1,2)=B(3)
A(2,1)=A(1,2)
A(3,1)=A(1,3)
A(3,2)=A(2,3)
A(3,3)=FLOAT(ND)
N=3
CALL INVB(N,A)
DO 210 I=1,N
X(I)=0.0
DO 220 J=1,N
X(I)=X(I)+A(I,J)*B(J)
220 CONTINUE
210 CONTINUE
ANTA=X(1)
ANTC=-X(2)
ANTB=ANTA*ANTC-X(3)
N=3
X(1)=ANTA
X(2)=ANTB
X(3)=ANTC
M=ND
CALL SM(N,M,X)
ANTA=X(1)
ANTB=X(2)
ANTC=X(3)
RDP=0.0
DO 230 I=1,ND
PC(I)=EXP((ANTA-ANTB/(TE(I)+ANTC))*2.302585)
DP(I)=PE(I)-PC(I)
RDP=RDP+ABS(DP(I))/PE(I)
230 CONTINUE
RDP=RDP*100.0/FLOAT(ND)
WRITE(*,'(1X////1X)')
WRITE(*,*) 'Input the file name of result file'
READ(*,'(A12)') RESFILE
OPEN(4,FILE=ResFile,STATUS='NEW')
WRITE(4,240)
240 FORMAT(5X,'Correlation of vapor pressure for component')
WRITE(4,250) SYSNAME,SYSFORM
250 FORMAT(/5X,A30,3X,A10)
WRITE(4,260) ANTA,ANTB,ANTC,TMIN,TMAX
260 FORMAT(/5X,'ANTA=',F7.5,2X,'ANTB=',F9.3,2X,'ANTC=',F7.3,
1 2X,'TMIN=',F6.2,2X,'TMAX=',F6.2)
WRITE(4,270)
270 FORMAT(/7X,'T/K',5X,'Pe/kPa',3X,'Pc/kPa',2X,'DP/kPa')
WRITE(4,280) (I,TE(I),PE(I),PC(I),DP(I),I=1,ND)
280 FORMAT(2X,I2,2X,F6.2,2X,F7.3,2X,F7.3,2X,F6.3)
WRITE(4,290) RDP
290 FORMAT(/5X,'RDP=',F6.2,'%')
CLOSE(4)
STOP
END
SUBROUTINE VFRE(X,N,M,Y)
DIMENSION X(5),Y(60)
COMMON /TP/ TE(60),PE(60)
ANTA=X(1)
ANTB=X(2)
ANTC=X(3)
DO 10 I=1,M
PC=EXP((ANTA-ANTB/(TE(I)+ANTC))*2.302585)
Y(I)=1.0-PC/PE(I)
10 CONTINUE
RETURN
END
C***********************************************************************C
C C
C SUBROUTINE INVB(N,A) C
C C
C PURPOSE C
C C
C calculate the inverse matrix of matrix A C
C C
C A ( N , N ) C
C C
C The result is also put in A(N,N) C
C C
C C
C***********************************************************************C
C
SUBROUTINE INVB(N,A)
C
DIMENSION A(7,7)
IF (N-2) 11,11,22
11 AI=A(1,1)*A(2,2)-A(1,2)*A(2,1)
AI=1.0/AI
AJ=A(1,1)
A(1,1)=A(2,2)*AI
A(2,2)=AI*AJ
A(2,1)=-A(2,1)*AI
A(1,2)=-A(1,2)*AI
RETURN
22 DO 20 K=1,N
P=1.0/A(K,K)
A(K,K)=P
DO 30 I=1,N
IF (I.NE.K) A(I,K)=-P*A(I,K)
30 CONTINUE
DO 40 I=1,N
IF (I.EQ.K) GOTO 60
DO 50 J=1,N
IF (J.NE.K) A(I,J)=A(I,J)+A(I,K)*A(K,J)
50 CONTINUE
60 CONTINUE
40 CONTINUE
DO 70 J=1,N
IF (J.NE.K) A(K,J)=P*A(K,J)
70 CONTINUE
20 CONTINUE
RETURN
END
C***********************************************************************C
C C
C SUBROUTINE SM(N,M,X) C
C C
C PURPOSE C
C C
C calculating the optimum parameter values C
C C
C C
C N : the number of the paramrters required estimation C
C M : the number of the equation C
C X(N) : the parameters required estimation C
C put the start value at the first C
C put the optimum values in the last C
C SIGAMY : external function procedure for calculating the C
C objective function value C
C C
C***********************************************************************C
C
SUBROUTINE SM(N,M,X)
C
INTEGER T,A2,B2,C2,N,M,K,TT
DIMENSION X(5),R1(5),R2(5),R3(5),Q(6),P(6,5)
V=.01
MAXP=220
TT=0
AN=N
10 DO 20 T=1,N+1
DO 30 K=1,N
IF ((T-1).NE.K) THEN
R1(K)=X(K)
P(T,K)=X(K)
ELSE IF (ABS(X(K)).LT.1.E-6) THEN
R1(K)=1.E-6
P(T,K)=1.E-6
ELSE
R1(K)=(1.0+V)*X(K)
P(T,K)=(1.0+V)*X(K)
ENDIF
30 CONTINUE
Q(T)=SIGAMY(N,M,R1)
20 CONTINUE
40 A1=0.
B1=0.
C1=1.E18
DO 50 T=1,N+1
IF (A1.GT.Q(T)) GO TO 60
A1=Q(T)
A2=T
60 IF (C1.LT.Q(T)) GO TO 50
C1=Q(T)
C2=T
50 CONTINUE
DO 70 T=1,N+1
IF (T.NE.A2.AND.B1.LT.Q(T)) THEN
B1=Q(T)
B2=T
ENDIF
70 CONTINUE
DO 80 K=1,N
X(K)=P(C2,K)
80 R1(K)=0.
TT=TT+1
SMN=0.
DO 85 K=1,N+1
85 SMN=SMN+Q(K)/(AN+1.)
SD=0.
ANTA=0.
DO 87 K=1,N+1
SD=SD+(Q(K)/SMN-1.)**2
87 ANTA=ANTA+(Q(K)-SMN)**2
ANTA=SQRT(ANTA/(AN+1.))
SD=SQRT(SD/(AN+1.0))
C WRITE(*,88)ANTA,SD,TT
IF (SD.LE.1.E-7.OR.ANTA.LE.1.E-7) GO TO 170
IF (TT.GE.MAXP.OR.V.LT.0.0001) GO TO 170
88 FORMAT(1X,5HANTA=,E10.2,2X,3HSD=,E10.2,2X,3HTT=,I3)
IF (C1.LE.1.E-7.OR.ABS(1.-A1/C1).LE.1.E-7) GO TO 170
DO 100 T=1,N+1
DO 90 K=1,N
IF (T.NE.A2) R1(K)=R1(K)+P(T,K)/AN
90 CONTINUE
100 CONTINUE
DO 110 K=1,N
110 R2(K)=2.*R1(K)-P(A2,K)
S=SIGAMY(N,M,R2)
IF (S.GE.A1) GO TO 150
DO 120 K=1,N
R3(K)=3.*R1(K)-2.*P(A2,K)
120 P(A2,K)=R2(K)
Q(A2)=S
U=SIGAMY(N,M,R3)
IF (U.GE.S) GO TO 40
130 DO 140 K=1,N
140 P(A2,K)=R3(K)
Q(A2)=U
GO TO 40
150 DO 160 K=1,N
160 R3(K)=R1(K)/2.+P(A2,K)/2.
U=SIGAMY(N,M,R3)
IF (U.LT.B1) GO TO 130
V=V/2.0
GOTO 10
170 CONTINUE
RETURN
END
FUNCTION SIGAMY(N,M,X)
DIMENSION X(5),Y(60)
CALL VFRE(X,N,M,Y)
A=0.0
DO 10 I=1,M
10 A=A+Y(I)*Y(I)
A=SQRT(A/FLOAT(M))
C WRITE(*,5) A,(X(J),J=1,N)
5 FORMAT(2X,'S9=',E12.6,2X,'X(1)=',E12.6,2X,'X(2)=',E12.6,2X,
*'X(3)=',E12.6)
SIGAMY=A
RETURN
END
页:
[1]