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