10000 REM HGK Vapor-liquid equilibrium.
10050 DEFDBL A-H, M-Z
10100 DIM HGKG(40),II(40),JJ(40),BP(10),BQ(10)
10150 DIM ATZ(4),ADZ(4),AAT(4),AAD(4)
10200 DIM BV(10),A(8),C(18)
10250 DIM QR(11),QT(10),QZR(9),QZT(9)
10300 DIM FFD(2),FFP(5),FFH(5),NNT$(2),NND$(2),NNP$(5),NNH$(5)
10350 GOSUB *BLOCKDATA
10400 GOSUB *UNIT
10450 PRINT
10500 PRINT "Enter temperature"
10550 INPUT"Temperature";TT
10600 T=TT
10650 GOSUB *TTTT
10700 T=TTT
10750 IF T>647.126# THEN PRINT "Input temperature>critical temperature.":GOTO 10550
10800 RT=GASCON*T
10850 GOSUB *BBT
10900 DLL=0:DVV=0:DLIQ=0:DVAP=0
10950 GOSUB *PCORRTPDLDV
11000 PPP=P
11050 D=DL
11060 IF T>646.3# THEN DOUT=D:GOTO 11200
11100 GOSUB *DFINDDOUTPDTDPD
11150 D=DOUT
11200 GOSUB *THERMDT
11250 DD=DOUT/FD
11300 U=UD*RT*FH
11350 H=HD*RT*FH
11400 S=SD*GASCON*FH
11450 CP=CPD*GASCON*FH
11500 CV=CVDX*GASCON*FH
11550 DPDD=DPDD*FD*FP
11600 DPDT=DPDT*FP
11610 G=GD*RT*FH
11620 A=AD*RT*FH
11650 PRES=PPP*FP
11700 LPRINT"Units"
11750 LPRINT A1$;SPC(3);NT$
11800 LPRINT A2$;SPC(3);ND$
11850 LPRINT A3$;SPC(3);NP$
11900 LPRINT A4$;SPC(3);NH$
11950 LPRINT
12000 LPRINT"Liquid phase"
12050 LPRINT USING"T=####.####   P= +#.#####^^^^^   D= +#.#####^^^^^";TT,PRES,DD
12100 LPRINT USING"DP/DT= +#.#####^^^^^   DP/DD= +#.#####^^^^^";DPDT,DPDD
12150 LPRINT USING"CP=    +#.#####^^^^^   CV=    +#.#####^^^^^";CP,CV
12200 LPRINT USING"S=     +#.#####^^^^^   H=     +#.#####^^^^^  U= +#.#####^^^^^";S,H,U
12210 LPRINT USING"G=     +#.#####^^^^^   A=     +#.#####^^^^^";G,A
12300 D=DV
12310 IF T>646.3# THEN DOUT=D:GOTO 12450
12350 GOSUB *DFINDDOUTPDTDPD
12400 D=DOUT
12450 GOSUB *THERMDT
12500 DD=DOUT/FD
12550 U=UD*RT*FH
12600 H=HD*RT*FH
12650 S=SD*GASCON*FH
12700 CP=CPD*GASCON*FH
12750 CV=CVDX*GASCON*FH
12800 DPDD=DPDD*FD*FP
12850 DPDT=DPDT*FP
12860 G=GD*RT*FH
12870 A=AD*RT*FH
12900 LPRINT"Vapor phase"
12950 LPRINT USING"T=####.####   P= +#.#####^^^^^   D= +#.#####^^^^^";TT,PRES,DD
13000 LPRINT USING"DP/DT= +#.#####^^^^^   DP/DD= +#.#####^^^^^";DPDT,DPDD
13050 LPRINT USING"CP=    +#.#####^^^^^   CV=    +#.#####^^^^^";CP,CV
13100 LPRINT USING"S=     +#.#####^^^^^   H=     +#.#####^^^^^  U= +#.#####^^^^^";S,H,U
13110 LPRINT USING"G=     +#.#####^^^^^   A=     +#.#####^^^^^";G,A
13150 LPRINT:LPRINT:LPRINT
13200 INPUT"Will you continue the calculation?  Input Y(or y) or N(or n)";CAL$
13250 IF CAL$="Y" OR CAL$="y" THEN GOTO 10450
13300 END
13350 *BBT
13400 BV(1)=1#
13450 FOR I=2 TO 10
13500 BV(I)=BV(I-1)*TZ/T
13550 NEXT I
13600 B1=BP(1)+BP(2)*LOG(1#/BV(2))
13650 B2=BQ(1)
13700 B1T=BP(2)*BV(2)/TZ
13750 B2T=0
13800 B1TT=0
13850 B2TT=0
13900 FOR I=3 TO 10
13950 B1=B1+BP(I)*BV(I-1)
14000 B2=B2+BQ(I)*BV(I-1)
14050 B1T=B1T-CDBL((I-2))*BP(I)*BV(I-1)/T
14100 B2T=B2T-CDBL((I-2))*BQ(I)*BV(I-1)/T
14150 B1TT=B1TT+BP(I)*CDBL((I-2))*CDBL((I-2))*BV(I-1)/(T*T)
14200 B2TT=B2TT+BQ(I)*CDBL((I-2))*CDBL((I-2))*BV(I-1)/(T*T)
14250 NEXT I
14300 B1TT=B1TT-B1T/T
14350 B2TT=B2TT-B2T/T
14400 RETURN
14450 *BASEDT
14500 Y=.25#*B1*D
14550 XX=1#-Y
14600 Z0=(1#+ALPHA*Y+BETA*Y*Y)/(XX*XX*XX)
14650 Z=Z0+4#*Y*(B2/B1-GAMMA)
14700 DZ0=(ALPHA+2#*BETA*Y)/(XX*XX*XX)+3#*(1#+ALPHA*Y+BETA*Y*Y)/(XX*XX*XX*XX)
14750 DZB=DZ0+4#*(B2/B1-GAMMA)
14800 AB=(-1#)*LOG(XX)-(BETA-1#)/XX+28.16666667#/(XX*XX)+4#*Y*(B2/B1-GAMMA)+15.166666667#+LOG(D*RT/.101325)
14900 BASEF=Z
14950 BB2TT=T*T*B2TT
15000 UB=(-1#)*T*B1T*(Z-1#-D*B2)/B1-D*T*B2T
15100 CVB=2#*UB+(Z0-1#)*((T*B1T/B1)*(T*B1T/B1)-T*T*B1TT/B1)-D*(BB2TT-GAMMA*B1TT*T*T)-(T*B1T/B1)*(T*B1T/B1)*Y*DZ0
15150 DPDTB=BASEF/T+BASEF*D/Z*(DZB*B1T/4#+B2T-B2/B1*B1T)
15200 SB=UB-AB
15250 RETURN
15300 *QQTD
15350 QR(1)=0
15400 Q5=0
15450 Q=0
15500 AR=0
15550 DADT=0
15600 CVR=0
15650 DPDTR=0
15700 E=EXP((-1#)*AA*D)
15750 Q10=D*D*E
15800 Q20=1#-E
15850 QR(2)=Q10
15900 QV=TZ/T
15950 QT(1)=T/TZ
16000 FOR I=2 TO 10
16050 QR(I+1)=QR(I)*Q20
16100 QT(I)=QT(I-1)*QV
16150 NEXT I
16200 FOR I=1 TO INC
16250 K=II(I)+1
16300 L=JJ(I)
16350 QK=CDBL(K) : QL=CDBL(L)
16450 QZR(K-1)=QR(K+1):QZT(L)=QT(L+1):QZR(K)=QR(K+2):QZT(L+1)=QT(L+2)
16500 QP=HGKG(I)*AA*QZR(K-1)*QZT(L)
16550 Q=Q+QP
16600 Q5=Q5+AA*(2#/D-AA*(1#-E*(QK-1#)/Q20))*QP
16650 AR=AR+HGKG(I)*QZR(K)*QZT(L)/(Q10*QK*RT)
16700 DFDT=Q20^QK*(1#-QL)*QZT(L+1)/(TZ*QK)
16750 D2F=QL*DFDT
16800 DPT=DFDT*Q10*AA*QK/Q20
16850 DADT=DADT+HGKG(I)*DFDT
16900 DPDTR=DPDTR+HGKG(I)*DPT
16950 CVR=CVR+HGKG(I)*D2F/GASCON
17000 NEXT I
17050 QP=0
17100 Q2A=0
17150 FOR J=37 TO 40
17200 IF HGKG(J)=0 THEN GOTO 18450
17250 K=II(J)
17300 KM=JJ(J)
17350 QK=CDBL(K) : QKM=CDBL(KM)
17400 DDZ=ADZ(J-36)
17450 DEL=D/DDZ-1#
17500 IF ABS(DEL)<1D-010 THEN DEL=1D-010
17600 EX1=(-1#)*AAD(J-36)*DEL^QK
17650 DEX=EXP(EX1)*DEL^QKM
17700 ATT=AAT(J-36)
17750 TX=ATZ(J-36)
17800 TAU=T/TX-1#
17850 EX2=(-1#)*ATT*TAU*TAU
17900 TEX=EXP(EX2)
17950 Q10=DEX*TEX
18000 QM=QKM/DEL-QK*AAD(J-36)*DEL^(QK-1#)
18050 FCT=QM*D*D*Q10/DDZ
18100 Q5T=FCT*(2#/D+QM/DDZ)-(D/DDZ)*(D/DDZ)*Q10*(QKM/(DEL*DEL)+QK*(QK-1#)*AAD(J-36)*DEL^(QK-2#))
18150 Q5=Q5+Q5T*HGKG(J)
18200 QP=QP+HGKG(J)*FCT
18250 DADT=DADT-2#*HGKG(J)*ATT*TAU*Q10/TX
18300 DPDTR=DPDTR-2#*HGKG(J)*ATT*TAU*FCT/TX
18350 Q2A=Q2A+T*HGKG(J)*(4#*ATT*EX2+2#*ATT)*Q10/(TX*TX)
18400 AR=AR+Q10*HGKG(J)/RT
18450 NEXT J
18500 SR=(-1#)*DADT/GASCON
18550 UR=AR+SR
18600 CVR=CVR+Q2A/GASCON
18650 Q=Q+QP
18700 RETURN
18750 *DFINDDOUTPDTDPD
18800 DD=D
18950 LL=0
19000 LL=LL+1
19050 IF DD=<0 THEN DD=1D-008
19100 IF DD>1.9# THEN DD=1.9#
19150 D=DD
19200 GOSUB *QQTD
19250 Q0=Q
19300 GOSUB *BASEDT
19350 PP=RT*DD*BASEF+Q0
19400 DPD=RT*(Z+Y*DZB)+Q5:DQ=DPD
19450 IF DPD>0 THEN GOTO 19650
19500 IF D>=.2967# THEN DD=DD*1.02
19550 IF D<.2967# THEN DD=DD*.98#
19600 IF LL=<10 GOTO 19000
19650 DPDX=DPD*1.1#
19700 IF DPDX<.1# THEN DPDX=.1#
19750 DP=ABS(1#-PP/PPP)
19800 IF DP<1D-009 THEN GOTO 20200
19850 IF D>.3# AND DP<1D-008 THEN GOTO 20200
19900 IF D>.7# AND DP<1D-007 THEN GOTO 20200
19950 XP=(PPP-PP)/DPDX
20000 IF ABS(XP)>.1# THEN XP=XP*.1#/ABS(XP)
20050 DD=DD+XP
20100 IF DD=<0 THEN DD=1D-008
20150 IF LL=<30 THEN GOTO 19000
20200 DOUT=DD
20250 RETURN
20300 *THERMDT
20350 GOSUB *IDEALT
20400 GOSUB *BASEDT
20450 GOSUB *QQTD
20500 QPQ=Q:QDPQ=Q5
20550 Z=BASEF+QPQ/RT/D
20600 DPDD=RT*(BASEF+Y*DZB)+QDPQ
20650 AD=AB+AR+AI-UREF/T+SREF
20700 GD=AD+Z
20750 UD=UB+UR+UI-UREF/T
20800 DPDT=RT*D*DPDTB+DPDTR
20850 CVDX=CVB+CVR+CVIX
20900 CPD=CVDX+T*DPDT*DPDT/(D*D*DPDD*GASCON)
20950 HD=UD+Z
21000 SD=SB+SR+SI-SREF
21050 RETURN
21100 *PST
21150 IF T>314# THEN GOTO 21350
21200 PL=6.3573118#-8858.843#/T+607.56335*T^(-.6#)
21250 PS=.1#*EXP(PL)
21300 RETURN
21350 TR=T/647.25#
21400 W=ABS(1#-TR)
21450 BPST=0
21500 FOR I=1 TO 8
21550 ZPST=CDBL(I)
21600 BPST=BPST+A(I)*W^((ZPST+1#)/2#)
21650 NEXT I
21700 QPST=BPST/TR
21750 PS=22.093*EXP(QPST)
21800 RETURN
21850 *IDEALT
21900 TIDEAL=T/100
21950 TL=LOG(TIDEAL)
22000 GI=(-1#)*(C(1)/TIDEAL+C(2))*TL
22050 HI=(C(2)+C(1)*(1#-TL)/TIDEAL)
22100 CPI=C(2)-C(1)/TIDEAL
22150 FOR I=3 TO 18
22200 GI=GI-C(I)*TIDEAL^CDBL(I-6)
22250 HI=HI+C(I)*CDBL(I-6)*TIDEAL^CDBL(I-6)
22300 CPI=CPI+C(I)*CDBL(I-6)*CDBL(I-5)*TIDEAL^CDBL(I-6)
22350 NEXT I
22400 AI=GI-1#
22450 UI=HI-1#
22500 CVIX=CPI-1#
22550 SI=UI-AI
22600 RETURN
22650 *CORRTPDLDVDELG
22700 IF T>646.3# THEN GOTO 23600
22750 DLIQ=DLL
22800 IF DLL=<0 THEN DLIQ=1.11#-.0004*T
22850 DLL=DLIQ:D=DLIQ
22900 GOSUB *DFINDDOUTPDTDPD
22950 D=DOUT:DL=DOUT
23000 GOSUB *THERMDT
23050 GL=GD
23100 DVAP=DVV
23150 IF DVV=<0 THEN DVAP=PPP/RT
23200 D=DVAP:DVV=DVAP
23250 GOSUB *DFINDDOUTPDTDPD
23300 IF DOUT<5D-007 THEN DOUT=5D-007
23350 D=DOUT:DV=DOUT
23400 GOSUB *THERMDT
23450 GV=GD
23500 DELG=GL-GV
23550 RETURN
23600 PPP=0
23650 IF T>647.126# THEN RETURN
23700 DELG=0
23750 TAUC=.657128#*(1#-T/647.126#)^.325#
23800 DL=.322#+TAUC
23850 DV=.322#-TAUC
23900 D=DV
23950 GOSUB *BASEDT
24000 GOSUB *QQTD
24100 PPP=RT*DV*BASEF+Q
24150 RETURN
24200 *PCORRTPDLDV
24250 GOSUB *PST
24300 PPP=PS
24350 GOSUB *CORRTPDLDVDELG
24400 DP=0
24450 DP=DELG*RT/(1#/DV-1#/DL)
24500 PPP=PPP+DP
24550 IF ABS(DELG)<.00001# THEN GOTO 24650
24600 DLL=DL:DVV=DV:GOTO 24350
24650 P=PPP
24700 RETURN
24750 *UNIT
24800 PRINT"*******************"
24850 PRINT"* Enter units     *"
24900 PRINT"*******************"
24950 PRINT A1$
25000 PRINT"Choose from 1=deg K, 2=deg C"
25050 INPUT IT
25100 IF IT<1 OR IT>2 THEN GOTO 25000
25150 NT$=NNT$(IT)
25200 PRINT A2$
25250 PRINT"Choose from 1=kg/m3, 2=g/cm3"
25300 INPUT ID
25350 IF ID>2 OR ID<1 THEN GOTO 25250
25400 ND$=NND$(ID)
25450 FD=FFD(ID)
25500 PRINT A3$
25550 PRINT"Choose from 1=MPa, 2=bar"
25600 INPUT IP
25650 IF IP>2 OR IP<1 THEN GOTO 25550
25700 NP$=NNP$(IP)
25750 FP=FFP(IP)
25800 PRINT A4$
25850 PRINT"Choose from 1=kJ/kg, 2=J/g, 3=J/mol"
25900 INPUT IH
25950 IF IH>3 OR IH<1 THEN GOTO 25850
26000 NH$=NNH$(IH)
26050 FH=FFH(IH)
26100 RETURN
26150 *TTTT
26200 ON IT GOTO 26250, 26400
26250 TTT=T
26350 GOTO 26500
26400 TTT=T+273.15#
26500 RETURN
26550 *BLOCKDATA
26600 FOR I=1 TO 4:READ ATZ(I):NEXT I
26650 DATA 640#,640#,641.6#,270#
26700 FOR I=1 TO 4:READ ADZ(I):NEXT I
26750 DATA 0.319#,0.319#,0.319#,1.55#
26800 FOR I=1 TO 4:READ AAT(I):NEXT I
26850 DATA 2.0D+004,2.0D+004,4.0D+004,25.0#
26900 FOR I=1 TO 4:READ AAD(I):NEXT I
26950 DATA 34.0#,40.0#,30.0#,1.05D+003
27000 WM=18.0152:GASCON=.461522#:TZ=647.073:AA=1#:INC=36
27050 UREF=-4328.454977#:SREF=7.6180720#
27100 ALPHA=11#:BETA=44.333333333333#:GAMMA=3.5#
27150 FOR I=1 TO 10:READ BP(I):NEXT I
27200 DATA 0.7478629#,-0.3540782#,0.0#,0.0#,0.007159876#,0.0#,-0.003528426#,0.0#,0.0#,0.0#
27250 FOR I=1 TO 10:READ BQ(I):NEXT I
27300 DATA 1.1278334#,0.0#,-0.5944001#,-5.010996#,0.0#,0.63684256#,0.0#,0.0#,0.0#,0.0#
27350 FOR I=1 TO 40:READ HGKG(I):NEXT I
27400 DATA -5.3062968529023D+002,2.2744901424408D+003,7.8779333020687D+002
27450 DATA -6.9830527374994D+001,1.7863832875422D+004,-3.9514731563338D+004
27500 DATA 3.3803884280753D+004,-1.3855050202703D+004,-2.5637436613260D+005
27550 DATA 4.8212575981415D+005,-3.4183016969660D+005,1.2223156417448D+005
27600 DATA 1.1797433655832D+006,-2.1734810110373D+006,1.0829952168620D+006
27650 DATA -2.5441998064049D+005,-3.1377774947767D+006,5.2911910757704D+006
27700 DATA -1.3802577177877D+006,-2.5109914369001D+005,4.6561826115608D+006
27750 DATA -7.2752773275387D+006,4.1774246148294D+005,1.4016358244614D+006
27800 DATA -3.1555231392127D+006,4.7929666384584D+006,4.0912664781209D+005
27850 DATA -1.3626369388386D+006,6.9625220862664D+005,-1.0834900096447D+006
27900 DATA -2.2722827401688D+005,3.8365486000660D+005,6.8833257944332D+003
27950 DATA 2.1757245522644D+004,-2.6627944829770D+003,-7.0730418082074D+004
28000 DATA -0.225#,-1.68#,0.055#,-93.0#
28050 FOR I=1 TO 40:READ II(I):NEXT I
28100 DATA 0,0,0,0,1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,5,6,6,6,6,8,8,8,8,2,2,0,4,2,2,2,4
28150 FOR I=1 TO 40:READ JJ(I):NEXT I
28200 DATA 2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,2,3,5,7,1,4,4,4,0,2,0,0
28250 FOR I=1 TO 8:READ A(I):NEXT I
28300 DATA -7.8889166#,2.5514255#,-6.716169#,33.239495#
28350 DATA -105.38479#,174.35319#,-148.39348#,48.631602#
28400 FOR I=1 TO 18:READ C(I):NEXT I
28450 DATA 1.9730271018D+001,2.09662681977D+001,-4.83429455355D-001,6.05743189245D+000
28500 DATA 2.256023885D+001,-9.87532442D+000,-4.3135538513D+000,4.58155781D-001
28550 DATA -4.7754901883D-002,4.1238460633D-003,-2.7929052852D-004
28600 DATA 1.4481695261D-005,-5.6473658748D-007,1.6200446D-008,-3.303822796D-010
28650 DATA 4.51916067368D-012,-3.70734122708D-014,1.37546068238D-016
28700 FOR I=1 TO 2:READ FFD(I):NEXT I
28750 DATA 1.0D-003,1.0#
28800 FOR I=1 TO 2:READ FFP(I):NEXT I
28850 DATA 1.0#,10.0#
28900 FOR I=1 TO 3:READ FFH(I):NEXT I
28950 DATA 1.0#,1.0#,18.0152#
29000 FOR I=1 TO 2:READ NNT$(I):NEXT I
29050 DATA " K"," C"
29100 FOR I=1 TO 2:READ NND$(I):NEXT I
29150 DATA " kg/m3"," g/cm3"
29200 FOR I=1 TO 2:READ NNP$(I):NEXT I
29250 DATA "   MPa","   bar"
29300 FOR I=1 TO 3:READ NNH$(I):NEXT I
29350 DATA " kJ/kg","   J/g"," J/mol"
29400 A1$="TEMPERATURE":A2$="DENSITY":A3$="PRESSURE":A4$="ENERGY"
29450 RETURN