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