program main !主程序开始
REAL t,p,k,eps,ych4,yh2s,yh2,ys,ycs2,x,G,L !定义实型变量,单精度
external f
common p,L
data t/1000/ !使用data命令赋值
data eps/1e-6/ !精度确定
data G1/41.17/,G2/17.45/,G3/19.38/,R/8.314e-03/ !使用DATA命令批量输入数值
G=G2+G1-G3 !定义数值,反应吉布斯自由能
L=exp(-G/R/t) !定义数值,反应标准平衡常数
p=152000 !赋值
a=0.0 !赋值
b=1.0 !赋值
eps=1.0e-6 !精度确定
call half(a,b,eps,f,x) !调用二分法子程序
ych4=(1-x)/(3-x) !计算数值,平衡组成甲烷含量
yh2s=(x)/(3-x) !计算数值,平衡组成硫化氢含量
ys=(2-3*x)/(3-x) !计算数值,平衡组成硫含量
ycs2=(x)/(3-x) !计算数值,平衡组成二硫化碳含量
yh2=x/(3-x) !计算数值,平衡组成氢气含量
Write(*,*)'T=:'
write(*,*)t !输出
Write(*,*)'P=:'
write(*,*)P !输出
Write(*,*)'x='
write(*,*)x !输出结果
Write(*,*)'ych4=:'
write(*,*)ych4 !输出结果
Write(*,*)'ys=:'
write(*,*)ys !输出结果
Write(*,*)'ycs2=:'
write(*,*)ycs2 !输出结果
Write(*,*)'yh2s=:'
write(*,*)yh2s !输出结果
Write(*,*)'yh2=:'
write(*,*)yh2 !输出结果
end
function f(x) !定义函数
real x,p,L !定义实型变量,单精度
common p,L
f=x*x*x*(3-x)*(1.013e5)/(1-x)/(2-3*x)/(2-3*x)/(2-3*x)/p-L !函数表达式
return
end
SUBROUTINE HALF(A,B,EPS,F,X) !二分法子程序开始
REAL A,B,EPS !定义实型变量,单精度
Y0=F(A) !变量赋值
A1=A
B1=B
K=0
10 X=(A1+B1)*0.5
Y=F(X)
IF(Y*Y0.GT.0.0)A1=X !判断是否异号
IF(Y*Y0.LE.0.0)B1=X !判断是否异号
IF(ABS(B1-A1)/X.GT.EPS)GOTO 10 !迭代次数控制
K=K+1
IF(K.GT.50) THEN
WRITE(*,*)'NO RESOLUTION'
GOTO 30
END IF
X=(A1+B1)*0.5 !输出返回x的值
30 RETURN
END !二分法子程序结束