Help!Help!!一个关于求解安托尼方程参数的程序

zhuxiaobao

小白是化学系研究生，现跟着老板做计算，老板让学fortran,然后就甩了这么个程序给我。小白折腾了很久这个程序也只看懂了皮毛，有好多问题希望各路大神能帮忙解答，小白在此先行谢过！！
这个程序是用fortran77编的，语法比较古老，还望海涵！
Q1：求出AntA,AntB,AntC的算法是什么？
Q2：那个求逆的算法是什么？小白感觉和一般的高斯消元不一样
Q3：下面那个计算最优参数的那个子程序，我就觉得跟天书一样，taila

[Fortran] 纯文本查看 复制代码

查看代码

打印?

001	C      This programme is used to correlate the Antoine constants       C

002	C       for vapor pressure of pure components                           C

003	C                                                                       C

004	C       The Antoine equation of vapor pressure is                       C

005	C                                                                       C

006	C       lg(p)=AntA-AntB/(t+AntC)                                        C

007	C       where                                                           C

008	C       p in kPa and t in C

009	PROGRAM LIU53

010	CHARACTER DATFILE*12,RESFILE*12,SYSNAME*30,SYSFORM*10

011	DIMENSION A(7,7),B(7),X(7),PC(60),DP(60)

012	COMMON /TP/ TE(60),PE(60)

013	WRITE(*,'(1X////1X)')

014	WRITE(*,*) 'Input the file name of data file'

015	READ(*,'(A12)') DATFILE

016	OPEN(4,FILE=DatFile,STATUS='OLD')

017	READ(4,*) SYSNAME,SYSFORM

018	READ(4,*) ND,IPUNIT,ITUNIT

019

C

020	C       ITUNIT  indicating the unit of temperature

021	C               =0      for K           =1      for C

022

C

023	C       IPUNIT  indicating the unit of pressure

024	C               =0      for  kPa

025	C               =1      for  mmHg

026	C               =2      for  atm

027	C               =3      for  bar

028	C               =4      for  MPa

029

C

030	READ(4,*) (TE(I),PE(I),I=1,ND)

031

CLOSE(4)

032	TMIN=TE(1)

033	TMAX=TE(ND)

034	DO 10 I=1,ND

035	IF (ITUNIT.EQ.0) TE(I)=TE(I)-273.15

036	IF (IPUNIT.EQ.1) PE(I)=PE(I)/760.*101.325

037	IF (IPUNIT.EQ.2) PE(I)=PE(I)*101.325

038	IF (IPUNIT.EQ.3) PE(I)=PE(I)*100.0

039	IF (IPUNIT.EQ.4) PE(I)=PE(I)*1000.0

040	10      CONTINUE

041	DO 180 I=1,7

042	DO 190 J=1,7

043	A(I,J)=0.0

044	190     CONTINUE

045

B(I)=0.0

046	180     CONTINUE

047	DO 200 I=1,ND

048	PLG=ALOG(PE(I))/2.302585

049

T=TE(I)

050	B(1)=B(1)+T*T*PLG

051	B(2)=B(2)+T*PLG*PLG

052	B(3)=B(3)+T*PLG

053	A(1,1)=A(1,1)+T*T

054	A(1,3)=A(1,3)+T

055	A(2,2)=A(2,2)+PLG*PLG

056	A(2,3)=A(2,3)+PLG

057	200     CONTINUE

058	A(1,2)=B(3)

059	A(2,1)=A(1,2)

060	A(3,1)=A(1,3)

061	A(3,2)=A(2,3)

062	A(3,3)=FLOAT(ND)

063

N=3

064	CALL INVB(N,A)

065	DO 210 I=1,N

066

X(I)=0.0

067	DO 220 J=1,N

068	X(I)=X(I)+A(I,J)*B(J)

069	220     CONTINUE

070	210     CONTINUE

071

ANTA=X(1)

072	ANTC=-X(2)

073	ANTB=ANTA*ANTC-X(3)

074

N=3

075

X(1)=ANTA

076

X(2)=ANTB

077

X(3)=ANTC

078

M=ND

079	CALL SM(N,M,X)

080

ANTA=X(1)

081

ANTB=X(2)

082

ANTC=X(3)

083

RDP=0.0

084	DO 230 I=1,ND

085	PC(I)=EXP((ANTA-ANTB/(TE(I)+ANTC))*2.302585)

086	DP(I)=PE(I)-PC(I)

087	RDP=RDP+ABS(DP(I))/PE(I)

088	230     CONTINUE

089	RDP=RDP*100.0/FLOAT(ND)

090	WRITE(*,'(1X////1X)')

091	WRITE(*,*) 'Input the file name of result file'

092	READ(*,'(A12)') RESFILE

093	OPEN(4,FILE=ResFile,STATUS='NEW')

094	WRITE(4,240)

095	240     FORMAT(5X,'Correlation of vapor pressure for component')

096	WRITE(4,250) SYSNAME,SYSFORM

097	250     FORMAT(/5X,A30,3X,A10)

098	WRITE(4,260) ANTA,ANTB,ANTC,TMIN,TMAX

099	260     FORMAT(/5X,'ANTA=',F7.5,2X,'ANTB=',F9.3,2X,'ANTC=',F7.3,

100	1          2X,'TMIN=',F6.2,2X,'TMAX=',F6.2)

101	WRITE(4,270)

102	270     FORMAT(/7X,'T/K',5X,'Pe/kPa',3X,'Pc/kPa',2X,'DP/kPa')

103	WRITE(4,280) (I,TE(I),PE(I),PC(I),DP(I),I=1,ND)

104	280     FORMAT(2X,I2,2X,F6.2,2X,F7.3,2X,F7.3,2X,F6.3)

105	WRITE(4,290) RDP

106	290     FORMAT(/5X,'RDP=',F6.2,'%')

107

CLOSE(4)

108

STOP

109

END

110

111	SUBROUTINE VFRE(X,N,M,Y)

112	DIMENSION X(5),Y(60)

113	COMMON /TP/ TE(60),PE(60)

114

ANTA=X(1)

115

ANTB=X(2)

116

ANTC=X(3)

117	DO 10 I=1,M

118	PC=EXP((ANTA-ANTB/(TE(I)+ANTC))*2.302585)

119	Y(I)=1.0-PC/PE(I)

120	10      CONTINUE

121

RETURN

122

END

123

124

125	C***********************************************************************C

126	C                                                                       C

127	C       SUBROUTINE INVB(N,A)                                            C

128	C                                                                       C

129	C       PURPOSE                                                         C

130	C                                                                       C

131	C       calculate the inverse matrix of matrix A                        C

132	C                                                                       C

133	C                       A ( N , N )                                     C

134	C                                                                       C

135	C       The result is also put in A(N,N)                                C

136	C                                                                       C

137	C                                 C

138	C***********************************************************************C

139

C

140	SUBROUTINE INVB(N,A)

141

C

142	DIMENSION A(7,7)

143	IF (N-2) 11,11,22

144	11    AI=A(1,1)*A(2,2)-A(1,2)*A(2,1)

145

AI=1.0/AI

146

AJ=A(1,1)

147	A(1,1)=A(2,2)*AI

148	A(2,2)=AI*AJ

149	A(2,1)=-A(2,1)*AI

150	A(1,2)=-A(1,2)*AI

151

RETURN

152	22    DO 20 K=1,N

153	P=1.0/A(K,K)

154

A(K,K)=P

155	DO 30 I=1,N

156	IF (I.NE.K) A(I,K)=-P*A(I,K)

157	30    CONTINUE

158	DO 40 I=1,N

159	IF (I.EQ.K) GOTO 60

160	DO 50 J=1,N

161	IF (J.NE.K) A(I,J)=A(I,J)+A(I,K)*A(K,J)

162	50    CONTINUE

163	60    CONTINUE

164	40    CONTINUE

165	DO 70 J=1,N

166	IF (J.NE.K) A(K,J)=P*A(K,J)

167	70    CONTINUE

168	20    CONTINUE

169

RETURN

170

END

171

172	C***********************************************************************C

173	C                                                                       C

174	C       SUBROUTINE SM(N,M,X)                                            C

175	C                                                                       C

176	C       PURPOSE                                                         C

177	C                                                                       C

178	C       calculating the optimum parameter values                        C

179	C                                                C

180	C                                                                       C

181	C        N : the number of the paramrters required estimation           C

182	C        M : the number of the equation                                 C

183	C        X(N) : the parameters required estimation                      C

184	C               put the start value at the first                        C

185	C               put the optimum values in the last                      C

186	C       SIGAMY : external function procedure for calculating the        C

187	C                 objective function value                              C

188	C                                                                       C

189	C***********************************************************************C

190

C

191	SUBROUTINE SM(N,M,X)

192

C

193	INTEGER T,A2,B2,C2,N,M,K,TT

194	DIMENSION X(5),R1(5),R2(5),R3(5),Q(6),P(6,5)

195

V=.01

196

MAXP=220

197

TT=0

198

AN=N

199	10     DO 20 T=1,N+1

200	DO 30 K=1,N

201	IF ((T-1).NE.K) THEN

202	R1(K)=X(K)

203	P(T,K)=X(K)

204	ELSE IF (ABS(X(K)).LT.1.E-6) THEN

205	R1(K)=1.E-6

206	P(T,K)=1.E-6

207

ELSE

208	R1(K)=(1.0+V)*X(K)

209	P(T,K)=(1.0+V)*X(K)

210

ENDIF

211	30     CONTINUE

212	Q(T)=SIGAMY(N,M,R1)

213	20     CONTINUE

214	40     A1=0.

215

B1=0.

216

C1=1.E18

217	DO 50 T=1,N+1

218	IF (A1.GT.Q(T)) GO TO 60

219

A1=Q(T)

220

A2=T

221	60    IF (C1.LT.Q(T)) GO TO 50

222

C1=Q(T)

223

C2=T

224	50     CONTINUE

225	DO 70 T=1,N+1

226	IF (T.NE.A2.AND.B1.LT.Q(T)) THEN

227

B1=Q(T)

228

B2=T

229

ENDIF

230	70     CONTINUE

231	DO 80 K=1,N

232	X(K)=P(C2,K)

233	80     R1(K)=0.

234

TT=TT+1

235

SMN=0.

236	DO 85 K=1,N+1

237	85     SMN=SMN+Q(K)/(AN+1.)

238

SD=0.

239

ANTA=0.

240	DO 87 K=1,N+1

241	SD=SD+(Q(K)/SMN-1.)**2

242	87     ANTA=ANTA+(Q(K)-SMN)**2

243	ANTA=SQRT(ANTA/(AN+1.))

244	SD=SQRT(SD/(AN+1.0))

245	C        WRITE(*,88)ANTA,SD,TT

246	IF (SD.LE.1.E-7.OR.ANTA.LE.1.E-7) GO TO 170

247	IF (TT.GE.MAXP.OR.V.LT.0.0001) GO TO 170

248	88     FORMAT(1X,5HANTA=,E10.2,2X,3HSD=,E10.2,2X,3HTT=,I3)

249	IF (C1.LE.1.E-7.OR.ABS(1.-A1/C1).LE.1.E-7) GO TO 170

250	DO 100 T=1,N+1

251	DO 90 K=1,N

252	IF (T.NE.A2) R1(K)=R1(K)+P(T,K)/AN

253	90     CONTINUE

254	100    CONTINUE

255	DO 110 K=1,N

256	110    R2(K)=2.*R1(K)-P(A2,K)

257	S=SIGAMY(N,M,R2)

258	IF (S.GE.A1) GO TO 150

259	DO 120 K=1,N

260	R3(K)=3.*R1(K)-2.*P(A2,K)

261	120    P(A2,K)=R2(K)

262

Q(A2)=S

263	U=SIGAMY(N,M,R3)

264	IF (U.GE.S) GO TO 40

265	130    DO 140 K=1,N

266	140    P(A2,K)=R3(K)

267

Q(A2)=U

268

GO TO 40

269	150    DO 160 K=1,N

270	160    R3(K)=R1(K)/2.+P(A2,K)/2.

271	U=SIGAMY(N,M,R3)

272	IF (U.LT.B1) GO TO 130

273

V=V/2.0

274

GOTO 10

275	170    CONTINUE

276

RETURN

277

END

278

279	FUNCTION SIGAMY(N,M,X)

280	DIMENSION X(5),Y(60)

281	CALL VFRE(X,N,M,Y)

282

A=0.0

283	DO 10 I=1,M

284	10     A=A+Y(I)*Y(I)

285	A=SQRT(A/FLOAT(M))

286	C        WRITE(*,5) A,(X(J),J=1,N)

287	5     FORMAT(2X,'S9=',E12.6,2X,'X(1)=',E12.6,2X,'X(2)=',E12.6,2X,

288	*  'X(3)=',E12.6)

289

SIGAMY=A

290

RETURN

291

END