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