Attribute VB_Name = "HeProp"
Option Explicit
Option Base 1

Private Const NCF As Integer = 76
Private Const NOV As Integer = 39
Private Const NIV As Integer = 10
Private Const NIP As Integer = 15
Private Const gintPrecis As Integer = 3
Private E(0 To 15) As Double

'The following parameters are originally from Common SubRT
Private SubRT(13) As Double
Private JR As Integer
            
Private Const UnivGasConst As Double = 0.0083143   'kJ/g-mole-K
Private Const HeIEoSGamma  As Double = -0.0033033259
Private HeIEoS(32)     As Double
Private HeIIEoS(41, 2) As Double
Private Prop(0 To 41, 0 To 2) As Double
Private Kunits(0 To 40) As Integer
Private Units(NCF) As String
Private ConvF(NCF) As Double
Private gstrLabels(NOV + 6) As String
Private gintPropToInput(0 To NOV + 6) As Integer
Private gstrInputLabels(NIP) As String
Private gintLabelToProp(NOV + 6) As Integer
Public gstrPairs As String

'The following parameters are originally from Common PRInit
Private Const HeGasConst   As Double = 2077.23476   'J/kg-K
Private Const GMWT   As Double = 4.0026             'g/g-mole
Private Const Pcrit  As Double = 227462.3           'Pa
Private Const Tcrit  As Double = 5.1953             'K
Private Const Dcrit  As Double = 69.6412374         'kg/m3
Private Const Hcrit  As Double = 21932#             'J/kg
Private Const Ucrit  As Double = 18682#             'J/kg
Private Const Scrit  As Double = 5768.6             'J/kg-K
Private Const PLamL  As Double = 5041.8             'Pa       PTR
Private Const TLamL  As Double = 2.1768             'K        TTR, T0
Private Const DLamLf As Double = 146.15             'kg/m3    DTTR
Private Const DLamLg As Double = 1.1923727          'kg/m3    DTRG
Private Const HLamLf As Double = 2961.33            '/kg      HTRF
Private Const HLamLg As Double = 25757.05           'J/kg     HTRG
Private Const ULamLf As Double = 2926.83            'J/kg     UTRF
Private Const ULamLg As Double = 21528.67           'J/kg     UTRG
Private Const SLamLf As Double = 1579.97            'J/kg-K   STRF
Private Const SLamLg As Double = 11991.95           'J/kg-K   STRG
Private Const THI    As Double = 300#               'K
Private Const DHI    As Double = 0.1604             'kg/m3
Private Const HHI    As Double = 1573500#           'J/kg
Private Const UHI    As Double = 950100#            'J/kg
Private Const SHI    As Double = 31612#             'J/kg-K
Private Const CpHI   As Double = 5193#              'J/kg-K
Private Const CvHI   As Double = 3117#              'J/kg-K
Private Const Hmax2P As Double = 30744#             'J/kg, max two-phase H
Private Const Umax2P As Double = 24724#             'J/kg, max two-phase U
Private Const Pmax   As Double = 2028000000#        'Pa
Private Const Pmin   As Double = 0.1                'Pa
Private Const Tmax   As Double = 1500#              'K
Private Const Tmin   As Double = 0.8                'K

Private Const VLamLf As Double = 1 / DLamLf * 1000  'cc/g     V0
Private Const Pref As Double = 0.101325             'MPa      P0

Private Const TNRC As Double = 5.106                'K
Private Const Tlow As Double = 4.804                'K
Private Const DelT = TNRC - Tlow                    'K

'The following parameters are originally from Common SubLam
Private TLam As Double                              'K
Private dTLdV As Double                             'K-gm/cm3
Private d2TLdV2 As Double                           'K-gm2/cm6
Private Vsave As Double

Private Const Zero   As Double = 0#
Private Const One3RD As Double = 1# / 3#
Private Const Two3RD As Double = 2# / 3#
Private Const Two7TH As Double = 2# / 7#

Private Const gstrLicensee As String = _
   "Mr. Cecile Guillot" & vbCrLf & _
   "Air Liquide - Sassenage, FRANCE" & vbCrLf
Private Const gstrCryodata As String = _
   "Cryodata, Inc." & vbCrLf & _
   "P.O. Box 173" & vbCrLf & _
   "Louisville, CO 80027-0173" & vbCrLf & _
   "Phone/Fax  303.664.1354" & vbCrLf & _
   "e-mail  partners@cryodata.com" & vbCrLf & _
   "www.cryodata.com/partners" & vbCrLf
Private Const gstrCopyright As String = _
   "Copyright (C) Cryodata Inc."
Private Const gdblVersion As String = 3.4


'Private Sub Class_Initialize()
'' Used only in class implementation
'   Call Initialize
'End Sub


Static Sub Initialize()
   MsgBox "HePak  Version " & gdblVersion & vbCrLf & vbCrLf & _
          gstrCryodata & vbCrLf & _
          "Licensed to: " & vbCrLf & _
          gstrLicensee & vbCrLf & _
          gstrCopyright, vbOKOnly, "HePak"
   HeIEoS(1) = 4.558980227431E-05
   HeIEoS(2) = 0.001260692007853
   HeIEoS(3) = -0.007139657549318
   HeIEoS(4) = 0.009728903861441
   HeIEoS(5) = -0.01589302471562
   HeIEoS(6) = 1.454229259623E-06
   HeIEoS(7) = -4.708238429298E-05
   HeIEoS(8) = 0.001132915223587
   HeIEoS(9) = 0.002410763742104
   HeIEoS(10) = -5.093547838381E-09
   HeIEoS(11) = 2.6997269279E-06
   HeIEoS(12) = -3.954146691114E-05
   HeIEoS(13) = 1.551961438127E-09
   HeIEoS(14) = 1.050712335785E-08
   HeIEoS(15) = -5.50115836675E-08
   HeIEoS(16) = -1.037673478521E-10
   HeIEoS(17) = 6.446881346448E-13
   HeIEoS(18) = 3.298960057071E-11
   HeIEoS(19) = -3.555585738784E-13
   HeIEoS(20) = -0.00688540136769
   HeIEoS(21) = 0.009166109232806
   HeIEoS(22) = -6.544314242937E-06
   HeIEoS(23) = -3.315398880031E-05
   HeIEoS(24) = -2.067693644676E-08
   HeIEoS(25) = 3.850153114958E-08
   HeIEoS(26) = -1.399040626999E-11
   HeIEoS(27) = -1.888462892389E-12
   HeIEoS(28) = -4.595138561035E-15
   HeIEoS(29) = 6.872567403738E-15
   HeIEoS(30) = -6.097223119177E-19
   HeIEoS(31) = -7.636186157005E-18
   HeIEoS(32) = 3.848665703556E-18
   
   HeIIEoS(1, 1) = -0.55524231:     HeIIEoS(1, 2) = -1.4629809
   HeIIEoS(2, 1) = 0.18333872:      HeIIEoS(2, 2) = 0.76365652
   HeIIEoS(3, 1) = -0.042819388:    HeIIEoS(3, 2) = -0.11943389
   HeIIEoS(4, 1) = -0.049962336:    HeIIEoS(4, 2) = 0.0030525707
   HeIIEoS(5, 1) = -0.083818656:    HeIIEoS(5, 2) = -0.10405828
   HeIIEoS(6, 1) = -0.14081863:     HeIIEoS(6, 2) = -0.14748127
   HeIIEoS(7, 1) = 0.899013:        HeIIEoS(7, 2) = -0.96308013
   HeIIEoS(8, 1) = 6.6841845:       HeIIEoS(8, 2) = 0.67943869
   HeIIEoS(9, 1) = 9.8899347:       HeIIEoS(9, 2) = -0.5527632
   HeIIEoS(10, 1) = 7.3876336:      HeIIEoS(10, 2) = 0.019535989
   HeIIEoS(11, 1) = 2.0130513:      HeIIEoS(11, 2) = -0.04792884
   HeIIEoS(12, 1) = -0.025251119:   HeIIEoS(12, 2) = 0.075410775
   HeIIEoS(13, 1) = -1.0649342:     HeIIEoS(13, 2) = -0.08060907
   HeIIEoS(14, 1) = -3.5520547:     HeIIEoS(14, 2) = -0.766418
   HeIIEoS(15, 1) = -5.846516:      HeIIEoS(15, 2) = -3.5260824
   HeIIEoS(16, 1) = -4.2352097:     HeIIEoS(16, 2) = -5.1006607
   HeIIEoS(17, 1) = -1.2206228:     HeIIEoS(17, 2) = -2.0845765
   HeIIEoS(18, 1) = 0.016905918:    HeIIEoS(18, 2) = -0.012991643
   HeIIEoS(19, 1) = 0.21835833:     HeIIEoS(19, 2) = -0.014392344
   HeIIEoS(20, 1) = 0.74548499:     HeIIEoS(20, 2) = 0.72407645
   HeIIEoS(21, 1) = 1.1932777:      HeIIEoS(21, 2) = 3.5487925
   HeIIEoS(22, 1) = 0.86121784:     HeIIEoS(22, 2) = 4.9508203
   HeIIEoS(23, 1) = 0.25519882:     HeIIEoS(23, 2) = 2.0014366
   HeIIEoS(24, 1) = -0.0013224602:  HeIIEoS(24, 2) = 0.0012086022
   HeIIEoS(25, 1) = -0.020161157:   HeIIEoS(25, 2) = 0.01089437
   HeIIEoS(26, 1) = -0.065799115:   HeIIEoS(26, 2) = -0.21841979
   HeIIEoS(27, 1) = -0.10330654:    HeIIEoS(27, 2) = -1.0697343
   HeIIEoS(28, 1) = -0.075388298:   HeIIEoS(28, 2) = -1.454107
   HeIIEoS(29, 1) = -0.023788871:   HeIIEoS(29, 2) = -0.56844527
   HeIIEoS(30, 1) = 0.000052810077: HeIIEoS(30, 2) = -0.000082701229
   HeIIEoS(31, 1) = 0.00066823412:  HeIIEoS(31, 2) = -0.0010295381
   HeIIEoS(32, 1) = 0.0021297216:   HeIIEoS(32, 2) = 0.02077704
   HeIIEoS(33, 1) = 0.0032541554:   HeIIEoS(33, 2) = 0.10085164
   HeIIEoS(34, 1) = 0.002442111:    HeIIEoS(34, 2) = 0.13267276
   HeIIEoS(35, 1) = 0.00082971274:  HeIIEoS(35, 2) = 0.04970154
   HeIIEoS(36, 1) = -1.9221178:     HeIIEoS(36, 2) = 0.95829687
   HeIIEoS(37, 1) = 1.1203003:      HeIIEoS(37, 2) = -1.4025436
   HeIIEoS(38, 1) = -0.16430199:    HeIIEoS(38, 2) = 0.63551829
   HeIIEoS(39, 1) = -0.0034805651:  HeIIEoS(39, 2) = -0.055978553
   HeIIEoS(40, 1) = 3.63345:        HeIIEoS(40, 2) = -0.044009
   HeIIEoS(41, 1) = -4.325:         HeIIEoS(41, 2) = -0.78777
    
    Units(1) = " [-] ":        ConvF(1) = 1#
    Units(2) = " [kPa] ":      ConvF(2) = 1000#
    Units(3) = " [MPa] ":      ConvF(3) = 1000000#
    Units(4) = " [bar] ":      ConvF(4) = 100000#
    Units(5) = " [atmos] ":    ConvF(5) = 101325#
    Units(6) = " [psi] ":      ConvF(6) = 6894.757
    Units(7) = " [1/kPa] ":    ConvF(7) = 0.001
    Units(8) = " [1/MPa] ":    ConvF(8) = 0.000001
    Units(9) = " [1/bar] ":    ConvF(9) = 0.00001
   Units(10) = " [1/atmos] ":  ConvF(10) = 0.00000986923
   Units(11) = " [1/psi] ":    ConvF(11) = 0.000145037
   Units(12) = " [K] ":        ConvF(12) = 1#
   Units(13) = " [R] ":        ConvF(13) = 0.555555556
   Units(14) = " [1/K] ":      ConvF(14) = 1#
   Units(15) = " [1/F] ":      ConvF(15) = 1.8
   Units(16) = " [kg/m3] ":    ConvF(16) = 1#
   Units(17) = " [g/cm3] ":    ConvF(17) = 1000#
   Units(18) = " [mol/L] ":    ConvF(18) = GMWT
   Units(19) = " [lb/ft3] ":   ConvF(19) = 16.01846
   Units(20) = " [m3/kg] ":    ConvF(20) = 1#
   Units(21) = " [cm3/g] ":    ConvF(21) = 0.001
   Units(22) = " [L/mol] ":    ConvF(22) = 0.2498376
   Units(23) = " [ft3/lb] ":   ConvF(23) = 0.0624279
   Units(24) = " [J/kg-K] ":   ConvF(24) = 1#
   Units(25) = " [J/g-K] ":    ConvF(25) = 1000#
   Units(26) = " [J/mol-K] ":  ConvF(26) = 249.8376
   Units(27) = " [BTU/lb-R] ": ConvF(27) = 4186.8
   Units(28) = " [J/kg] ":     ConvF(28) = 1#
   Units(29) = " [J/g] ":      ConvF(29) = 1000#
   Units(30) = " [J/mol] ":    ConvF(30) = 249.8376
   Units(31) = " [BTU/lb] ":   ConvF(31) = 2326#
   Units(32) = " [m/s] ":      ConvF(32) = 1#
   Units(33) = " [ft/s] ":     ConvF(33) = 0.3048
   Units(34) = " [K/kPa] ":    ConvF(34) = 0.001
   Units(35) = " [K/bar] ":    ConvF(35) = 0.00001
   Units(36) = " [K/Pa] ":     ConvF(36) = 1#
   Units(37) = " [R/psi] ":    ConvF(37) = 0.0000805765
   Units(38) = " [W/m-K] ":    ConvF(38) = 1#
   Units(39) = " BTU/hrftR ":  ConvF(39) = 1.730735
   Units(40) = " [uPa-s] ":    ConvF(40) = 0.000001
   Units(41) = " [upoise] ":   ConvF(41) = 0.0000001
   Units(42) = " [lbm/ft-s] ": ConvF(42) = 1.488164
   Units(43) = " [m2/s] ":     ConvF(43) = 1#
   Units(44) = " [cm2/s] ":    ConvF(44) = 0.0001
   Units(45) = " [ft2/s] ":    ConvF(45) = 0.09290304
   Units(46) = " [N/m] ":      ConvF(46) = 1#
   Units(47) = " [dyne/cm] ":  ConvF(47) = 0.001
   Units(48) = " [lbF/ft] ":   ConvF(48) = 14.5939
   Units(49) = " [m-s/kg] ":   ConvF(49) = 1#
   Units(50) = " [W3/m5-K] ":  ConvF(50) = 1#
   Units(51) = " [cm-s/g] ":   ConvF(51) = 10#
   Units(52) = " [W3/cm5-K] ": ConvF(52) = 10000000000#
   Units(53) = " [1/Pa] ":     ConvF(53) = 1#
   Units(54) = " [V/Vc] ":     ConvF(54) = 1# / Dcrit
   Units(55) = " [F] ":        ConvF(55) = 0.555555556
   Units(56) = " [C] ":        ConvF(56) = 1#
   Units(57) = " [Pa] ":       ConvF(57) = 1#
   Units(58) = " [Pa-s] ":     ConvF(58) = 1#
   Units(59) = " [P/Pc] ":     ConvF(59) = 227462.3
   Units(60) = " [T/Tc] ":     ConvF(60) = Tcrit
   Units(61) = " [D/Dc] ":     ConvF(61) = Dcrit
   Units(62) = " [TdP/PdT] ":  ConvF(62) = 1#
   Units(63) = " [../R] ":     ConvF(63) = HeGasConst
   Units(64) = " [../RTc] ":   ConvF(64) = HeGasConst * Tcrit
   Units(65) = " [TdV/VdT] ":  ConvF(65) = 1#
   Units(66) = " [PdD/DdP] ":  ConvF(66) = 1#
   Units(67) = " [PdT/TdP] ":  ConvF(67) = 1#
   Units(68) = " [psi/R] ":    ConvF(68) = 12410.56
   Units(69) = " [Pa/K] ":     ConvF(69) = 1#
   Units(70) = " [kPa/K] ":    ConvF(70) = 1000#
   Units(71) = " [MPa/K] ":    ConvF(71) = 1000000#
   Units(72) = " [bar/K] ":    ConvF(72) = 100000#
   Units(73) = " [Pa-m3/kg] ": ConvF(73) = 1#
   Units(74) = " [kPaL/mol] ": ConvF(74) = 249.8376
   Units(75) = " [kPam3/kg] ": ConvF(75) = 1000#
   Units(76) = " psi-ft3/lb ": ConvF(76) = 430.4257
    
   'gstrLabels is an array of property labels.
   'gintLabelToProp maps the labels in gstrLabels() to the Prop() index.
   'Added last six elements for VBA implementation.
   'gintLabelToProp(label#) = Prop() index #
   gstrLabels(1) = " Pressure ":           gintLabelToProp(1) = 1
   gstrLabels(2) = " Temperature ":        gintLabelToProp(2) = 2
   gstrLabels(3) = " Density ":            gintLabelToProp(3) = 3
   gstrLabels(4) = " Specific Volume ":    gintLabelToProp(4) = 4
   gstrLabels(5) = " Enthalpy ":           gintLabelToProp(5) = 9
   gstrLabels(6) = " Entropy ":            gintLabelToProp(6) = 8
   gstrLabels(7) = " Internal Energy ":    gintLabelToProp(7) = 11
   gstrLabels(8) = " Quality ":            gintLabelToProp(8) = 0
   gstrLabels(9) = " Latent Heat ":        gintLabelToProp(9) = 7
   gstrLabels(10) = " DPTSat ":            gintLabelToProp(10) = 6
   gstrLabels(11) = " Z = PV/RT ":         gintLabelToProp(11) = 5
   gstrLabels(12) = " Cp ":                gintLabelToProp(12) = 14
   gstrLabels(13) = " Cv ":                gintLabelToProp(13) = 15
   gstrLabels(14) = " Gamma ":             gintLabelToProp(14) = 16
   gstrLabels(15) = " Expansivity ":       gintLabelToProp(15) = 17
   gstrLabels(16) = " Gruneisen ":         gintLabelToProp(16) = 18
   gstrLabels(17) = " Compressibility ":   gintLabelToProp(17) = 19
   gstrLabels(18) = " Sound Velocity ":    gintLabelToProp(18) = 20
   gstrLabels(19) = " JT Coefficient ":    gintLabelToProp(19) = 21
   gstrLabels(20) = " V*dHdV|P ":          gintLabelToProp(20) = 24
   gstrLabels(21) = " dPdD|T ":            gintLabelToProp(21) = 22
   gstrLabels(22) = " dPdT|D ":            gintLabelToProp(22) = 23
   gstrLabels(23) = " Conductivity ":      gintLabelToProp(23) = 26
   gstrLabels(24) = " Viscosity ":         gintLabelToProp(24) = 25
   gstrLabels(25) = " Prandtl # ":         gintLabelToProp(25) = 27
   gstrLabels(26) = " Thermal Diff ":      gintLabelToProp(26) = 28
   gstrLabels(27) = " Surface Tension ":   gintLabelToProp(27) = 29
   gstrLabels(28) = " Dielectric - 1 ":    gintLabelToProp(28) = 30
   gstrLabels(29) = " Refraction - 1 ":    gintLabelToProp(29) = 31
   gstrLabels(30) = " RhoS/Rho ":          gintLabelToProp(30) = 34
   gstrLabels(31) = " 2nd Sound Vel. ":    gintLabelToProp(31) = 35
   gstrLabels(32) = " 4th Sound Vel. ":    gintLabelToProp(32) = 36
   gstrLabels(33) = " Gorter-Mellink ":    gintLabelToProp(33) = 37
   gstrLabels(34) = " SFTC ":              gintLabelToProp(34) = 38
   gstrLabels(35) = " Res. for Fugacity ": gintLabelToProp(35) = 13
   gstrLabels(36) = " Gibbs Energy ":      gintLabelToProp(36) = 12
   gstrLabels(37) = " dT(V) ":             gintLabelToProp(37) = 33
   gstrLabels(38) = " dT(P) ":             gintLabelToProp(38) = 32
   gstrLabels(39) = " Helmholtz ":         gintLabelToProp(39) = 10
   gstrLabels(40) = " (T-Tlambda) ":       gintLabelToProp(40) = 39
   gstrLabels(41) = " Lambda line ":       gintLabelToProp(41) = 41
   gstrLabels(42) = " Melting line ":      gintLabelToProp(42) = 42
   gstrLabels(43) = " Saturated liquid ":  gintLabelToProp(43) = 43
   gstrLabels(44) = " Saturated vapor ":   gintLabelToProp(44) = 44
   gstrLabels(45) = " Liquid and Vapor ":  gintLabelToProp(45) = 0

   gintPropToInput(0) = 9
   gintPropToInput(1) = 1
   gintPropToInput(2) = 2
   gintPropToInput(3) = 3
   gintPropToInput(4) = 4
   gintPropToInput(5) = -2
   gintPropToInput(6) = -2
   gintPropToInput(7) = -2
   gintPropToInput(8) = 5
   gintPropToInput(9) = 6
   gintPropToInput(10) = -2
   gintPropToInput(11) = 7
   gintPropToInput(12) = 8
   gintPropToInput(13) = -2
   gintPropToInput(14) = -2
   gintPropToInput(15) = -2
   gintPropToInput(16) = -2
   gintPropToInput(17) = -2
   gintPropToInput(18) = -2
   gintPropToInput(19) = -2
   gintPropToInput(20) = -2
   gintPropToInput(21) = -2
   gintPropToInput(22) = -2
   gintPropToInput(23) = -2
   gintPropToInput(24) = -2
   gintPropToInput(25) = -2
   gintPropToInput(26) = -2
   gintPropToInput(27) = -2
   gintPropToInput(28) = -2
   gintPropToInput(29) = -2
   gintPropToInput(30) = -2
   gintPropToInput(31) = -2
   gintPropToInput(32) = 10
   gintPropToInput(33) = 10
   gintPropToInput(34) = -2
   gintPropToInput(35) = -2
   gintPropToInput(36) = -2
   gintPropToInput(37) = -2
   gintPropToInput(38) = -2
   gintPropToInput(39) = 10
   gintPropToInput(40) = -2
   gintPropToInput(41) = 15
   gintPropToInput(42) = 11
   gintPropToInput(43) = 13
   gintPropToInput(44) = 14
   gintPropToInput(45) = -2

   gstrPairs = _
   "Array of valid pairs of input parameters (marked x)" & vbCrLf & _
   "     P  T  D  V  S  H  U  G  X dT  M SL SV  L" & vbCrLf & _
   "     1  2  3  4  8  9 11 12  0 39 42 43 44 41" & vbCrLf & _
   "     ----------------------------------------" & vbCrLf & _
   " P|     x  x  x  x  x  x  x  x  x  x  x  x  x" & vbCrLf & _
   " T|  x     x  x  x        x  x     x  x  x  x" & vbCrLf & _
   " D|  x  x        x  x  x        x     x  x  x" & vbCrLf & _
   " V|  x  x        x  x  x        x     x  x  x" & vbCrLf & _
   " S|  x  x  x  x     x                 x  x   " & vbCrLf & _
   " H|  x     x  x  x                    x      " & vbCrLf & _
   " U|  x     x  x                       x      " & vbCrLf & _
   " G|  x  x                             x      " & vbCrLf & _
   " X|  x  x                                    " & vbCrLf & _
   "dT|  x     x  x                       x      " & vbCrLf & _
   " M|  x  x                                   x" & vbCrLf & _
   "SL|  x  x  x  x  x  x  x  x     x           x" & vbCrLf & _
   "SV|  x  x  x  x  x                           " & vbCrLf & _
   " L|  x  x  x  x                    x  x      "
   Call Prcis(gintPrecis)
End Sub


Private Function RootFunc(I As Integer, X As Double) As Double
   Select Case I
      Case 1
         RootFunc = TLfD(X)
      Case 2
         RootFunc = PLfT(X)
      Case 3
         RootFunc = DMfT(X)
      Case 4
         RootFunc = PMfT(X)
      Case 5
         RootFunc = Psat(X)
      Case 6
         RootFunc = SatS(X)
      Case 7
         RootFunc = SatD1(X)
      Case 8
         RootFunc = SatD(X)
      Case 9
         RootFunc = SatLY(X)
      Case 10
         RootFunc = SatVY(X)
      Case 11
         RootFunc = DPT(X)
      Case 12
         RootFunc = DST(X)
      Case 13
         RootFunc = DTHS(X)
      Case 14
         RootFunc = TDP(X)
      Case 15
         RootFunc = YJfDT(X)
      Case 16
         RootFunc = DTG(X)
      Case Else
         RootFunc = -1
      End Select
End Function
      
      
Private Function Min(First As Double, Second As Double) As Double
   Min = First
   If First > Second Then Min = Second
End Function


Private Function Max(First As Double, Second As Double) As Double
   Max = First
   If First < Second Then Max = Second
End Function


Private Function IMin(First As Integer, Second As Integer) As Integer
   IMin = First
   If First > Second Then IMin = Second
End Function


Private Function IMax(First As Integer, Second As Integer) As Integer
   IMax = First
   If First < Second Then IMax = Second
End Function


Private Function Sign(A As Double, B As Double) As Double
'VB version of the FORTRAN function
   If B >= 0# Then
      Sign = Abs(A)
   Else
      Sign = -Abs(A)
   End If
End Function


Public Function HeConst(Index As Integer) As Double
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Returns a parameter associated with the state equation, as determined
' by the Index parameter.
' See global declaration section for values and units.
   Select Case Index
      Case 0                   'User wants the error status code,
         HeConst = 0#          'but at this point we can have no error.
      Case 1
         HeConst = Pmax        'Maximum valid pressure in HePak
      Case 2
         HeConst = Tmax        'Maximum valid temperature in HePak
      Case 3
         HeConst = Tmin        'Minimum valid temperature in HePak
      Case 4
         HeConst = HeGasConst  'Helium gas constant
      Case 5
         HeConst = Pcrit       'Critical pressure
      Case 6
         HeConst = Tcrit       'Critical temperature
      Case 7
         HeConst = Dcrit       'Critical density
      Case 8
         HeConst = PLamL       'Lambda point pressure
      Case 9
         HeConst = TLamL       'Lambda point temperature
      Case 10
         HeConst = DLamLg      'Vapor density at the lambda point
      Case 11
         HeConst = DLamLf      'Liquid density at the lambda point
      Case 12
         HeConst = GMWT        'Molecular weight
      Case 13
         HeConst = Pref        'Reference pressure for entropy scale
      Case 14
         HeConst = 0#          'reserved
      Case 15
         HeConst = gdblVersion 'HePak properties module version #
      Case Else
         HeConst = 0#          'Index out of range
   End Select
End Function


Public Function HeCalc(Index As Integer, Phase As Integer, _
                       Input1 As Variant, Value1 As Double, _
                       Input2 As Variant, Value2 As Double, _
                       Unit As Integer) As Variant
' (C) Copyright (1993), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' HeProp function, calls HeProp's Calc subroutine.
' Valid for HeProp v3.23 (and later);
' This function contains no code specific to spreadsheets.
'
' The calling parameter Index specifies the HeProp property function;
'     valid range is 0 to 39.  See USER.FOR for definitions.
'
' The calling parameter Phase combines with the HeProp's output
'     parameter IDID to specify the K in the output properties array
'     Prop(Index,K) from which HeCalc is evaluated.  Some combinations
'     of (positive) IDID and Phase force HeCalc to an invalid answer
'     (see function JPCalc).
'
' Input1 and Input2 specifies the input property type and has been
' changed to be the consistent with the Prop() array index, with the
' exception that 0 to 44 is supported to add the melting line, lambda
' line and saturated liquid and vapor. The alphabetic representation
' is still supported ("P", "sL" etc.).
'
' Version Dec. 08, 1998
'
' Set precision = 2; not controlled by the calling program.
'     This choice can be changed (1, 2, or 3, as desired).
'
   Const Four7 As Double = -400# / 7#
   Const Five7 As Double = -500# / 7#
   Const JAll As Integer = 11111
   Const JRefr As Integer = 10
   Dim Repeat As Boolean
   Dim IDID As Integer
   Dim J As Integer, JP As Integer
   Dim In1 As Integer, In2 As Integer
   Static In1S As Integer
   Static Valu1S As Double
   Static In2S As Integer
   Static Valu2S As Double
   Static JUnits As Integer
   Static IDIDS As Integer
   Static WaveSv As Double
   Dim WaveLn As Double
'
' Error codes:
' IDID = output indicator from HeProp's Calc; will be < 0 if
'        Calc cannot do the requested calculation.
' IErr = error indicator from HeCalc; = 0 for no error;
'        otherwise, = 36 if IDID < 0 or if Index, Phase, or Unit
'        are out of range
'        otherwise will = 42 if a phase error occurs (see JPCalc)
'
   HeCalc = 0#
'
   In1 = UserToIn(Input1) 'Convert to original (1-15) integer input
   In2 = UserToIn(Input2) 'Convert to original (1-15) integer input
   If (((Index + 1) * (Index - 39) > 0) Or _
       (Unit * (Unit - 5) >= 0) Or _
       (Phase * (Phase - 5) > 0)) Then
' Index specifies the desired property, i.e.,
' answer = Prop(Index, k)
' Index = -1 requests HeCalc to return the error code (IDID)
' Phase controls the second Prop Index, i.e., k in Prop(i, k)
      IDIDS = -13
      HeCalc = "NUM" 'IErr = 36
      Valu2S = Five7
      Exit Function
   End If
'
' Is this a repeat of a previous data point?
'
   Repeat = False
   If (Value2 = Valu2S) Then
      If (Value1 = Valu1S) Then
         If ((In1 = In1S) And (In2 = In2S) And _
             (Unit = JUnits) And _
             (Prop(1, 0) + Prop(1, 1) + Prop(1, 2)) >= 0) And _
             IDIDS > 0 Then Repeat = True
      End If
   End If
'
' Note that Repeat can be true even if Phase changes.
'
' Special case: If the wavelength changes, and refractive Index is
'     requested, but otherwise the calculation is a repeat,
'     then a unique call to Calc is all that is required.
'
   If Prop(40, 0) <= 0# Then Prop(40, 0) = 0.55 'default wavelength
   WaveLn = Prop(40, 0)
   If (Repeat And (Index = 31) And (WaveLn <> WaveSv)) Then
      WaveSv = WaveLn
      Prop(41, 2) = -262144#
      Call Calc(IDIDS, Prop, JRefr, _
                In1S, Valu1S, In2S, Valu2S, _
                gintPrecis, JUnits)
      Prop(41, 2) = 0#
      JP = JPCalc(Phase, IDIDS)
      If (JP <> -42) Then
         HeCalc = Prop(Index, JP)
      Else
         HeCalc = "N/A" 'IErr = 42 IDID -301
      End If
      Exit Function
   End If
'
   If (Not Repeat) Then
      Call Calc(IDID, Prop, JAll, _
                In1, Value1, In2, Value2, _
                gintPrecis, Unit)
      In1S = In1
      In2S = In2
      Valu1S = Value1
      Valu2S = Value2
      JUnits = Unit
      WaveSv = WaveLn
      IDIDS = IDID
   End If
'
' Obtain requested output
'
   JP = JPCalc(Phase, IDIDS)
   If (Index = -1) Then
      HeCalc = CDbl(IDIDS)
      Exit Function
   End If
   If ((IDIDS >= 0) And (JP < 0)) Then
      HeCalc = -15#
   ElseIf (IDIDS < 0) Then
'
' When IDID < 0, HeProp found a problem.
' Either the indices In1 and/or In2 were invalid
' (-99 < IDID < 0) or the state point was out of
' the valid thermodynamic range (IDID < -100).
      HeCalc = "NUM" 'IErr = 36
   Else
'
' JP defines the desired storage location in Prop (Index, JP)
' JP = -42 if IDID and Phase are incompatible.
'
      If (JP <> -42) Then
         HeCalc = Prop(Index, JP)
      Else
         HeCalc = "N/A" 'IErr = 42
      End If
   End If
End Function


Public Function HeRefrac(Wavelength As Double, Phase As Integer, _
                         Input1 As Variant, Value1 As Double, _
                         Input2 As Variant, Value2 As Double, _
                         Unit As Integer) As Double
' Converted to Visual Basic by Jay Theilacker (1998)
' Wavelength is the input wavelength of light [micrometers], needed for
' refractive index calculation.  In HeProp's Calc, it is Prop(40,0).
   Prop(40, 0) = Wavelength 'First set wavelength in Prop array
   HeRefrac = HeCalc(31, Phase, Input1, Value1, Input2, Value2, Unit)
End Function


Public Function HeUnit(Index As Integer, Unit As Integer) As String
' Converted to Visual Basic by Jay Theilacker (1998)
' Output: HeUnit = units label
' Input:  Index  = Property index in Prop array (Prop(Index,k))
'         Unit   = choice of units (1 - 4)
   If (Unit > 4 Or Unit < 0 Or Index > 44 Or Index < 0) Then
      HeUnit = "[error]"
   Else
      If (Unit <> Kunits(40)) Then Call UniLib(Unit)
      If Index > 39 Then
         HeUnit = Units(Kunits(0))
      Else
         HeUnit = Units(Kunits(Index))
      End If
   End If
End Function


Public Function HeProperty(Index As Integer) As String
' Used to return the property label given by gstrLabels(Index).
   Dim I As Integer
   For I = 1 To 45
      If gintLabelToProp(I) = Index Then Exit For
   Next I
   If I > 44 Then
      HeProperty = ""
   Else
      HeProperty = gstrLabels(I)
   End If
End Function


Public Function HeMsg(ErrorID As Integer) As String
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
' Error messages, when the program returns with ErrorID < zero
'-----Version 3.23
' Input values of ErrorID:
'  -10 TO  -99: external coding errors.  User should investigate.
' -100 to -199: input is outside of the valid EOS range.
' -200 to -299: unexpected iteration convergence failure.  Probably
'               due to input out of range (which should have been
'               caught elsewhere) or, possible programming error.
'               Cryodata should investigate.
   Const Nlines As Integer = 55
   Const NLM1 As Integer = Nlines - 1
   Dim Msg(Nlines) As String * 60
   Dim JID(NLM1) As Integer, J As Integer
    Msg(1) = " Invalid input parameter(s).                                "
    Msg(2) = " Invalid combination of input parameters.                   "
    Msg(3) = " Input Pressure <= zero.                                    "
    Msg(4) = " Input Pressure too high; out of range.                     "
    Msg(5) = " Input Temperature < 0.8 K; out of range.                   "
    Msg(6) = " Input Temperature > 1500 K; out of range.                  "
    Msg(7) = " Input Density <= zero.                                     "
    Msg(8) = " Input Density is outside of valid range.                   "
    Msg(9) = " Solid phase encountered.                                   "
   Msg(10) = " Entropy out of range.                                      "
   Msg(11) = " Enthalpy out of range.                                     "
   Msg(12) = " Internal Energy out of range.                              "
   Msg(13) = " Input Gibbs energy must be in liquid between 0.8 and 3.0 K."
   Msg(14) = " Gibbs energy iteration requires 0.8 <= Temperature <=3.0 K."
   Msg(15) = " Pressure out of range for Gibbs energy iteration.          "
   Msg(16) = " Saturation pressure out of range (1.48 to 226463 Pa).      "
   Msg(17) = " Saturation temperature out of range (0.8 to 5.1953 K).     "
   Msg(18) = " Saturation density out of range (0.000888 to 146.16 Kg/m3)."
   Msg(19) = " Saturation entropy out of range (0.00471 to 23.936 J/g-K). "
   Msg(20) = " Saturation enthalpy out of range (0.00187 to 30.744 J/g).  "
   Msg(21) = " Saturation Helmholz E out of range (-.00191 to -11.287 J/g)"
   Msg(22) = " Saturation energy out of range (0.00186 to 24.724 J/g).    "
   Msg(23) = " Saturation Gibbs E out of range (-0.00189 to -8.0212 J/g). "
   Msg(24) = " Saturation quality must be between 0 and 1 (inclusive).    "
   Msg(25) = " Melting pressure out of range (25.328 to 1013.25 bars)     "
   Msg(26) = " Melting temperature > 14.0 K, pressure exceeds valid range."
   Msg(27) = " Lambda pressure out of range (.050418 to 30.134 bars).     "
   Msg(28) = " Lambda temperature out of range (2.1768 to 1.7673 K).      "
   Msg(29) = " Lambda density out of range (146.15 to 179.83 Kg/m3).      "
   Msg(30) = " On input, absolute value of delta-T must be < 0.5 K        "
   Msg(31) = " Liquid-vapor mixture state point.                          "
   Msg(32) = " Calculated temperature < 0.8 K                             "
   Msg(33) = " Calculated temperature > 1500. K                           "
   Msg(34) = " Calculated pressure > maximum valid pressure.              "
   Msg(35) = " Calculated density <= zero.                                "
   Msg(36) = " Input specific volume <= zero.                             "
   Msg(37) = " Indeterminant function in compressed liquid.               "
   Msg(38) = " Iteration failure with (P,T) input.  Out of range?         "
   Msg(39) = " Iteration failure with (H,S) input.  Out of range?         "
   Msg(40) = " Iteration failure with (P,G) input.  Out of range?         "
   Msg(41) = " Iteration failure with (D,P) input; compressed liquid?     "
   Msg(42) = " Failure with (D,S), (D,H) or (D,U) input.  Out of range?   "
   Msg(43) = " Failure with (P,S), (P,H) or (P,U) input.  Out of range?   "
   Msg(44) = " Unexpected iteration failure near the lambda line.         "
   Msg(45) = " (C) Copyright (1992), Cryodata, Inc.; all rights reserved. "
   Msg(46) = " Successful calculation: single phase fluid.                "
   Msg(47) = " Successful calculation: sat. liquid and vapor.             "
   Msg(48) = " Successful calculation: liquid-vapor mixture.              "
   Msg(49) = " Successful calculation: sat. liquid only.                  "
   Msg(50) = " Successful calculation: sat. vapor only.                   "
   Msg(51) = " Invalid property or units or phase index in HeCalc.        "
   Msg(52) = " HePak has not performed any calculations.                  "
   Msg(53) = " Calculated state point was not in specified fluid phase.   "
   Msg(54) = " Phase and IDID are inconsistent.                           "
   Msg(55) = " Cryodata program error: invalid ErrorID index.             "
'-----Version 3.23
    JID(1) = -11:   JID(2) = -21:   JID(3) = -101:  JID(4) = -102
    JID(5) = -103:  JID(6) = -104:  JID(7) = -105:  JID(8) = -106
    JID(9) = -107: JID(10) = -108: JID(11) = -109: JID(12) = -110
   JID(13) = -111: JID(14) = -112: JID(15) = -113: JID(16) = -121
   JID(17) = -122: JID(18) = -123: JID(19) = -124: JID(20) = -125
   JID(21) = -126: JID(22) = -127: JID(23) = -128: JID(24) = -129
   JID(25) = -131: JID(26) = -132: JID(27) = -141: JID(28) = -142
   JID(29) = -143: JID(30) = -144: JID(31) = -145: JID(32) = -161
   JID(33) = -162: JID(34) = -163: JID(35) = -164: JID(36) = -170
   JID(37) = -180: JID(38) = -201: JID(39) = -202: JID(40) = -203
   JID(41) = -204: JID(42) = -205: JID(43) = -206: JID(44) = -207
   JID(45) = 0:    JID(46) = 1:    JID(47) = 2:    JID(48) = 3
   JID(49) = 4:    JID(50) = 5:    JID(51) = -13:  JID(52) = -20
   JID(53) = -15:  JID(54) = -301
' Note: HEPAK does not assign ErrorID values of -13, -15, or -20.
' However, function HEFUN can return these values.
   For J = 1 To NLM1
      If (JID(J) = ErrorID) Then
         HeMsg = Msg(J)
         Exit Function
      End If
   Next
   HeMsg = Msg(Nlines)
End Function


Public Function HeValidate(P1 As Integer, P2 As Integer) As Integer
' Used to validate an input property pair.
' Differes from JM in that the two inputs are Prop() index numbers.
   Dim In1 As Integer, In2 As Integer
   If (P1 >= 0 And P1 <= 45) And (P2 >= 0 And P2 <= 45) Then
      In1 = gintPropToInput(P1)
      In2 = gintPropToInput(P2)
      HeValidate = JM(In1, In2)
   Else
      HeValidate = 0
   End If
   If HeValidate = 0 Then
      frmPairs.lblPairs.Caption = gstrPairs
      frmPairs.Show
   End If
End Function
 
 
Private Function UserToIn(InASCII As Variant) As Integer
' Converted ITrans to Visual Basic by Jay Theilacker (1998)
' This function maps Prop() indices or ASCII representation to
' the original input indices.
   Dim ASCII As String * 2
   Dim InChar As String * 1
   Dim InList As String * 40
   If TypeName(InASCII) = "String" Then 'String input
      InList = "XPTDV222SH2UG2222222" & _
               "22222222222222222222"
'     Search for two character inputs.
'     Saturation line (satLV, SatL, SatV) or dT to lambda line
      ASCII = Left(CStr(InASCII), 2)
      ASCII = UCase(ASCII)
      Select Case ASCII
         Case "LV"
            UserToIn = -2
         Case "SL"
            UserToIn = 13
         Case "SV"
            UserToIn = 14
         Case "DT"
            UserToIn = 10
         Case Else
            InChar = Left(ASCII, 1)
            Select Case InChar
               Case "M"
                  UserToIn = 11
               Case "L"
                  UserToIn = 15
               Case Else
                  UserToIn = CInt(InStr(InList, InChar))
                  If (UserToIn = 0) Then
                     UserToIn = -2
                  Else
                     UserToIn = UserToIn - 1 'reference to 0 instead of 1
                  End If
            End Select
      End Select
   Else
      UserToIn = gintPropToInput(CInt(InASCII)) 'Integer input
   End If
End Function
 
 
Private Function JPCalc(IPhase As Integer, IDID As Integer) As Integer
' Converted to Visual Basic by Jay Theilacker (1998)

' Calculate JPCalc as a function of IPhase, and IDID
'     IPhase identifies what fluid phase the user wants.
'     IDID identifies all the phases which Calc has calculated.
' JPCalc is the phase index in Prop(i, JP).
'     Calling subroutine must input valid values of IPhase and IDID.
'     0 <= IPhase <= 5  ;   1 <= IDID <= 5
' Errors: JPCalc = -42 when IPhase and IDID are inconsistent. The calling
'     routine will be unable to find a value for the desired parameter.
'
'    JPCalc (IDID, IPhase) is found from the array:
'                   IDID =  1      2      3      4      5
'        IPhase   ----------------------------------------
'          0      |         0      1      0      1      2
'          1      |         0      err    0      err    err
'          2      |         0      err    err    err    err
'          3      |         err    err    0      err    err
'          4      |         err    1      1      1      err
'          5      |         err    2      2      err    2
' Note:
'     IPhase=0 always returns a result, controlled by IDID:
'              single phase when IDID=1, mixture when IDID=3,
'              satL when IDID=2 or 4, satV when IDID=5;
'     IPhase=1 returns either single phase or mixture properties,
'              depending on whether IDID = 1 or IDID=3;
'     IPhase=2 returns single phase fluid only;
'     IPhase=3 returns mixture properties only.
'     IPhase=4 returns saturated liquid properties only;
'     IPhase=5 returns saturated vapor properties only.
'
   Dim M(5, 0 To 5) As Integer
   If (IPhase < 0 Or IPhase > 5 Or IDID < 1 Or IDID > 5) Then
      JPCalc = -42
      Exit Function
   End If
   M(1, 0) = 0
   M(1, 1) = 0
   M(1, 2) = 0
   M(1, 3) = -42
   M(1, 4) = -42
   M(1, 5) = -42
   M(2, 0) = 1
   M(2, 1) = -42
   M(2, 2) = -42
   M(2, 3) = -42
   M(2, 4) = 1
   M(2, 5) = 2
   M(3, 0) = 0
   M(3, 1) = 0
   M(3, 2) = -42
   M(3, 3) = 0
   M(3, 4) = 1
   M(3, 5) = 2
   M(4, 0) = 1
   M(4, 1) = -42
   M(4, 2) = -42
   M(4, 3) = -42
   M(4, 4) = 1
   M(4, 5) = -42
   M(5, 0) = 2
   M(5, 1) = -42
   M(5, 2) = -42
   M(5, 3) = -42
   M(5, 4) = -42
   M(5, 5) = 2
   JPCalc = M(IDID, IPhase)
End Function


Private Function ITrans(InASCII As Variant) As Integer
' Converted to Visual Basic by Jay Theilacker (1998)
   Dim ASCII As String * 2
   Dim InList As String * 15
   If TypeName(InASCII) = "String" Then 'String input
      InList = "PTDVSHUGXtM245L"
'     Search for two character inputs.
'     Saturation line (satLV, SatL, SatV) or dT to lambda line
      ASCII = Left(CStr(InASCII), 2)
      ASCII = UCase(ASCII)
      Select Case ASCII
         Case "LV"
            ITrans = 12
         Case "SL"
            ITrans = 13
         Case "SV"
            ITrans = 14
         Case "DT"
            ITrans = 10
         Case Else
            ITrans = CInt(InStr(InList, Left(ASCII, 1)))
            If (ITrans = 0) Then ITrans = -2
      End Select
   Else
      ITrans = CInt(InASCII) 'Integer input
   End If
End Function


Private Function He(OutParam As Integer, _
            InValue1 As Double, _
            Optional InValue2 As Double) As Double
   Dim Error As Integer
   Dim Dum1 As Double, Dum2 As Double, Dum3 As Double, Dum4 As Double
   Dim D As Double, Q As Double, T As Double, P As Double
   Dim S As Double, H As Double, A As Double, U As Double
   Dim G As Double, Cp As Double, Cv As Double
   Dim DL As Double, DV As Double, HL As Double, HV As Double
   Dim Gamma As Double, Alpha As Double, Grun As Double
   Dim Kt As Double, C As Double, JT As Double
   Dim dPdD As Double, xdPdT As Double, VdHdV As Double
   Select Case OutParam
      Case 0 'Find D as a function of P and T
         Call DfPT(Error, D, Q, InValue1, InValue2)
         He = D
      Case 1 'Find T as a function of P and S
         Call DTfPX(Error, D, T, Q, DL, DV, xdPdT, _
                    InValue1, InValue2, "PS")
         He = T
      Case 2 'Find T as a function of P and H
         Call DTfPX(Error, D, T, Q, DL, DV, xdPdT, _
                    InValue1, InValue2, "PH")
         He = T
      Case 3 'Find T as a function of P and U
         Call DTfPX(Error, D, T, Q, DL, DV, xdPdT, _
                    InValue1, InValue2, "PU")
         He = T
      Case 4 'Find H as a function of P and T
         Call DfPT(Error, D, Q, InValue1, InValue2)
         Call SHAUG(InValue1, InValue2, D, _
                    Dum1, Dum2, Dum3, Dum4, S, H, A, U, G)
         He = H
      Case 5 'Find S as a function of D and T
         He = Entrop(InValue1, InValue2)
      Case 6 'Find Cp as a function of D and T
         Call Deriv(Cp, Cv, Gamma, Alpha, Grun, Kt, C, JT, _
                    dPdD, xdPdT, VdHdV, InValue1, InValue2)
         He = Cp
      Case 7 'Find Cv as a function of D and T
         Call Deriv(Cp, Cv, Gamma, Alpha, Grun, Kt, C, JT, _
                    dPdD, xdPdT, VdHdV, InValue1, InValue2)
         He = Cv
      Case 8 'Find gamma, Cp/Cv, as a function of D and T
         Call Deriv(Cp, Cv, Gamma, Alpha, Grun, Kt, C, JT, _
                    dPdD, xdPdT, VdHdV, InValue1, InValue2)
         He = Gamma
      Case 9 'Find alpha, T/V*dV/dT at constant P,
             'as a function of D and T
         Call Deriv(Cp, Cv, Gamma, Alpha, Grun, Kt, C, JT, _
                    dPdD, xdPdT, VdHdV, InValue1, InValue2)
         He = Alpha
      Case 10 'Find Grun, V/Cv*dP/dT at constant V,
              'as a function of D and T
         Call Deriv(Cp, Cv, Gamma, Alpha, Grun, Kt, C, JT, _
                    dPdD, xdPdT, VdHdV, InValue1, InValue2)
         He = Grun
      Case 11 'Find Kt, 1/D*dDdP at constant T,
              'as a function of D and T
         Call Deriv(Cp, Cv, Gamma, Alpha, Grun, Kt, C, JT, _
                    dPdD, xdPdT, VdHdV, InValue1, InValue2)
         He = Kt
      Case 12 'Find velocity of sound as a function of D and T
         Call Deriv(Cp, Cv, Gamma, Alpha, Grun, Kt, C, JT, _
                    dPdD, xdPdT, VdHdV, InValue1, InValue2)
         He = C
      Case 13 'Find JT, dT/dP at constant H, as a function of D and T
         Call Deriv(Cp, Cv, Gamma, Alpha, Grun, Kt, C, JT, _
                    dPdD, xdPdT, VdHdV, InValue1, InValue2)
         He = JT
      Case 14 'Find dPdD at constant T as a function of D and T
         Call Deriv(Cp, Cv, Gamma, Alpha, Grun, Kt, C, JT, _
                    dPdD, xdPdT, VdHdV, InValue1, InValue2)
         He = dPdD
      Case 15 'Find dPdT at constant D as a function of D and T
         Call Deriv(Cp, Cv, Gamma, Alpha, Grun, Kt, C, JT, _
                    dPdD, xdPdT, VdHdV, InValue1, InValue2)
         He = xdPdT
      Case 16 'Find viscosity as a function of D and T
         He = Viscos(InValue1, InValue2)
      Case 17 'Find thermal conductivity as a function of D and T
         He = TCon(InValue1, InValue2)
      Case 18 'Find Tsat as a function of P
         He = TsatfP(InValue1)
      Case 19 'Find Psat as a function of T
         He = Psat(InValue1)
      Case 20 'Find density of saturated liquid as a function of T
         He = DFsat(InValue1)
      Case 21 'Find density of saturated vapor as a function of T
         He = DGsat(InValue1)
      Case 22 'Find P as a function of D and T
         He = Press(InValue1, InValue2)
      Case 23 'Find superfluid density fraction
              'as a function of P and T
         He = RSRfPT(InValue1, InValue2)
      Case 24 'Find He II viscosity as a function of D and T
         He = VisLam(InValue1, InValue2)
      Case 25 'Find vapor quality as a function of D and P
         Call TfDP(Error, T, Q, DL, DV, InValue1, InValue2)
         He = Q
      Case 26 'Find H as a function of P and T
         Call DfPT(Error, D, Q, InValue1, InValue2)
         Call SHAUG(InValue1, InValue2, D, _
                    Dum1, Dum2, Dum3, Dum4, S, H, A, U, G)
         He = H
      Case 27 'Find S as a function of P and T
         Call DfPT(Error, D, Q, InValue1, InValue2)
         Call SHAUG(InValue1, InValue2, D, _
                    Dum1, Dum2, Dum3, Dum4, S, H, A, U, G)
         He = S
      Case 28 'Find H as a function of P and T
         Call DfPT(Error, D, Q, InValue1, InValue2)
         Call SHAUG(InValue1, InValue2, D, _
                    Dum1, Dum2, Dum3, Dum4, S, H, A, U, G)
         He = H
      Case 29 'Find density of saturated liquid as a function of T
         He = DFsat(InValue1)
      Case 30 'Find vapor quality as a function of D and P
         Call TfDP(Error, T, Q, DL, DV, InValue1, InValue2)
         He = Q
      Case 31 'Find T as a function of D and P
         Call TfDP(Error, T, Q, DL, DV, InValue1, InValue2)
         He = T
      Case 32 'Find vapor quality as a function of P and S
         Call DTfPX(Error, D, T, Q, DL, DV, xdPdT, _
                    InValue1, InValue2, "PS")
         He = Q
      Case 33 'Find vapor quality as a function of P and H
         Call DTfPX(Error, D, T, Q, DL, DV, xdPdT, _
                    InValue1, InValue2, "PH")
         He = Q
      Case 34 'Find vapor quality as a function of P and U
         Call DTfPX(Error, D, T, Q, DL, DV, xdPdT, _
                    InValue1, InValue2, "PU")
         He = Q
      Case 35 'Find U as a function of P and T
         Call DfPT(Error, D, Q, InValue1, InValue2)
         Call SHAUG(InValue1, InValue2, D, _
                    Dum1, Dum2, Dum3, Dum4, S, H, A, U, G)
         He = U
      Case 36 'Find V*dH/dV as a function of D and T
         Call Deriv(Cp, Cv, Gamma, Alpha, Grun, Kt, C, JT, _
                    dPdD, xdPdT, VdHdV, InValue1, InValue2)
         He = VdHdV
      Case 37 'Find Cp as a function of P and T
         Call DfPT(Error, D, Q, InValue1, InValue2)
         Call Deriv(Cp, Cv, Gamma, Alpha, Grun, Kt, C, JT, _
                    dPdD, xdPdT, VdHdV, D, InValue2)
         He = Cp
      Case 38 'Find Cp as a function of D and T
         Call Deriv(Cp, Cv, Gamma, Alpha, Grun, Kt, C, JT, _
                    dPdD, xdPdT, VdHdV, InValue1, InValue2)
         He = Cp
      Case 39 'Find surface tension as a function of T
         He = STen(InValue1) 'N/m
      Case 40 'Find latent heat as a function of T
         He = 0#
         If (T < Tcrit) Then
            P = Psat(InValue1)
            DL = DFsat(InValue1)
            DV = DGsat(InValue1)
            Call SHAUG(P, InValue1, DL, _
                       Dum1, Dum2, Dum3, Dum4, S, HL, A, U, G)
            Call SHAUG(P, InValue1, DV, _
                       Dum1, Dum2, Dum3, Dum4, S, HV, A, U, G)
            He = HV - HL
         End If
      Case Else
         He = -1
   End Select
End Function
      
      
Private Sub Calc(IDID As Integer, PRP() As Double, Jout As Integer, _
                 J1 As Integer, Valu1 As Double, _
                 J2 As Integer, Valu2 As Double, _
                 JPrecs As Integer, JUnits As Integer)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'-----Version 3.30
'
' Please update the Heprop version # as new
'     versions are developed!
'
' This subroutine performs all the calculations.
' !! See file USER.FOR for properties definitions and locations !!
' -------------------------------- AFTER ANY SUCCESSFUL CALL TO CALC
' the following properties will have been calculated:
'    PRP (0:5) + (if IDID > 1) PRP(6:7)
' Additional output is controlled by the 5-digit parameter Jout.
' Each non-zero digit specifies additional parameters.
'   If the 10000's digit > 0:   PRP(8:13)  = State Variables
'   If the 1000's digit > 0 :   PRP(14:24) = Derivatives
'   If the 100's digit > 0  :   PRP(25:26) = Thermal conductivity and
'                               Viscosity and (if IDID > 1) PRP(29) =
'                               surface tension.
'                               Also, if Derivatives have been calculated,
'                               PRP(27:28) = Diffusivity & Prandtl No.
'                               will be calculated.
'   If the 10's digit > 0   :   PRP(30:31) = Dielectric constant and
'                               refractive index.
'       Note! PRP(40,0) = wavelength (micrometers) is a required input for
'       the refractive index calculation!
'   If the 1's digit > 0    :  PRP(8:24) + PRP(32:39) =  HeII properties
'       State and derivative properties will automatically be calculated
'       when HeII properties are requested.
' ------------------------------------------------------
'-----Initialization of constants:
   Const DLamH As Double = 179.83
   Const PLamH As Double = 3013400#
   Const TLamH As Double = 1.7673
   
   Dim Kin(NIV) As Integer
   Dim Ain(NIV) As Double
   Dim M As Integer, K1 As Integer, J As Integer
   Dim JO1 As Integer, JO2 As Integer, JO3 As Integer
   Dim JO4 As Integer, JO5 As Integer
   Dim Dum1 As Double, Dum2 As Double, QAL As Double, OMQ As Double
   Dim CVFJ As Double, WaveLn As Double
   Dim TminM As Double, TmaxA As Double
   Dim LBL1 As String * 1
   Kin(1) = 1: Kin(2) = 2:  Kin(3) = 3:  Kin(4) = 4:  Kin(5) = 8
   Kin(6) = 9: Kin(7) = 11: Kin(8) = 12: Kin(9) = 5: Kin(10) = 33
'
'
' Check for valid input parameters
'
   If (((J1 - 1) * (J1 - NIP) > 0) Or ((J2 - 1) * (J2 - NIP) > 0)) Then
      IDID = -11
      Exit Sub
   End If
'
' Check for valid thermodynamic parameter pair
'
   M = JM(J1, J2)
   If (M <= 0) Then
      IDID = -21
      Exit Sub
   End If
'
' Reset precision and units as needed
'
   If (JPrecs <> CInt(E(0))) Then Call Prcis(JPrecs)
   If (JUnits <> Kunits(40)) Then Call UniLib(JUnits)
'
' Place input values into the array and convert to HEPAK units
' (Note: input temperatures must be K or R, not C or F)
'
   If (J1 <= NIV) Then Ain(J1) = Valu1 * ConvF(Kunits(Kin(J1)))
   If (J2 <= NIV) Then Ain(J2) = Valu2 * ConvF(Kunits(Kin(J2)))
'
'-----Version 3.21B: Add one line for use by spreadsheets.
   If (PRP(41, 2) = -262144#) Then GoTo 200
'
' If Volume was an input parameter, calculate density.
'
   If (J1 = 4) Then
      If (Valu1 <= 0#) Then
         IDID = -170
         Exit Sub
      End If
      Ain(3) = 1# / Valu1
   ElseIf (J2 = 4) Then
      If (Valu2 <= 0#) Then
         IDID = -170
         Exit Sub
      End If
      Ain(3) = 1# / Valu2
   End If
'
' Set the output array to signify "not a valid answer"
'     Properties 40 & 41 are not used for output (in this version),
'     but may be used for various internal communications, or
'     for input.  Do not erase.
   For K1 = 0 To 39
      PRP(K1, 0) = 0#
      PRP(K1, 1) = 0#
      PRP(K1, 2) = 0#
   Next
   IDID = 0
'-----Obtain temperature, density and quality from input parameters
  Select Case M
'-----Pressure, temperature input
   Case 1
      Call DfPT(IDID, PRP(3, 0), PRP(0, 0), Ain(1), Ain(2))
      PRP(1, 0) = Ain(1)
      PRP(2, 0) = Ain(2)
'-----(D,T) Input
   Case 2
      Call FDT(IDID, PRP(0, 0), PRP(1, 0), PRP(6, 1), PRP(3, 1), _
               PRP(3, 2), Ain(3), Ain(2))
      PRP(2, 0) = Ain(2)
      PRP(3, 0) = Ain(3)
'-----Pressure and density input
   Case 3
      Call TfDP(IDID, PRP(2, 0), PRP(0, 0), PRP(3, 1), PRP(3, 2), _
                Ain(3), Ain(1))
      If (IDID = 3) Then Call PsatfT(PRP(1, 0), PRP(6, 1), PRP(2, 0))
      PRP(3, 0) = Ain(3)
'-----Pressure and enthalpy input
   Case 4
      Call DTfPX(IDID, PRP(3, 0), PRP(2, 0), PRP(0, 0), PRP(3, 1), _
                 PRP(3, 2), Dum1, Ain(1), Ain(6), "PH")
      If (IDID = 3) Then Call PsatfT(PRP(1, 0), PRP(6, 1), PRP(2, 0))
'-----Density and enthalpy input
   Case 5
      Call TfDY(IDID, PRP(2, 0), PRP(0, 0), PRP(1, 0), PRP(6, 1), _
                PRP(3, 1), PRP(3, 2), Ain(3), Ain(6), "H")
      PRP(3, 0) = Ain(3)
'-----Pressure and entropy input
   Case 6
      Call DTfPX(IDID, PRP(3, 0), PRP(2, 0), PRP(0, 0), PRP(3, 1), _
                 PRP(3, 2), Dum1, Ain(1), Ain(5), "PS")
      If (IDID = 3) Then Call PsatfT(PRP(1, 0), PRP(6, 1), PRP(2, 0))
'-----Temperature and entropy input
   Case 7
      Call DfST(IDID, PRP(3, 0), PRP(0, 0), PRP(3, 1), PRP(3, 2), _
                PRP(1, 0), PRP(6, 1), Ain(5), Ain(2))
      PRP(2, 0) = Ain(2)
'-----Density and entropy input
   Case 8
      Call TfDY(IDID, PRP(2, 0), PRP(0, 0), PRP(1, 0), PRP(6, 1), _
                PRP(3, 1), PRP(3, 2), Ain(3), Ain(5), "S")
      PRP(3, 0) = Ain(3)
'-----Enthalpy and entropy input
   Case 9
      Call DTfHS(IDID, PRP(3, 0), PRP(2, 0), PRP(0, 0), PRP(3, 1), _
                 PRP(3, 2), PRP(1, 0), PRP(6, 1), Ain(6), Ain(5))
'-----Pressure and internal energy input
   Case 10
      Call DTfPX(IDID, PRP(3, 0), PRP(2, 0), PRP(0, 0), PRP(3, 1), _
                 PRP(3, 2), Dum1, Ain(1), Ain(7), "PU")
      If (IDID = 3) Then Call PsatfT(PRP(1, 0), PRP(6, 1), PRP(2, 0))
'-----Density and internal energy input
   Case 11
      Call TfDY(IDID, PRP(2, 0), PRP(0, 0), PRP(1, 0), PRP(12, 0), _
                PRP(3, 1), PRP(3, 2), Ain(3), Ain(7), "U")
      PRP(3, 0) = Ain(3)
'-----Saturation pressure input, possibly quality and pressure
   Case 12
      Ain(2) = TsatfP(Ain(1))
      If (Ain(2) < 0#) Then
         IDID = -121
         GoTo 200
      End If
      GoTo 1135
'-----delta T to the lambda line along the saturation line
   Case 27
      Ain(2) = Ain(10) + TLamL
      GoTo 2113
'-----Saturation temperature input, possibly quality and temperature
   Case 13
2113  If ((Ain(2) < 0.7999) Or (Ain(2) > Tcrit)) Then
         IDID = -122
         GoTo 200
      End If
1135  PRP(2, 0) = Ain(2)
      Call PsatfT(PRP(1, 0), PRP(6, 1), Ain(2))
      PRP(3, 1) = D2LfPT(PRP(1, 0), Ain(2))
      PRP(3, 2) = D2VfPT(PRP(1, 0), Ain(2))
' If quality has been specified, check that it is between 0 and 1;
' Else with satL or satV specification, set quality = 0
' so that saturated superfluid properties will be calculated correctly.
      If ((J1 = 9) Or (J2 = 9)) Then
         If ((Ain(9) < 0#) Or (Ain(9) > 1#)) Then
            IDID = -129
            GoTo 200
         End If
         IDID = 3
      Else
         Ain(9) = 0#
         IDID = 2
      End If
      PRP(3, 0) = 1# / ((1# - Ain(9)) / PRP(3, 1) + Ain(9) / PRP(3, 2))
      PRP(0, 0) = Ain(9)
'-----Saturation S, H, U, or G input
   Case 14
      If ((J1 = 5) Or (J2 = 5)) Then
         LBL1 = "S"
         Dum1 = Ain(5)
      ElseIf ((J1 = 6) Or (J2 = 6)) Then
         LBL1 = "H"
         Dum1 = Ain(6)
      ElseIf ((J1 = 7) Or (J2 = 7)) Then
         LBL1 = "U"
         Dum1 = Ain(7)
      Else
         LBL1 = "G"
         Dum1 = Ain(8)
      End If
      Call SatfY(IDID, PRP(0, 0), PRP(1, 0), PRP(6, 1), PRP(2, 0), _
                 PRP(3, 1), PRP(3, 2), Dum1, LBL1)
'-----Saturation density input (double-valued region excluded).
   Case 15
      Call SatfD(IDID, PRP(0, 0), PRP(1, 0), PRP(6, 1), PRP(2, 0), _
                 PRP(3, 1), PRP(3, 2), Ain(3))
      PRP(3, 0) = Ain(3)
      If (IDID > 0) Then IDID = 2
'-----Melting pressure input
   Case 16
      PRP(2, 0) = TMfP(Ain(1))
      If (PRP(2, 0) <= 0#) Then
         IDID = -131
      Else
         Call DfPT(IDID, PRP(3, 0), Dum1, Ain(1), PRP(2, 0))
         PRP(0, 0) = -1#
      End If
'-----Melting temperature input
   Case 17
      If ((Ain(2) - 0.8) * (Ain(2) - 14.05) > 0#) Then
         IDID = -132
      Else
         Ain(1) = PMfT(Ain(2))
         Call DfPT(IDID, PRP(3, 0), Dum1, Ain(1), Ain(2))
         PRP(1, 0) = Ain(1)
         PRP(2, 0) = Ain(2)
         PRP(0, 0) = -1#
      End If
'-----Pressure and Gibbs energy input
   Case 18
      Call DTfPG(IDID, PRP(3, 0), PRP(2, 0), Ain(1), Ain(8))
      PRP(0, 0) = -1#
'-----Temperature and Gibbs energy input
   Case 19
      Call DfTG(IDID, PRP(3, 0), Ain(2), Ain(8))
      PRP(2, 0) = Ain(2)
      PRP(0, 0) = -1#
'-----Delta-T to the lambda line at a given density
   Case 20
      If (Abs(Ain(10)) >= 0.5) Then
         IDID = -144
      ElseIf ((Ain(3) > DLamH) Or (Ain(3) < 145.56)) Then
         IDID = -143
      Else
         PRP(1, 0) = PressA(Ain(3), Ain(10))
         PRP(2, 0) = TLam + Ain(10)
         PRP(3, 0) = Ain(3)
' Quality will be adjusted following statement 230
         PRP(0, 0) = -1#
         IDID = 1
' Possible out-or-range near the saturation line
         Call PsatfT(Dum1, Dum2, PRP(2, 0))
         If (PRP(1, 0) < Dum1) Then IDID = -145
      End If
'-----Delta-T to the lambda line at a given pressure
   Case 21
      If (Abs(Ain(10)) >= 0.5) Then
         IDID = -144
      ElseIf ((Ain(1) < 0#) Or (Ain(1) > PLamH)) Then
         IDID = -141
      Else
         IDID = 1
         PRP(2, 0) = TLfP(Ain(1)) + Ain(10)
         Call Df2PT(IDID, PRP(3, 0), Ain(1), PRP(2, 0))
         PRP(1, 0) = Ain(1)
         PRP(0, 0) = -1#
      End If
'---- LAMBDA line pressure input
   Case 22
      If ((Ain(1) < PLamL - 2#) Or (Ain(1) > PLamH * 1.0001)) Then
         IDID = -141
      Else
         PRP(1, 0) = Ain(1)
         PRP(2, 0) = TLfP(Ain(1))
         PRP(3, 0) = DLfT(PRP(2, 0))
         PRP(0, 0) = -2#
         IDID = 1
      End If
'-----LAMBDA line temperature input
   Case 23
      If ((Ain(2) > TLamL + 0.0001) Or (Ain(2) < TLamH - 0.001)) Then
         IDID = -142
      Else
         PRP(1, 0) = PLfT(Ain(2))
         PRP(2, 0) = Ain(2)
         PRP(3, 0) = DLfT(Ain(2))
         PRP(0, 0) = -2#
         IDID = 1
      End If
'-----LAMBDA line density input
   Case 24
      If ((Ain(3) < DLamLf - 0.1) Or (Ain(3) > DLamH + 0.1)) Then
         IDID = -143
      Else
         PRP(2, 0) = TLfD(Ain(3))
         PRP(1, 0) = PLfT(PRP(2, 0))
         PRP(3, 0) = Ain(3)
         PRP(0, 0) = -2#
         IDID = 1
      End If
'---- Lower LAMBDA point
   Case 25
      PRP(1, 0) = PLamL
      PRP(2, 0) = TLamL
      PRP(0, 0) = 0#
      IDID = 2
      PRP(3, 1) = DLamLf
' PRP(3, 2) = D2VFPT(PRP(1, 0), PRP(2, 0))
      PRP(3, 2) = 1.192378
' Following point is defined by the T76 temperature scale
      PRP(12, 0) = 12407.9
'-----Upper LAMBDA point
   Case 26
      PRP(1, 0) = PLamH
      PRP(2, 0) = TLamH
      PRP(3, 0) = DLamH
      PRP(0, 0) = -2#
      IDID = 1
   End Select
'-----All input goes through this point;  check for error.
'     Valid properties at this point are:
'     If IDID=1: density, temperature, and quality
'                = prp(3,0), prp(2,0), prp(0,0)
'                  (In HeII, prp(0,0) will be reset from -1 to -2 later)
'                If M >= 20, pressure = prp(1,0) has been calculated
'     If IDID=3: Mixture, liquid, vapor densities, temperature,
'                saturation pressure, mixture quality
'                = prp(3,0), prp(3,1), prp(3,2), prp(2,0),
' PRP(1, 0), PRP(0, 0)
'     If IDID=2: Liquid and Vapor densities, temperature, pressure,
'                = prp(3,1), prp(3,2), prp(2,0), prp(1,0)
'                  (note: quality = prp(0,0) will be set = 0)
200 'Continue point
   If (IDID >= 1) Then
      If (PRP(2, 0) < 0.7999) Then
         IDID = -161
      ElseIf (PRP(2, 0) > 1500.1) Then
         IDID = -162
      ElseIf (M <= 11) Then
'-----Check for valid density if not already checked
         If (PRP(3, 0) > DenMax(PRP(2, 0))) Then
            IDID = -163
         Else
            Dum2 = Max(PRP(3, 0), PRP(3, 1))
            If (Dum2 <= 0#) Then IDID = -164
         End If
      End If
   End If
   If (IDID = 0) Then IDID = -300
   If (IDID <= 0) Then Exit Sub
'-----Calculate pressure if not already done.
   If (IDID = 1) Then
      If (M < 20) Then PRP(1, 0) = Press(PRP(3, 0), PRP(2, 0))
   Else
'-----In two-phase or mixture, transfer P and T to 2-phase locations
      PRP(1, 1) = PRP(1, 0)
      PRP(1, 2) = PRP(1, 0)
      PRP(2, 1) = PRP(2, 0)
      PRP(2, 2) = PRP(2, 0)
'-----Version 3.23: one line added
      PRP(6, 0) = PRP(6, 1)
      PRP(6, 2) = PRP(6, 1)
      PRP(0, 1) = 0#
      PRP(0, 2) = 1#
   End If
'
'-----P, T, RHO, and X have been evaluated.
'     Also calculate specific volume, PV/RT, and latent heat if at saturation.
   If (IDID <> 2) Then
      PRP(4, 0) = 1# / PRP(3, 0)
      PRP(5, 0) = PRP(1, 0) / (HeGasConst * PRP(2, 0) * PRP(3, 0))
   End If
   If (IDID >= 2) Then
      PRP(5, 1) = PRP(1, 1) / (HeGasConst * PRP(2, 1) * PRP(3, 1))
      PRP(5, 2) = PRP(1, 2) / (HeGasConst * PRP(2, 2) * PRP(3, 2))
      PRP(4, 1) = 1# / PRP(3, 1)
      PRP(4, 2) = 1# / PRP(3, 2)
'-----Version 3.211:
' Clausius-Clapyron value in PRP(7,1) and PRP(7,2);  PRP(7,2) may be
' overwritten later with delH (if delH is calculated).
      PRP(7, 1) = PRP(6, 1) * PRP(2, 0) * (PRP(4, 2) - PRP(4, 1))
      PRP(7, 2) = PRP(7, 1)
'-----Version 3.23: following line
      PRP(7, 0) = PRP(7, 1)
   End If
   QAL = PRP(0, 0)
   OMQ = 1# - QAL
'-----Set output mode
   If (IDID = 2) Then
      If ((J1 = 13) Or (J2 = 13)) Then
         IDID = 4
      ElseIf ((J1 = 14) Or (J2 = 14)) Then
         IDID = 5
      End If
   End If
   JO1 = Jout / 10000
   JO2 = (Jout Mod 10000) / 1000
   JO3 = (Jout Mod 1000) / 100
   JO4 = (Jout Mod 100) / 10
   JO5 = (Jout Mod 10)
   If (JO5 >= 1) Then
      JO1 = 1
      JO2 = 1
   End If
'-----OUTPUT
'-----Version 3.23: non-functional numbered CONTINUEs removed below.
'-----State properties
   If (JO1 > 0) Then
      If (IDID >= 2) Then
         If (IDID <> 5) Then Call SHAUG(PRP(1, 1), PRP(2, 1), PRP(3, 1), _
                                        PRP(4, 1), PRP(5, 1), PRP(6, 1), _
                                        PRP(7, 1), PRP(8, 1), PRP(9, 1), _
                                        PRP(10, 1), PRP(11, 1), PRP(12, 1))
         If (IDID <> 4) Then Call SHAUG(PRP(1, 2), PRP(2, 2), PRP(3, 2), _
                                        PRP(4, 2), PRP(5, 2), PRP(6, 2), _
                                        PRP(7, 2), PRP(8, 2), PRP(9, 2), _
                                        PRP(10, 2), PRP(11, 2), PRP(12, 2))
'-----Version 3.23:
         If (IDID <= 3) Then
            PRP(7, 0) = PRP(9, 2) - PRP(9, 1)
            PRP(7, 2) = PRP(7, 0)
         End If
         PRP(8, 0) = QAL * PRP(8, 2) + OMQ * PRP(8, 1)
         PRP(9, 0) = QAL * PRP(9, 2) + OMQ * PRP(9, 1)
         PRP(10, 0) = QAL * PRP(10, 2) + OMQ * PRP(10, 1)
         PRP(11, 0) = QAL * PRP(11, 2) + OMQ * PRP(11, 1)
         PRP(12, 0) = QAL * PRP(12, 2) + OMQ * PRP(12, 1)
      Else
         Call SHAUG(PRP(1, 0), PRP(2, 0), PRP(3, 0), _
                    PRP(4, 0), PRP(5, 0), PRP(6, 0), _
                    PRP(7, 0), PRP(8, 0), PRP(9, 0), _
                    PRP(10, 0), PRP(11, 0), PRP(12, 0))
      End If
   End If
'-----Derivative properties
   If (JO2 > 0) Then
      If (IDID >= 2) Then
         If (IDID <> 5) Then Call Deriv(PRP(14, 1), PRP(15, 1), PRP(16, 1), _
                                        PRP(17, 1), PRP(18, 1), PRP(19, 1), _
                                        PRP(20, 1), PRP(21, 1), PRP(22, 1), _
                                        PRP(23, 1), PRP(24, 1), _
                                        PRP(3, 1), PRP(2, 1))
         If (IDID <> 4) Then Call Deriv(PRP(14, 2), PRP(15, 2), PRP(16, 2), _
                                        PRP(17, 2), PRP(18, 2), PRP(19, 2), _
                                        PRP(20, 2), PRP(21, 2), PRP(22, 2), _
                                        PRP(23, 2), PRP(24, 2), _
                                        PRP(3, 2), PRP(2, 2))
      Else
         Call Deriv(PRP(14, 0), PRP(15, 0), PRP(16, 0), _
                    PRP(17, 0), PRP(18, 0), PRP(19, 0), _
                    PRP(20, 0), PRP(21, 0), PRP(22, 0), _
                    PRP(23, 0), PRP(24, 0), _
                    PRP(3, 0), PRP(2, 0))
      End If
   End If
'-----Transport properties
   If (JO3 > 0) Then
      If (IDID >= 2) Then
         PRP(29, 1) = STen(PRP(2, 1))
         PRP(29, 2) = PRP(29, 1)
         If (IDID <> 5) Then
            PRP(25, 1) = Viscos(PRP(3, 1), PRP(2, 1))
            PRP(26, 1) = TCon(PRP(3, 1), PRP(2, 1))
            If (PRP(26, 1) > 0#) Then
               PRP(27, 1) = PRP(25, 1) * PRP(14, 1) / PRP(26, 1)
            End If
            If (PRP(14, 1) > 0#) Then
               PRP(28, 1) = PRP(26, 1) / (PRP(3, 1) * PRP(14, 1))
            End If
         End If
         If (IDID <> 4) Then
            PRP(25, 2) = Viscos(PRP(3, 2), PRP(2, 2))
            PRP(26, 2) = TCon(PRP(3, 2), PRP(2, 2))
            If (PRP(26, 2) > 0#) Then
               PRP(27, 2) = PRP(25, 2) * PRP(14, 2) / PRP(26, 2)
            End If
            If (PRP(14, 2) > 0#) Then
               PRP(28, 2) = PRP(26, 2) / (PRP(3, 2) * PRP(14, 2))
            End If
         End If
      ElseIf (PRP(1, 0) < 101000000#) Then
' The transport functions are not fitted at pressures above 100 MPa.
         PRP(25, 0) = Viscos(PRP(3, 0), PRP(2, 0))
         PRP(26, 0) = TCon(PRP(3, 0), PRP(2, 0))
         If (PRP(26, 0) > 0#) Then
            PRP(27, 0) = PRP(25, 0) * PRP(14, 0) / PRP(26, 0)
         End If
         If (PRP(14, 0) > 0#) Then
            PRP(28, 0) = PRP(26, 0) / (PRP(3, 0) * PRP(14, 0))
         End If
      End If
   End If
'-----Dielectric contant and/or refractive index
   If (JO4 > 0) Then
'-----Version 3.23: following IF logic changed.
      WaveLn = PRP(40, 0)
      If ((IDID - 1) * (IDID - 3) = 0) Then
         PRP(30, 0) = DielHe(PRP(3, 0))
         PRP(31, 0) = Refr(PRP(3, 0), WaveLn)
      End If
      If ((IDID - 2) * (IDID - 4) <= 0) Then
         PRP(30, 1) = DielHe(PRP(3, 1))
         PRP(31, 1) = Refr(PRP(3, 1), WaveLn)
      End If
      If ((IDID - 2) * (IDID - 3) * (IDID - 5) = 0) Then
         PRP(30, 2) = DielHe(PRP(3, 2))
         PRP(31, 2) = Refr(PRP(3, 2), WaveLn)
      End If
   End If
'-----Unique superfluid properties
' The array HII(1:5) is evaluated by SuperF
'-----Version 3.23: minor change in Fortran; no functional change.
   If ((JO5 > 0) And (IDID < 5)) Then
      J = Min(1, (IDID - 1))
' Proceed only if A-equations have been used.
      Call AMLap(TminM, TmaxA, PRP(3, J))
      If (PRP(2, J) >= TmaxA) Then GoTo 400
      PRP(39, J) = TLam
      PRP(33, J) = PRP(2, J) - TLam
      If (Abs(PRP(33, J)) < 0.000001) Then PRP(33, J) = 0#
      PRP(32, J) = PRP(2, J) - TLfP(PRP(1, J))
      If (Abs(PRP(32, J)) < 0.000001) Then PRP(32, J) = 0#
      If (PRP(33, J) <= -0.000001) Then
         PRP(0, J) = -2#
'-----Version 3.23: accidently created an index error in following line
'-----Version 3.24: index error corrected: PRP(2,j) corrected to PRP(1,j)
         Call SuperF(PRP(34, J), PRP(35, J), PRP(36, J), PRP(37, J), _
                     PRP(38, J), _
                     PRP(1, J), PRP(3, J), PRP(2, J), _
                     PRP(33, J), PRP(8, J), _
                     PRP(14, J), PRP(15, J), PRP(16, J), PRP(17, J), _
                     PRP(18, J), PRP(19, J), PRP(20, J), PRP(21, J), _
                     PRP(22, J), PRP(23, J), PRP(24, J))
      End If
   End If
400   'Continuation point
'-----Units conversion.
   If (Kunits(40) > 1) Then
      For J = 1 To 38
         CVFJ = ConvF(Kunits(J))
         PRP(J, 0) = PRP(J, 0) / CVFJ
         PRP(J, 1) = PRP(J, 1) / CVFJ
         PRP(J, 2) = PRP(J, 2) / CVFJ
      Next
' Note: PRP(39,k), = Tlambda, remains in Kelvin regardless of units.
   End If
   For J = 0 To 41
      Prop(J, 0) = PRP(J, 0)
      Prop(J, 1) = PRP(J, 1)
      Prop(J, 2) = PRP(J, 2)
   Next J
End Sub


Private Function JM(J1 As Integer, J2 As Integer) As Integer
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Version 3.21, Jan. 5, 1992
'
' Symmetric matrix of response to input variables
   Dim K(15, 15) As Integer
   If (J1 >= 1 And J1 <= 15) And (J2 >= 1 And J2 <= 15) Then
'  PP  TT  DD  VV  SS  HH  UU  GG  XX  dt  MM   2  SL  SV  LL
'   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15
      K(1, 1) = 0:    K(2, 1) = 1:    K(3, 1) = 3    'PP PT PD
      K(4, 1) = 3:    K(5, 1) = 6:    K(6, 1) = 4    'PV PS PH
      K(7, 1) = 10:   K(8, 1) = 18:   K(9, 1) = 12   'PU PG PX
      K(10, 1) = 21:  K(11, 1) = 16:  K(12, 1) = 0   'Pd PM P2
      K(13, 1) = 12:  K(14, 1) = 12:  K(15, 1) = 22  'Pf Pg PL
      K(1, 2) = 1:    K(2, 2) = 0:    K(3, 2) = 2    'TP TT TD
      K(4, 2) = 2:    K(5, 2) = 7:    K(6, 2) = 0    'TV TS TH
      K(7, 2) = 0:    K(8, 2) = 19:   K(9, 2) = 13   'TU TG TX
      K(10, 2) = 0:   K(11, 2) = 17:  K(12, 2) = 0   'Td TM T2
      K(13, 2) = 13:  K(14, 2) = 13:  K(15, 2) = 23  'Tf Tg TL
      K(1, 3) = 3:    K(2, 3) = 2:    K(3, 3) = 0    'DP DT DD
      K(4, 3) = 0:    K(5, 3) = 8:    K(6, 3) = 5    'DV DS DH
      K(7, 3) = 11:   K(8, 3) = 0:    K(9, 3) = 0    'DU DG DX
      K(10, 3) = 20:  K(11, 3) = 0:   K(12, 3) = 0   'Dd DM D2
      K(13, 3) = 15:  K(14, 3) = 15:  K(15, 3) = 24  'Df Dg DL
      K(1, 4) = 3:    K(2, 4) = 2:    K(3, 4) = 0    'VP VT VD
      K(4, 4) = 0:    K(5, 4) = 8:    K(6, 4) = 5    'VV VS VH
      K(7, 4) = 11:   K(8, 4) = 0:    K(9, 4) = 0    'VU VG VX
      K(10, 4) = 20:  K(11, 4) = 0:   K(12, 4) = 0   'Vd VM V2
      K(13, 4) = 15:  K(14, 4) = 15:  K(15, 4) = 24  'Vf Vg VL
      K(1, 5) = 6:    K(2, 5) = 7:    K(3, 5) = 8    'SP ST SD
      K(4, 5) = 8:    K(5, 5) = 0:    K(6, 5) = 9    'SV SS SH
      K(7, 5) = 0:    K(8, 5) = 0:    K(9, 5) = 0    'SU SG SX
      K(10, 5) = 0:   K(11, 5) = 0:   K(12, 5) = 0   'Sd SM S2
      K(13, 5) = 14:  K(14, 5) = 14:  K(15, 5) = 0   'Sf Sg SL
      K(1, 6) = 4:    K(2, 6) = 0:    K(3, 6) = 5    'HP HT HD
      K(4, 6) = 5:    K(5, 6) = 9:    K(6, 6) = 0    'HV HS HH
      K(7, 6) = 0:    K(8, 6) = 0:    K(9, 6) = 0    'HU HG HX
      K(10, 6) = 0:   K(11, 6) = 0:   K(12, 6) = 0   'Hd HM H2
      K(13, 6) = 14:  K(14, 6) = 14:  K(15, 6) = 0   'Hf Hg HL
      K(1, 7) = 10:   K(2, 7) = 0:    K(3, 7) = 11   'UP UT UD
      K(4, 7) = 11:   K(5, 7) = 0:    K(6, 7) = 0    'UV US UH
      K(7, 7) = 0:    K(8, 7) = 0:    K(9, 7) = 0    'UU UG UX
      K(10, 7) = 0:   K(11, 7) = 0:   K(12, 7) = 0   'Ud UM U2
      K(13, 7) = 14:  K(14, 7) = 14:  K(15, 7) = 0   'Uf Ug UL
      K(1, 8) = 18:   K(2, 8) = 19:   K(3, 8) = 0    'GP GT GD
      K(4, 8) = 0:    K(5, 8) = 0:    K(6, 8) = 0    'GV GS GH
      K(7, 8) = 0:    K(8, 8) = 0:    K(9, 8) = 0    'GU GG GX
      K(10, 8) = 0:   K(11, 8) = 0:   K(12, 8) = 0   'Gd GM G2
      K(13, 8) = 14:  K(14, 8) = 14:  K(15, 8) = 0   'Gf Gg GL
      K(1, 9) = 12:   K(2, 9) = 13:   K(3, 9) = 0    'XP XT XD
      K(4, 9) = 0:    K(5, 9) = 0:    K(6, 9) = 0    'XV XS XH
      K(7, 9) = 0:    K(8, 9) = 0:    K(9, 9) = 0    'XU XG XX
      K(10, 9) = 0:   K(11, 9) = 0:   K(12, 9) = 0   'Xd XM X2
      K(13, 9) = 0:   K(14, 9) = 0:   K(15, 9) = 0   'Xf Xg XL
      K(1, 10) = 21:  K(2, 10) = 0:   K(3, 10) = 20  'dP dT dD
      K(4, 10) = 20:  K(5, 10) = 0:   K(6, 10) = 0   'dV dS dH
      K(7, 10) = 0:   K(8, 10) = 0:   K(9, 10) = 0   'dU dG dX
      K(10, 10) = 0:  K(11, 10) = 0:  K(12, 10) = 0  'dd dM d2
      K(13, 10) = 27: K(14, 10) = 0:  K(15, 10) = 0  'df dg dL
      K(1, 11) = 16:  K(2, 11) = 17:  K(3, 11) = 0   'MP MT MD
      K(4, 11) = 0:   K(5, 11) = 0:   K(6, 11) = 0   'MV MS MH
      K(7, 11) = 0:   K(8, 11) = 0:   K(9, 11) = 0   'MU MG MX
      K(10, 11) = 0:  K(11, 11) = 0:  K(12, 11) = 0  'Md MM M2
      K(13, 11) = 0:  K(14, 11) = 0:  K(15, 11) = 26 'Mf Mg ML
      K(1, 12) = 0:   K(2, 12) = 0:   K(3, 12) = 0   '2P 2T 2D
      K(4, 12) = 0:   K(5, 12) = 0:   K(6, 12) = 0   '2V 2S 2H
      K(7, 12) = 0:   K(8, 12) = 0:   K(9, 12) = 0   '2U 2G 2X
      K(10, 12) = 0:  K(11, 12) = 0:  K(12, 12) = 0  '2d 2M 22
      K(13, 12) = 0:  K(14, 12) = 0:  K(15, 12) = 0  '2f 2g 2L
      K(1, 13) = 12:  K(2, 13) = 13:  K(3, 13) = 15  'fP fT fD
      K(4, 13) = 15:  K(5, 13) = 14:  K(6, 13) = 14  'fV fS fH
      K(7, 13) = 14:  K(8, 13) = 14:  K(9, 13) = 0   'fU fG fX
      K(10, 13) = 27: K(11, 13) = 0:  K(12, 13) = 0  'fd fM f2
      K(13, 13) = 0:  K(14, 13) = 0:  K(15, 13) = 25 'ff fg fL
      K(1, 14) = 12:  K(2, 14) = 13:  K(3, 14) = 15  'gP gT gD
      K(4, 14) = 15:  K(5, 14) = 14:  K(6, 14) = 14  'gV gS gH
      K(7, 14) = 14:  K(8, 14) = 14:  K(9, 14) = 0   'gU gG gX
      K(10, 14) = 0:  K(11, 14) = 0:  K(12, 14) = 0  'gd gM g2
      K(13, 14) = 0:  K(14, 14) = 0:  K(15, 14) = 0  'gf gg gL
      K(1, 15) = 22:  K(2, 15) = 23:  K(3, 15) = 24  'LP LT LD
      K(4, 15) = 24:  K(5, 15) = 0:   K(6, 15) = 0   'LV LS LH
      K(7, 15) = 0:   K(8, 15) = 0:   K(9, 15) = 0   'LU LG LX
      K(10, 15) = 0:  K(11, 15) = 26: K(12, 15) = 0  'Ld LM L2
      K(13, 15) = 25: K(14, 15) = 0:  K(15, 15) = 0  'Lf Lg LL
      JM = K(J1, J2)
   Else
      JM = 0
   End If
End Function


Private Function Roun01(X) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Correction of roundoff errors on X in the range 0 to 1
   If (X < 0.00003) Then
      Roun01 = 0#
   ElseIf (X > 0.9999) Then
      Roun01 = 1#
   Else
      Roun01 = X
   End If
End Function


Private Sub Prcis(J As Integer)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Set precision for use in iterative routines.
   Const NF As Integer = 3
   Const NE As Integer = 15
   Dim N As Integer, K As Integer
   Dim A(0 To NE) As Double
   Dim Fact(NF) As Double, R As Double
   A(0) = 0#:          A(1) = 0.03
   A(2) = 0.01:        A(3) = 0.003
   A(4) = 0.001:       A(5) = 0.0003
   A(6) = 0.0001:      A(7) = 0.00003
   A(8) = 0.00001:     A(9) = 0.000003
   A(10) = 0.000001:   A(11) = 0.0000003
   A(12) = 0.0000001:  A(13) = 0.00000003
   A(14) = 0.00000001: A(15) = 0.000000003
   Fact(1) = 20#
   Fact(2) = 0.5
   Fact(3) = 0.002
   If (J > NF) Then
      N = NF
   ElseIf (J < 1) Then
      N = 1
   Else
      N = J
   End If
   R = Fact(N)
   E(0) = N
   For K = 1 To NE
      E(K) = R * A(K)
   Next
End Sub


Private Sub UniLib(Unit As Integer)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
' Output is the global array Kunits
'    the library array U (j,Unit) is copied onto Kunits.
' Kunits is used in subroutine Calc.  Sequence is that of Prop array.
' Kunits may be overwritten in order to customize units.
   Const Unitmax As Integer = 4
   Dim N As Integer, J As Integer
   Dim U(0 To 40, 1 To Unitmax) As Integer
   U(0, 1) = 1:   U(1, 1) = 57:  U(2, 1) = 12:  U(3, 1) = 16:  U(4, 1) = 20
   U(5, 1) = 1:   U(6, 1) = 69:  U(7, 1) = 28:  U(8, 1) = 24:  U(9, 1) = 28
   U(10, 1) = 28: U(11, 1) = 28: U(12, 1) = 28: U(13, 1) = 60: U(14, 1) = 24
   U(15, 1) = 24: U(16, 1) = 1:  U(17, 1) = 65: U(18, 1) = 1:  U(19, 1) = 53
   U(20, 1) = 32: U(21, 1) = 36: U(22, 1) = 73: U(23, 1) = 69: U(24, 1) = 28
   U(25, 1) = 58: U(26, 1) = 38: U(27, 1) = 1:  U(28, 1) = 43: U(29, 1) = 46
   U(30, 1) = 1:  U(31, 1) = 1:  U(32, 1) = 12: U(33, 1) = 12: U(34, 1) = 1
   U(35, 1) = 32: U(36, 1) = 32: U(37, 1) = 49: U(38, 1) = 50: U(39, 1) = 12
   U(40, 1) = 1
   
   U(0, 2) = 1:   U(1, 2) = 2:   U(2, 2) = 12:  U(3, 2) = 16:  U(4, 2) = 21
   U(5, 2) = 1:   U(6, 2) = 70:  U(7, 2) = 29:  U(8, 2) = 25:  U(9, 2) = 29
   U(10, 2) = 29: U(11, 2) = 29: U(12, 2) = 29: U(13, 2) = 4:  U(14, 2) = 25
   U(15, 2) = 25: U(16, 2) = 1:  U(17, 2) = 65: U(18, 2) = 1:  U(19, 2) = 7
   U(20, 2) = 32: U(21, 2) = 34: U(22, 2) = 75: U(23, 2) = 70: U(24, 2) = 29
   U(25, 2) = 40: U(26, 2) = 38: U(27, 2) = 1:  U(28, 2) = 44: U(29, 2) = 47
   U(30, 2) = 1:  U(31, 2) = 1:  U(32, 2) = 12: U(33, 2) = 12: U(34, 2) = 1
   U(35, 2) = 32: U(36, 2) = 32: U(37, 2) = 51: U(38, 2) = 52: U(39, 2) = 12
   U(40, 2) = 2
   
   U(0, 3) = 1:   U(1, 3) = 2:   U(2, 3) = 12:  U(3, 3) = 18:  U(4, 3) = 22
   U(5, 3) = 1:   U(6, 3) = 70:  U(7, 3) = 30:  U(8, 3) = 26:  U(9, 3) = 30
   U(10, 3) = 30: U(11, 3) = 30: U(12, 3) = 30: U(13, 3) = 2:  U(14, 3) = 26
   U(15, 3) = 26: U(16, 3) = 1:  U(17, 3) = 65: U(18, 3) = 1:  U(19, 3) = 7
   U(20, 3) = 32: U(21, 3) = 34: U(22, 3) = 74: U(23, 3) = 70: U(24, 3) = 30
   U(25, 3) = 40: U(26, 3) = 38: U(27, 3) = 1:  U(28, 3) = 43: U(29, 3) = 47
   U(30, 3) = 1:  U(31, 3) = 1:  U(32, 3) = 12: U(33, 3) = 12: U(34, 3) = 1
   U(35, 3) = 32: U(36, 3) = 32: U(37, 3) = 51: U(38, 3) = 52: U(39, 3) = 12
   U(40, 3) = 3
     
   U(0, 4) = 1:   U(1, 4) = 6:   U(2, 4) = 13:  U(3, 4) = 19:  U(4, 4) = 23
   U(5, 4) = 1:   U(6, 4) = 68:  U(7, 4) = 31:  U(8, 4) = 27:  U(9, 4) = 31
   U(10, 4) = 31: U(11, 4) = 31: U(12, 4) = 31: U(13, 4) = 6:  U(14, 4) = 27
   U(15, 4) = 27: U(16, 4) = 1:  U(17, 4) = 65: U(18, 4) = 1:  U(19, 4) = 11
   U(20, 4) = 33: U(21, 4) = 37: U(22, 4) = 76: U(23, 4) = 68: U(24, 4) = 31
   U(25, 4) = 42: U(26, 4) = 39: U(27, 4) = 1:  U(28, 4) = 45: U(29, 4) = 48
   U(30, 4) = 1:  U(31, 4) = 1:  U(32, 4) = 13: U(33, 4) = 13: U(34, 4) = 1
   U(35, 4) = 33: U(36, 4) = 33: U(37, 4) = 49: U(38, 4) = 50: U(39, 4) = 12
   U(40, 4) = 4
' Note: (j,40) is the set number, equivalent to Unit
' If Unit > Unitmax, customized units have been set in an external program:
'     do not overwrite.
   If (Unit > Unitmax) Then Exit Sub
   N = IMax(1, Unit)
   For J = 0 To 40
      Kunits(J) = U(J, N)
   Next
End Sub


Private Sub FDT(IDID As Integer, Q As Double, PS As Double, xdPdT As Double, _
                DL As Double, DV As Double, DD As Double, TT As Double)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
' Check of (density, temperature) input variables for out-of-range
' or liquid-vapor mixtures.
' Converted to Visual Basic by Jay Theilacker (1998)
   Dim D As Double, T As Double
   Dim XDL As Double, XDV As Double
   D = DD
   T = TT
   IDID = 1
   If (T < Tmin) Then
      IDID = -103
   ElseIf (T > Tmax) Then
      IDID = -104
   ElseIf (D <= 0#) Then
      IDID = -105
   ElseIf (D > DenMax(T)) Then
      IDID = -106
   End If
   If (IDID < 0) Then Exit Sub
   If (D > 69.64) Then
      Q = -1#
   Else
      Q = 2#
   End If
   If ((T < Tcrit) And (D < DLamLf)) Then
      XDL = DFsat(T)
      If (D < XDL * 1.05) Then
         XDV = DGsat(T)
         If (D > XDV * 0.95) Then
' Close to or within 2-phase region
            Call PsatfT(PS, xdPdT, T)
            DL = D2LfPT(PS, T)
            If (D < DL) Then
               DV = D2VfPT(PS, T)
               If (D > DV) Then
                  IDID = 3
                  If (DL = DV) Then
                     Q = 0.5
                  Else
                     Q = (DL / D - 1#) / (DL / DV - 1#)
                  End If
               End If
            End If
         End If
      End If
   End If
End Sub


Private Function Press(D As Double, T As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' This function written by V. Arp.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Pressure [Pa] as a function of density [kg/m3] and temperature [K]
' Valid for temperatures 0.8 to 1500 K, but excluding the 2-phase region.
   Dim TL As Double, TH As Double
   Dim PA As Double, WT As Double
   Call AMLap(TL, TH, D)
   If (T > TL) Then
'  He I equation (above 2.5 K)
      Press = PressM((D / GMWT), T) * 1000000#
      If (T >= TH) Then Exit Function
   End If
'  He II and lambda line vicinity
   PA = PressA(D, T)
   If (T <= TL) Then
      Press = PA
   Else
      WT = (T - TL) / (TH - TL)
      Press = WT * Press + (1# - WT) * PA
   End If
End Function


Private Function dPdT(D As Double, T As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' dP/dT [kG/K] at constant V as a function of density [kg/m3]
' and temperature [K]
' Valid for temperatures 0.8 to 1500 K, but excluding the 2-phase region.
   Dim TL As Double, TH As Double
   Dim DA As Double, WT As Double
   Call AMLap(TL, TH, D)
   If (T > TL) Then
      dPdT = dPdTM((D / GMWT), T) * 1000000#
      If (T >= TH) Then Exit Function
   End If
   DA = dPdTA(D, T)
   If (T <= TL) Then
      dPdT = DA
   Else
      WT = (T - TL) / (TH - TL)
      dPdT = WT * dPdT + (1# - WT) * DA
   End If
End Function


Private Function dPdD(D As Double, T As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' dP/dD [Pa-m3/kg] at constant T as a function of density [kg/m3]
' and temperature [K]
' Valid for temperatures 0.8 to 1500 K, but excluding the 2-phase region.
   Dim TL As Double, TH As Double
   Dim DA As Double, WT As Double
   Call AMLap(TL, TH, D)
   If (T > TL) Then
      dPdD = dPdDM((D / GMWT), T) * 249837.6
      If (T >= TH) Then Exit Function
   End If
   If (T >= TH) Then Exit Function
   DA = dPdDA(D, T)
   If (T <= TL) Then
      dPdD = DA
   Else
      WT = (T - TL) / (TH - TL)
      dPdD = WT * dPdD + (1# - WT) * DA
   End If
End Function


Private Function Cv(D As Double, T As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Cv [J/(kg-K)] as a function of density [kg/m3] and temperature [K]
' Valid for temperatures 0.8 to 1500 K, but excluding the 2-phase region.
   Dim TL As Double, TH As Double
   Dim CA As Double, WT As Double
   Call AMLap(TL, TH, D)
   If (T > TL) Then
      Cv = 1000# * CvM((D / GMWT), T) / GMWT
      If (T >= TH) Then Exit Function
   End If
   CA = CvA(D, T)
   If (T <= TL) Then
      Cv = CA
   Else
      WT = (T - TL) / (TH - TL)
      Cv = WT * Cv + (1# - WT) * CA
   End If
End Function


Private Function Entrop(D As Double, T As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Entropy [J/(kg-K)] as a function of density [kg/m3] and temperature [K]
' Valid for temperatures 0.8 to 1500 K, but excluding the 2-phase region.
   Dim TL As Double, TH As Double
   Dim SA As Double, WT As Double
   Call AMLap(TL, TH, D)
   If (T > TL) Then
      Entrop = 1000# * EntrM((D / GMWT), T) / GMWT
      If (T >= TH) Then Exit Function
   End If
   SA = EntrA(D, T)
   If (T <= TL) Then
      Entrop = SA
   Else
      WT = (T - TL) / (TH - TL)
      Entrop = WT * Entrop + (1# - WT) * SA
   End If
End Function


Private Sub SHAUG(P As Double, _
                  T As Double, _
                  D As Double, _
                  Dummy1 As Double, _
                  Dummy2 As Double, _
                  Dummy3 As Double, _
                  Dummy4 As Double, _
                  S As Double, _
                  H As Double, _
                  A As Double, _
                  U As Double, _
                  G As Double)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Output: S, H, A, U, G;  Input: P, D, T ; all SI units
'    PRP: 1=P, 2=T, 3=D, 8=S, 9=H, 10=A, 11=U, 12=G, 13=Fugacity
' Valid for temperatures 0.8 to 1500 K, but excluding the 2-phase region.
   Dim PV As Double, TL As Double, TH As Double
   Dim WT As Double, OMWT As Double
   Dim SA As Double, AA As Double, UA As Double
   Dim HA As Double, GA As Double
   Call AMLap(TL, TH, D)
   PV = P / D
   If (T > TL) Then
      S = 1000# * EntrM((D / GMWT), T) / GMWT
      U = 1000# * EnerM((D / GMWT), T) / GMWT
      H = U + PV
      A = U - T * S
      G = A + PV
      If (T >= TH) Then Exit Sub
   End If
   SA = EntrA(D, T)
   AA = HelmA(D, T)
   UA = AA + T * SA
   HA = UA + PV
   GA = AA + PV
   If (T <= TL) Then
      S = SA
      H = HA
      A = AA
      U = UA
      G = GA
   Else
      WT = (T - TL) / (TH - TL)
      OMWT = 1# - WT
      S = WT * S + OMWT * SA
      H = WT * H + OMWT * HA
      A = WT * A + OMWT * AA
      U = WT * U + OMWT * UA
      G = WT * G + OMWT * GA
   End If
End Sub


Private Function CvM(Rho As Double, T As Double) As Double
' (C) Copyright (1994), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Output: CvM [J/mol-K] = Cv = specific heat at constant volume
' Input:  RHO [mol/L]   = density (2-phase region excluded)
'         T [K]         = temperature, range 2 to 1500 K
'
'-----Version 3.3; function rewritten for greater efficiency
   Dim TT As Double, TI As Double, TI2 As Double
   Dim D As Double, D2 As Double, D4 As Double
   Dim F As Double, GMI As Double
   Dim GG2 As Double, GG3 As Double, GG4 As Double, GG5 As Double
   Dim GG6 As Double, GD1 As Double, GD2 As Double, GD3 As Double
   Dim GD4 As Double, GD5 As Double, GD6 As Double, CIN As Double
   Dim P1 As Double, P2 As Double, P3 As Double, P4 As Double
   Dim P5 As Double, P6 As Double, P7 As Double, P8 As Double
   Const CPIG0 As Double = 2.5 * 0.0083143
   Const Six7TH As Double = 6# / 7#
   Const F1OP3 As Double = 1# / 0.3
   TT = T
   D = Rho
   D2 = D * D
   D4 = D2 * D2
   TI = 1# / TT
   TI2 = TI * TI
   F = Exp(HeIEoSGamma * D2)
   GMI = 0.5 / HeIEoSGamma
   GG2 = -2# * GMI * GMI
   GG3 = -4# * GG2 * GMI
   GG4 = -6# * GG3 * GMI
   GG5 = -8# * GG4 * GMI
   GG6 = -10# * GG5 * GMI
   GD1 = F * GMI
   GD2 = (F * D2 - 2# * GD1) * GMI
   GD3 = (F * D4 - 4# * GD2) * GMI
   GD4 = (F * D4 * D2 - 6# * GD3) * GMI
   GD5 = (F * D4 * D4 - 8# * GD4) * GMI
   GD6 = (F * D4 * D4 * D2 - 10# * GD5) * GMI
   P1 = (GD6 - GG6) * (HeIEoS(30) _
                       + TI * (2# * HeIEoS(31) _
                               + F1OP3 * TI * HeIEoS(32)))
   P2 = (GD5 - GG5) * (HeIEoS(28) + 2# * TI * HeIEoS(29))
   P3 = (GD4 - GG4) * (HeIEoS(26) + F1OP3 * TI2 * HeIEoS(27))
   P4 = (GD3 - GG3) * (HeIEoS(24) + 2# * TI * HeIEoS(25))
   P5 = (GD2 - GG2) * (HeIEoS(22) + F1OP3 * TI2 * HeIEoS(23))
   P6 = (GD1 - GMI) * (HeIEoS(20) + 2# * TI * HeIEoS(21))
   P7 = (0.75 * HeIEoS(19) * D + Six7TH * HeIEoS(18)) * TI _
        + Two7TH * HeIEoS(17)
   P8 = (P7 * D + One3RD * HeIEoS(16)) * D + 1.2 * TI * HeIEoS(15) _
        + 0.4 * HeIEoS(14)
   CIN = D * (-0.25 * HeIEoS(2) / Sqr(TT) _
              + TI2 * (2# * HeIEoS(4) + 6# * HeIEoS(5) * TI _
                       + D * (HeIEoS(8) + 3# * HeIEoS(9) * TI _
                              + Two3RD * HeIEoS(12) * D)) _
              + D4 * TI2 * P8) _
         + 6# * TI2 * TI * (P6 + P5 + P4 + P3 + P2 + P1)
   CvM = 1000# * (CPIG0 - (UnivGasConst + CIN))
End Function


Private Function dPdDM(RR As Double, TT As Double) As Double
' (C) Copyright (1994), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Output: dPdDM [MPa-L/mol] = dPressure/dDensity at constant temperature
' Input:  RR [mol/L]        = density (2-phase region excluded)
'         TT [K]            = temperature, range 2 to 1500 K
'
'-----Version 3.3; function rewritten for greater efficiency
   Const C53 As Double = 5# / 3#
   Const C73 As Double = 7# / 3#
   Const C113 As Double = 11# / 3#
   Const C133 As Double = 13# / 3#
   Dim D As Double, D2 As Double
   Dim T As Double, TI As Double, TI2 As Double
   Dim F As Double, F1 As Double, FZ As Double, FY As Double
   Dim P1 As Double, P2 As Double, P3 As Double, P4 As Double
   Dim P5 As Double, P6 As Double, P7 As Double, P8 As Double
   Dim P9 As Double, P10 As Double, P11 As Double
   D = RR
   T = TT
   D2 = D * D
   TI = 1# / TT
   TI2 = TI * TI
   F = Exp(HeIEoSGamma * D2)
   F1 = 2# * F * HeIEoSGamma * D
   FZ = F1 * D2 * D
   FY = 3# * F * D2
   P1 = (C133 * FY + FZ) * (HeIEoS(30) + TI * HeIEoS(31) _
                            + TI2 * HeIEoS(32))
   P2 = (C113 * FY + FZ) * (HeIEoS(28) + TI * HeIEoS(29)) + D2 * P1
   P3 = (3# * FY + FZ) * (HeIEoS(26) + TI2 * HeIEoS(27)) + D2 * P2
   P4 = (C73 * FY + FZ) * (HeIEoS(24) + TI * HeIEoS(25)) + D2 * P3
   P5 = (C53 * FY + FZ) * (HeIEoS(22) + TI2 * HeIEoS(23)) + D2 * P4
   P6 = 8# * (HeIEoS(17) + TI * HeIEoS(18)) _
        + 9# * TI * D * HeIEoS(19)
   P7 = 6# * (HeIEoS(14) + TI * HeIEoS(15)) _
        + D * 7# * HeIEoS(16) + D2 * P6
   P8 = 4# * (T * HeIEoS(10) + HeIEoS(11) + TI * HeIEoS(12) _
              + D * 1.25 * HeIEoS(13)) + D2 * TI * P7
   P9 = T * HeIEoS(6) + HeIEoS(7) + TI * (HeIEoS(8) + TI * HeIEoS(9))
   P10 = HeIEoS(3) + HeIEoS(2) * Sqr(T) + HeIEoS(1) * T + _
           TI * (HeIEoS(4) + TI * HeIEoS(5))
   P11 = (FY + FZ) * (HeIEoS(20) + TI * HeIEoS(21)) + D2 * P5
   dPdDM = D * (2# * P10 + D * (3# * P9 + D * P8)) _
           + TI2 * P11 + UnivGasConst * T
End Function


Private Function dPdTM(RR As Double, TT As Double) As Double
' (C) Copyright (1994), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Output: dPdTM [MPa/K] = dPressure/dTemperature at constant density
' Input:  RR [mol/L]    = density (2-phase region excluded)
'         TT [K]        = temperature, range 2 to 1500 K
'
'-----Version 3.3; function rewritten for greater efficiency
   Dim D As Double, D2 As Double
   Dim T As Double, TI As Double, TI2 As Double
   Dim P1 As Double, P2 As Double, P3 As Double, P4 As Double
   Dim P5 As Double, P6 As Double, P7 As Double, P8 As Double
   Dim P9 As Double, P10 As Double
   D = RR
   T = TT
   D2 = D * D
   TI = 1# / T
   TI2 = TI * TI
   P1 = 2# * HeIEoS(30) + 3# * TI * HeIEoS(31) + 4# * TI2 * HeIEoS(32)
   P2 = 2# * HeIEoS(28) + 3# * TI * HeIEoS(29) + D2 * P1
   P3 = 2# * HeIEoS(26) + 4# * TI2 * HeIEoS(27) + D2 * P2
   P4 = 2# * HeIEoS(24) + 3# * TI * HeIEoS(25) + D2 * P3
   P5 = 2# * HeIEoS(22) + 4# * TI2 * HeIEoS(23) + D2 * P4
   P6 = 2# * HeIEoS(20) + 3# * TI * HeIEoS(21) + D2 * P5
   P7 = 2# * HeIEoS(9) + Exp(HeIEoSGamma * D2) * P6
   P8 = HeIEoS(17) + 2# * TI * (HeIEoS(18) + D * HeIEoS(19))
   P9 = HeIEoS(14) + 2# * TI * HeIEoS(15) + D * HeIEoS(16) + D2 * P8
   P10 = HeIEoS(10) - TI2 * HeIEoS(12) - D2 * TI2 * P9
   dPdTM = D2 * (HeIEoS(1) + 0.5 * HeIEoS(2) / Sqr(T) _
                 - TI2 * (HeIEoS(4) + 2# * TI * HeIEoS(5)) _
                 + D * (HeIEoS(6) - TI2 * (HeIEoS(8) + TI * P7) _
                        + D * P10)) _
           + UnivGasConst * D
End Function


Private Function EnerM(RR As Double, TT As Double) As Double
' (C) Copyright (1994), Cryodata Inc.; all rights reserved.
'
'  Output: EnerM [J/mol] = Internal energy
'  Input:  RR [mol/L]    = density (2-phase region excluded)
'          TT [K]        = temperature, range 2 to 1500 K
'
'-----Version 3.3; function rewritten for greater efficiency
   Dim F As Double, GMI As Double, H0 As Double
   Dim T As Double, TI As Double, TI2 As Double
   Dim D As Double, D2 As Double, D3 As Double, D4 As Double
   Dim GG2 As Double, GG3 As Double, GG4 As Double, GG5 As Double
   Dim GG6 As Double, GD1 As Double, GD2 As Double, GD3 As Double
   Dim GD4 As Double, GD5 As Double, GD6 As Double, En As Double
   Dim E1 As Double, E2 As Double, E3 As Double
   Dim E4 As Double, E5 As Double, E6 As Double
   Dim E7 As Double, E8 As Double, E9 As Double
   Dim E10 As Double, E11 As Double
   Const Three7 As Double = 3# / 7#
   H0 = 0.0611315
   T = TT
   TI = 1# / TT
   TI2 = TI * TI
   D = RR
   D2 = D * D
   D3 = D2 * D
   D4 = D3 * D
   F = Exp(HeIEoSGamma * D2)
   GMI = 0.5 / HeIEoSGamma
   GG2 = -2# * GMI * GMI
   GG3 = -4# * GG2 * GMI
   GG4 = -6# * GG3 * GMI
   GG5 = -8# * GG4 * GMI
   GG6 = -10# * GG5 * GMI
   GD1 = F * GMI
   GD2 = (F * D2 - 2# * GD1) * GMI
   GD3 = (F * D4 - 4# * GD2) * GMI
   GD4 = (F * D4 * D2 - 6# * GD3) * GMI
   GD5 = (F * D4 * D4 - 8# * GD4) * GMI
   GD6 = (F * D4 * D4 * D2 - 10# * GD5) * GMI
   E1 = HeIEoS(17) * Two7TH + TI * (HeIEoS(18) * Three7 _
        + D * HeIEoS(19) * 0.375)
   E2 = HeIEoS(14) * 0.4 + TI * HeIEoS(15) * 0.6 _
        + D * (HeIEoS(16) * One3RD + D * E1)
   E3 = HeIEoS(11) * One3RD + TI * HeIEoS(12) * Two3RD _
        + D * (HeIEoS(13) * 0.25 + D * TI * E2)
   E4 = (GD1 - GMI) * (3# * HeIEoS(20) + TI * 4# * HeIEoS(21))
   E5 = (GD2 - GG2) * (3# * HeIEoS(22) + TI2 * 5# * HeIEoS(23))
   E6 = (GD3 - GG3) * (3# * HeIEoS(24) + TI * 4# * HeIEoS(25))
   E7 = (GD4 - GG4) * (3# * HeIEoS(26) + TI2 * 5# * HeIEoS(27))
   E8 = (GD5 - GG5) * (3# * HeIEoS(28) + TI * 4# * HeIEoS(29))
   E9 = (GD6 - GG6) * (3# * HeIEoS(30) + TI * (4# * HeIEoS(31) _
        + TI * 5# * HeIEoS(32)))
   E10 = HeIEoS(7) * 0.5 _
         + TI * (HeIEoS(8) + TI * HeIEoS(9) * 1.5) + D * E3
   E11 = HeIEoS(2) * 0.5 * Sqr(T) + HeIEoS(3) _
         + TI * (HeIEoS(4) * 2# + TI * HeIEoS(5) * 3#) + D * E10
   En = 2.5 * UnivGasConst * T + D * E11 _
        + TI2 * (E4 + E5 + E6 + E7 + E8 + E9)
   EnerM = (En - UnivGasConst * T + H0) * 1000#
End Function


Private Function EntrM(RR As Double, TT As Double) As Double
' (C) Copyright (1994), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
'
' Output: EntrM [J/mol-K] = entropy
' Input:  RR [mol/L]      = density (2-phase region excluded)
'         TT [K]          = temperature, range 2 to 1500 K
'
'-----Version 3.3; function rewritten for greater efficiency
   Const One6TH As Double = 1# / 6#
   Const One7TH As Double = 1# / 7#
   Dim F As Double, S0 As Double
   Dim T As Double, TI As Double, TI2 As Double
   Dim D As Double, D2 As Double, D4 As Double
   Dim GG1 As Double, GG2 As Double, GG3 As Double, GG4 As Double
   Dim GG5 As Double, GG6 As Double, GD1 As Double, GD2 As Double
   Dim GD3 As Double, GD4 As Double, GD5 As Double, GD6 As Double
   Dim S1 As Double, S2 As Double, S3 As Double, S4 As Double
   Dim S5 As Double, S6 As Double, S7 As Double, S8 As Double
   Dim S9 As Double, S10 As Double, S11 As Double
   Dim Sint As Double
   S0 = 0.0078616
   T = TT
   TI = 1# / T
   TI2 = TI * TI
   D = RR
   D2 = D * D
   D4 = D2 * D2
   F = Exp(HeIEoSGamma * D2)
   GG1 = 0.5 / HeIEoSGamma
   GG2 = -2# * GG1 * GG1
   GG3 = -4# * GG2 * GG1
   GG4 = -6# * GG3 * GG1
   GG5 = -8# * GG4 * GG1
   GG6 = -10# * GG5 * GG1
   GD1 = F * GG1
   GD2 = (F * D2 - 2# * GD1) * GG1
   GD3 = (F * D4 - 4# * GD2) * GG1
   GD4 = (F * D4 * D2 - 6# * GD3) * GG1
   GD5 = (F * D4 * D4 - 8# * GD4) * GG1
   GD6 = (F * D4 * D4 * D2 - 10# * GD5) * GG1
'vda Z(4) = 2.5
   S1 = 0.4 * HeIEoS(15) _
        + D2 * (Two7TH * HeIEoS(18) + D * 0.25 * HeIEoS(19))
   S2 = 2# * HeIEoS(5) + D * (HeIEoS(9) + D2 * D * S1)
   S3 = 0.2 * HeIEoS(14) + D * (One6TH * HeIEoS(16) _
          + D * One7TH * HeIEoS(17))
   S4 = 0.5 * HeIEoS(8) + D * (One3RD * HeIEoS(12) + D2 * S3)
   S5 = -HeIEoS(1) - 0.5 * HeIEoS(2) / Sqr(T) _
        - D * (0.5 * HeIEoS(6) + One3RD * D * HeIEoS(10)) _
        + TI2 * (HeIEoS(4) + D * S4 + TI * S2)
   S6 = (GD6 - GG6) * (HeIEoS(30) + HeIEoS(31) * TI * 1.5 _
          + HeIEoS(32) * TI2 * 2#)
   S7 = (GD5 - GG5) * (HeIEoS(28) + HeIEoS(29) * TI * 1.5)
   S8 = (GD4 - GG4) * (HeIEoS(26) + HeIEoS(27) * TI2 * 2#)
   S9 = (GD3 - GG3) * (HeIEoS(24) + HeIEoS(25) * TI * 1.5)
   S10 = (GD2 - GG2) * (HeIEoS(22) + HeIEoS(23) * TI2 * 2#)
   S11 = (F - 1#) * GG1 * (HeIEoS(20) + HeIEoS(21) * TI * 1.5)
   Sint = 2.5 * UnivGasConst * Log(T) + D * S5 _
          + 2# * TI2 * TI * (S11 + S10 + S9 + S8 + S7 + S6)
   EntrM = (Sint + S0 _
           - UnivGasConst * Log(D * UnivGasConst * TT / Pref)) * 1000#
End Function


Private Function PressM(RR As Double, TT As Double) As Double
' (C) Copyright (1994), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Output: PressM [MPa] = pressure
' Input:  RR [mol/L]   = density (2-phase region excluded)
'         TT [K]       = temperature, range 2 to 1500 K
'
'-----Version 3.3; function rewritten for greater efficiency
   Dim D As Double, D2 As Double
   Dim T As Double, TI As Double, TI2 As Double
   Dim P1 As Double, P2 As Double, P3 As Double
   Dim P4 As Double, P5 As Double, P6 As Double
   T = TT
   D = RR
   D2 = D * D
   TI = 1# / TT
   TI2 = TI * TI
   P1 = ((TI * HeIEoS(32) + HeIEoS(31)) * TI + HeIEoS(30)) * D2 _
        + HeIEoS(29) * TI
   P2 = ((P1 + HeIEoS(28)) * D2 + HeIEoS(27) * TI2 + HeIEoS(26)) * D2 _
        + HeIEoS(25) * TI
   P3 = ((P2 + HeIEoS(24)) * D2 + HeIEoS(23) * TI2 + HeIEoS(22)) * D2 _
        + HeIEoS(21) * TI
   P4 = ((P3 + HeIEoS(20)) * Exp(HeIEoSGamma * D2) + _
            HeIEoS(9)) * TI + HeIEoS(8)
   P5 = ((D * HeIEoS(19) + HeIEoS(18)) * TI + HeIEoS(17)) * D _
        + HeIEoS(16)
   P6 = ((P5 * D + TI * HeIEoS(15) + HeIEoS(14)) * D2 + HeIEoS(12)) * TI _
        + D * HeIEoS(13) + HeIEoS(11)
   PressM = ((P4 * TI + HeIEoS(7)) * D + P6 * D2 + _
            (TI * HeIEoS(5) + HeIEoS(4)) * TI + Sqr(T) * HeIEoS(2) + _
            ((D * HeIEoS(10) + HeIEoS(6)) * D + HeIEoS(1)) * T + _
            HeIEoS(3)) * D2 + UnivGasConst * D * T
End Function


Private Function PressA(D As Double, TT As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Pressure [Pa] as a function of density [kg/m3] and temperature [K]
' Valid in compressed liquid from 0.8 to about 3 K.
' If TT is less than 0.8 K, TT is assumed to be the isochoric distance
'    to the lambda line (negative in HeII, positive in HeI).
   Dim PVTCSA(6) As Double
   Dim V As Double, X As Double, T As Double, DT As Double
   Dim M As Integer
   V = 1000# / D
   X = V - VLamLf
   If (V <> Vsave) Then Call LamDer(V)
   If (TT > 0.7999) Then
      T = TT
      DT = T - TLam
   Else
      DT = TT
      T = TLam + DT
   End If
   If (DT > 0#) Then
      M = 1
   Else
      M = 2
   End If
   Call LogFun(PVTCSA, T, DT, M)
   PressA = (PVTCSA(1) + PresA2(X, T, M)) * 1000000#
End Function


Private Function PresA2(X As Double, T As Double, M As Integer) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' "background" pressure [MPa] as a function of V-VLamLf [cm3/gm], and T [K].
' M specifies HeI or HeII.  VLamLf is the volume at the lower lambda point.
' Called from PressA.
   Dim F(7) As Double, Q(5) As Double, A1 As Double
   Dim T2 As Double, J As Integer, K As Integer
   T2 = T * T
   F(1) = 1#
   For J = 1 To 5
      Q(J) = 0#
   Next
   For K = 1 To 6
      For J = 1 To 5
         Q(J) = Q(J) + F(K) * HeIIEoS(6 * J + K - 1, M)
      Next
      F(K + 1) = F(K) * X
   Next
   A1 = Q(3) + T2 * (Q(4) + T2 * Q(5))
   PresA2 = Q(1) + T2 * (Q(2) + T2 * A1)
End Function


Private Function dPdTA(D As Double, TT As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' dP/dT as a function of density and temperature [SI units]
' Valid in compressed liquid from 0.8 to about 3 K.
' If TT is less than 0.8 K, TT is assumed to be the isochoric distance
'    to the lambda line (negative in HeII, positive in HeI).
   Dim F(7) As Double, Q(4) As Double, PVTCSA(6) As Double
   Dim T As Double, T2 As Double, DT As Double
   Dim V As Double, X As Double, A1 As Double
   Dim M As Integer, J As Integer, K As Integer
   V = 1000# / D
   X = V - VLamLf
   If (V <> Vsave) Then Call LamDer(V)
   If (TT > 0.7999) Then
      T = TT
      DT = T - TLam
   Else
      DT = TT
      T = TLam + DT
   End If
   If (DT > 0#) Then
      M = 1
   Else
      M = 2
   End If
   Call LogFun(PVTCSA, T, DT, M)
   T2 = T * T
   F(1) = 1#
   For J = 1 To 4
      Q(J) = 0#
   Next
   For K = 1 To 6
      For J = 1 To 4
         Q(J) = Q(J) + F(K) * HeIIEoS(6 * J + K + 5, M)
      Next
      F(K + 1) = F(K) * X
   Next
   A1 = 4# * Q(2) + T2 * (6# * Q(3) + T2 * 8# * Q(4))
   dPdTA = (PVTCSA(3) + T * (2# * Q(1) + T2 * A1)) * 1000000#
End Function


Private Function dPdDA(D As Double, TT As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' dP/dD as a function of density and temperature [SI units]
' Valid in compressed liquid from 0.8 to about 3 K.
' If TT is less than 0.8 K, TT is assumed to be the isochoric distance
'    to the lambda line (negative in HeII, positive in HeI).
   Dim F(6) As Double, Q(5) As Double, PVTCSA(6) As Double
   Dim T As Double, T2 As Double, DT As Double, dPdV As Double
   Dim V As Double, X As Double, FF As Double, A1 As Double
   Dim M As Integer, J As Integer, K As Integer
   V = 1000# / D
   X = V - VLamLf
   If (V <> Vsave) Then Call LamDer(V)
   If (TT > 0.7999) Then
      T = TT
      DT = T - TLam
   Else
      DT = TT
      T = TLam + DT
   End If
    If (DT > 0#) Then
      M = 1
   Else
      M = 2
   End If
   Call LogFun(PVTCSA, T, DT, M)
   T2 = T * T
   F(1) = 1#
   For J = 1 To 5
      Q(J) = 0#
   Next
   For K = 1 To 5
      FF = F(K) * CDbl(K)
      For J = 1 To 5
         Q(J) = Q(J) + FF * HeIIEoS(6 * J + K, M)
      Next
      F(K + 1) = F(K) * X
   Next
   A1 = Q(3) + T2 * (Q(4) + T2 * Q(5))
   dPdV = PVTCSA(2) + Q(1) + T2 * (Q(2) + T2 * A1)
   dPdDA = -1000000000# * dPdV / (D * D)
End Function


Private Function CvA(D As Double, TT As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Cv as a function of density and temperature [SI units]
' Valid in compressed liquid from 0.8 to about 3 K.
' If TT is less than 0.8 K, TT is assumed to be the isochoric distance
'    to the lambda line (negative in HeII, positive in HeI).
   Dim F(7) As Double, Q(4) As Double, PVTCSA(6) As Double
   Dim T As Double, T2 As Double, DT As Double, Y1 As Double
   Dim V As Double, X As Double, FF As Double, A1 As Double
   Dim M As Integer, J As Integer, K As Integer
   V = 1000# / D
   X = V - VLamLf
   If (V <> Vsave) Then Call LamDer(V)
   If (TT > 0.7999) Then
      T = TT
      DT = T - TLam
   Else
      DT = TT
      T = TLam + DT
   End If
   If (DT > 0#) Then
      M = 1
   Else
      M = 2
   End If
   Call LogFun(PVTCSA, T, DT, M)
   T2 = T * T
   F(1) = X
   For J = 1 To 4
      Q(J) = 0#
   Next
   For K = 1 To 6
      FF = F(K) / CDbl(K)
      For J = 1 To 4
         Q(J) = Q(J) + FF * HeIIEoS(6 * J + K + 5, M)
      Next
      F(K + 1) = F(K) * X
   Next
   A1 = HeIIEoS(37, M) + T2 * (HeIIEoS(38, M) + T2 * HeIIEoS(39, M))
   Y1 = 12# * Q(2) + T2 * (30# * Q(3) + T2 * 56# * Q(4))
   CvA = (PVTCSA(4) + T * (2# * Q(1) + T2 * Y1) + _
         T * (HeIIEoS(36, M) + T2 * A1)) * 1000#
End Function


Private Function EntrA(D As Double, TT As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Entropy as a function of density and temperature [SI units]
' Valid in compressed liquid from 0.8 to about 3 K.
' If TT is less than 0.8 K, TT is assumed to be the isochoric distance
'    to the lambda line (negative in HeII, positive in HeI).
   Dim PVTCSA(6) As Double
   Dim T As Double, DT As Double, M As Integer
   Dim V As Double, X As Double
   V = 1000# / D
   X = V - VLamLf
   If (V <> Vsave) Then Call LamDer(V)
   If (TT > 0.7999) Then
      T = TT
      DT = T - TLam
   Else
      DT = TT
      T = TLam + DT
   End If
   If (DT > 0#) Then
      M = 1
   Else
      M = 2
   End If
   Call LogFun(PVTCSA, T, DT, M)
   EntrA = (PVTCSA(5) + EntrA2(X, T, M)) * 1000#
End Function


Private Function EntrA2(X As Double, T As Double, M As Integer) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' "background" entropy [J/gm-K] as a function of V-V0 [cm3/gm], and T [K].
' M specifies HeI or HeII.  V0 is the volume at the lower lambda point.
' Called from EntrA.
   Dim F(7) As Double, Q(4) As Double, A1 As Double
   Dim T2 As Double, T5 As Double, FF As Double
   Dim J As Integer, K As Integer
   T2 = T * T
   T5 = T2 * T2 * T
   F(1) = X
   For J = 1 To 4
      Q(J) = 0#
   Next
   For K = 1 To 6
      FF = F(K) / CDbl(K)
      For J = 1 To 4
         Q(J) = Q(J) + FF * HeIIEoS(6 * J + K + 5, M)
      Next
      F(K + 1) = F(K) * X
   Next
   A1 = 4# * Q(2) + T2 * (6# * Q(3) + T2 * 8# * Q(4))
   EntrA2 = T * (2# * Q(1) + T2 * A1) + HeIIEoS(36, M) * T + _
            HeIIEoS(37, M) * T2 * T / 3# + HeIIEoS(38, M) * T5 / 5# + _
            HeIIEoS(39, M) * T5 * T2 / 7# + HeIIEoS(40, M)
End Function


Private Function HelmA(D As Double, TT As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Helmholtz energy as a function of density and temperature [SI units]
' Valid in compressed liquid from 0.8 to about 3 K.
' If TT is less than 0.8 K, TT is assumed to be the isochoric distance
'    to the lambda line (negative in HeII, positive in HeI).
   Dim PVTCSA(6) As Double
   Dim T As Double, DT As Double, V As Double, X As Double
   Dim M As Integer
   V = 1000# / D
   X = V - VLamLf
   If (V <> Vsave) Then Call LamDer(V)
   If (TT > 0.7999) Then
      T = TT
      DT = T - TLam
   Else
      DT = TT
      T = TLam + DT
   End If
   If (DT > 0#) Then
      M = 1
   Else
      M = 2
   End If
   Call LogFun(PVTCSA, T, DT, M)
   HelmA = -1000# * (PVTCSA(6) + HelmA2(X, T, M))
End Function


Private Function HelmA2(X As Double, T As Double, M As Integer) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' "background" Helmholtz energy [J/gm] as a function of V-VLamLf [cm3/gm]
'  and T [K].
' M specifies HeI or HeII.  VLamLf is the volume at the lower lambda point.
' Called from HelmA.
   Dim F(7) As Double, Q(5) As Double, A1 As Double
   Dim T2 As Double, T4 As Double, FF As Double
   Dim J As Integer, K As Integer
   T2 = T * T
   T4 = T2 * T2
   F(1) = X
   For J = 1 To 5
      Q(J) = 0#
   Next
   For K = 1 To 6
      FF = F(K) / CDbl(K)
      For J = 1 To 5
         Q(J) = Q(J) + FF * HeIIEoS(6 * J + K - 1, M)
      Next
      F(K + 1) = F(K) * X
   Next
   A1 = Q(3) + T2 * (Q(4) + T2 * Q(5))
   HelmA2 = Q(1) + T2 * (Q(2) + T2 * A1) + HeIIEoS(36, M) * T2 / 2# + _
            HeIIEoS(37, M) * T4 / 12# + HeIIEoS(38, M) * T4 * T2 / 30# + _
            HeIIEoS(39, M) * T4 * T4 / 56# + HeIIEoS(40, M) * T + _
            HeIIEoS(41, M)
End Function


Private Sub LamDer(V As Double)
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
'-----OUTPUT to module variables
' TLam    = Lambda temperature [K] [T76 scale]
' dTLdV   = dT/dV              [(K-gm)/cm3]
' d2TLdV2 = d2T/dV2            [(K-gm2)/cm6]
'-----INPUT
' V = Specific volume [cm3/gm]
'
   Dim X As Double, Y1 As Double, Y2 As Double
   Dim A(5) As Double
   A(1) = 0.091672438
   A(2) = -0.082840336
   A(3) = 0.071832749
   A(4) = 0.04839517
   A(5) = 0.039159012
   X = V - VLamLf
   Y1 = A(3) + X * (A(4) + X * A(5))
   Y2 = 3# * A(3) + X * (4# * A(4) + X * 5# * A(5))
   TLam = TLamL + X * (A(1) + X * (A(2) + X * Y1))
   dTLdV = A(1) + X * (2# * A(2) + X * Y2)
   d2TLdV2 = 2# * A(2) + X * (6# * A(3) + X * (12# * A(4) + X * 20# * A(5)))
End Sub


Private Sub LogFun(PVTCSA() As Double, TT As Double, DTT As Double, _
                   M As Integer)
' (C) Copyright (1994), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' quasi-logarithmic thermodynamic functions
' Output:
' PVTCSA(1) = pressure [MPa]
' PVTCSA(2) = dP/dV [MPa-kg/m3]
' PVTCSA(3) = dP/dT [Mpa/K]
' PVTCSA(4) = Cv [J/kg -K]
' PVTCSA(5) = entropy [J/kg-K]
' PVTCSA(6) = Helmholz energy [J/kg]
' Input:
' T  = Temperature
' DT = T - Tlambda = isochoric distance to the lambda line
' (Note: global parameters must have evaluated prior to this call
'  to LogFun, by a call to subroutine LamDer.)
'
   Const Max As Integer = 9
   Dim Y(Max) As Double, CL5(5) As Double
   Dim T As Double, T2 As Double, DT As Double, Z As Double
   Dim Fact As Double, R0 As Double, CKM As Double
   Dim XP As Double, XT As Double, XD As Double
   Dim XS As Double, XA As Double, R As Double
   Dim EI As Double, EJ As Double
   Dim I As Integer, J As Integer
   Static AP As Double, AD As Double, AT As Double
   Static AC As Double, AR As Double, AA As Double
   Static Tsave As Double, DTsave As Double
   CL5(1) = 1#
   CL5(2) = -4#
   CL5(3) = 12#
   CL5(4) = -24#
   CL5(5) = 24#
   T = TT
   DT = DTT
   If ((DT <> DTsave) Or (T <> Tsave)) Then
' Calculation is skipped if input parameters are unchanged
      DTsave = DT
      Tsave = T
'
' Y(1) = quasi-logarithmic singularity
' Y(i), i>1, = successive indefinite integrals
'
      If (DT = Zero) Then
         Y(1) = -100#
         For I = 2 To Max
            Y(I) = Zero
         Next
      Else
         Z = Abs(DT)
         Y(1) = Log(Z)
         Fact = 1#
         For I = 2 To Max
            If ((Abs(Y(I - 1)) < 1E-25) And (I > 2)) Then
               Y(I) = Zero
            Else
               Y(I) = (DT * Y(I - 1) - DT ^ (I - 1) / Fact) / CDbl(I - 1)
            End If
            Fact = Fact * CDbl(I)
         Next
      End If
'
' Thermodynamic function evaluation
'-----Version 3.3; function rewritten for greater efficiency
'
      T2 = T * T
      R0 = T2 * T2
      AP = HeIIEoS(1, M) * Y(2) + HeIIEoS(2, M) * (T2 * Y(2) - _
           4# * T * Y(3) + 6# * Y(4))
      AT = HeIIEoS(1, M) * Y(1) + HeIIEoS(2, M) * (T2 * Y(1) - _
           2# * T * Y(2) + 2# * Y(3))
      AD = HeIIEoS(1, M) * Y(1) + HeIIEoS(2, M) * (T2 * Y(1) - _
           4# * T * Y(2) + 6# * Y(3))
      AC = HeIIEoS(1, M) * T * Y(1) + HeIIEoS(2, M) * T2 * T * Y(1)
      AR = HeIIEoS(1, M) * Y(2) + HeIIEoS(2, M) * (T2 * Y(2) - _
           2# * T * Y(3) + 2# * Y(4))
      AA = HeIIEoS(1, M) * Y(3) + HeIIEoS(2, M) * (T2 * Y(3) - _
           4# * T * Y(4) + 6# * Y(5))
      For I = 1 To 3
         CKM = HeIIEoS((I + 2), M)
         XP = R0 * Y(I + 1)
         XT = R0 * Y(I)
         XD = R0 * Y(I)
         XS = R0 * Y(I + 1)
         XA = R0 * Y(I + 2)
         R = R0 / T
         For J = 2 To 5
            EI = CL5(J) * R
            EJ = CDbl(J) * EI
            XP = XP + EJ * Y(I + J)
            XT = XT + EI * Y(I + J - 1)
            XD = XD + EJ * Y(I + J - 1)
            XS = XS + EI * Y(I + J)
            XA = XA + EJ * Y(I + J + 1)
            R = R / T
         Next
         AP = AP + CKM * XP
         AT = AT + CKM * XT
         AD = AD + CKM * XD
         AC = AC + CKM * T * R0 * Y(I)
         AR = AR + CKM * XS
         AA = AA + CKM * XA
      Next
      AD = dTLdV * dTLdV * AD - d2TLdV2 * AP
      AP = -dTLdV * AP
      AT = -dTLdV * AT
   End If
' Output
   PVTCSA(1) = AP
   PVTCSA(2) = AD
   PVTCSA(3) = AT
   PVTCSA(4) = AC
   PVTCSA(5) = AR
   PVTCSA(6) = AA
End Sub


Private Sub Deriv(Cp As Double, CCvo As Double, Gamma As Double, _
                  Alpha As Double, Grun As Double, Kt As Double, _
                  C As Double, JT As Double, xdPdD As Double, _
                  xdPdT As Double, VdHdV As Double, _
                  DI As Double, TI As Double)
' Converted to Visual Basic by Jay Theilacker (1998)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
'
'-----OUTPUT
' Cp    = Cp    = Specific heat at constant P    [J/(kG-K)]
' CCvo  = Cv    = Specific heat at constant V    [J/(kG-K)]
' Gamma = Gamma = Cp/Cv                          [-]
' Alpha = Alpha = (T/V)(dV/dT)  at constant P    [-]
' Grun  = Grun  = (V/Cv)(dP/dT) at constant V    [-]
' Kt    = Kt    = (1/D)(dD/dP)  at constant T    [1/Pa]
' C     = C     = velocity of sound              [m/s]
' JT    = JT    = Joule-Thomson coefficient      [K/Pa]
'                 (= dT/dP at constant H)
' xdPdD = dPdD  = dP/dD at constant T            [Pa-m3/kG]
' xdPdT = dPdT  = dP/dT at constant D            [Pa/K]
' VdHdV = V*(DH/DV) at constant P                [J/kG]
'-----INPUT
' D     = DI = Density                           [kG/M3]
' T     = TI = Temperature                       [K]
'
   Static D As Double, T As Double, A As Double
   Static B As Double, CCv As Double
' A, B, CCV, D, and T are local variables (Saved from the previous call)
'   If ((DI <> D) Or (TI <> T)) Then
      D = DI
      T = TI
      A = dPdT(D, T)
      B = dPdD(D, T)
      CCv = Cv(D, T)
'   End If
   CCvo = CCv
   Kt = 1# / (B * D)
   Alpha = T * A * Kt
   Grun = A / (D * CCv)
   Gamma = 1# + Alpha * Grun
   C = Max((Gamma * B), 1#)
   VdHdV = C / Grun
   C = Sqr(C)
   Cp = Gamma * CCv
   JT = (Alpha - 1#) / (D * Cp)
   xdPdD = B
   xdPdT = A
End Sub


Private Sub SuperF(RSR As Double, V2S As Double, V4S As Double, _
                   GM As Double, Cond As Double, _
                   PI As Double, D As Double, T As Double, _
                   DTV As Double, S As Double, _
                   Cp As Double, CCvo As Double, Gamma As Double, _
                   Alpha As Double, Grun As Double, Kt As Double, _
                   C As Double, JT As Double, xdPdD As Double, _
                   xdPdT As Double, VdHdV As Double)
' (C) Copyright (1990, 1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'  This subroutine returns parameters from the two-fluid model
'  of superfluidity:
'  RSR  = superfluid density fraction                   [-]
'  V2S  = velocity of second sound                      [m/s]
'  V4S  = velocity of fourth sound                      [m/s]
'  GM   = Gorter-Mellink mutual friction parameter      [m-s/kg]
'  Cond = conductivity parameter                        [(W/m2)3/(m/K)]
'         (The defining equation for COND is that used by
'          Srinivasan and Hofmann, Cryogenics 25, 652 (1985), i.e.,
'          Q**3 = Cond * (dT/dX)
'                 where Q     = heat flux               [W/m2]
'                       dT/dX = temperature gradient    [K/m]
'          ! This equation assumes turbulent counterflow.
'-----INPUT
'  P, D, T, S are pressure, density,  temperature & entropy (SI units)
'  (see subroutine Deriv for definitions of last 11 variables).
'  DTV is T - Tlambda(V), ie, the isochoric distance from Tlambda
'  Version July 12, 1992; for HEPAK v3.2
   Dim DS As Double, DN As Double, C2 As Double
   Dim ZZ As Double
' Validity very close to the lambda line is unknown.
   If (DTV < -0.0009) Then
      RSR = RSRfPT(PI, T)
      DS = D * RSR
      DN = D - DS
'  Maynard's equations (2), (3), and (5); assume small corrections.
      C2 = C * C
      V2S = (D / DN - 1#) * T * S * S / Cp
      V2S = V2S * (1# + Alpha * Grun * V2S / (C2 - V2S))
      ZZ = (V2S / C2) * (1# - 2# * Alpha * C2 / (Gamma * S * T))
      V4S = C2 * (1# - (DN / D) * (1# - ZZ))
      V2S = Sqr(V2S)
      V4S = Sqr(V4S)
'  GM is fitted to the data of Srinivasan and Hofmann
      GM = 383.7 + 10.649 * DN
'  Asymtotic form approaching TL is given by Frederking
      If (GM > (1413# + 12.4 * (D - 160#))) Then GM = 558.9 / RSR
      Cond = S * ((DS * S * T) ^ 3) / (GM * DN)
   Else
      RSR = 0#
      V2S = 0#
      V4S = 0#
      GM = 0#
      Cond = 0#
   End If
End Sub


Private Function RSRfPT(P As Double, T As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Rhos/Rho = superfluid density fraction, from Maynard, PR B14, 3868 (1976)
' as a function of P [Pascal] and T [K].
' Maynard's data includes the range T/Tlambda(P) > 0.55
   Dim A(3) As Double, X As Double
   A(1) = -2.08403569
   A(2) = 1.760235312
   A(3) = -1.469764627
   X = 1# - T / TLfP(P)
   If (X <= 0#) Then
      RSRfPT = 0#
   ElseIf (X <= 0.45) Then
      RSRfPT = A(1) * X * Log(X) + X * X * (A(2) + X * A(3))
   ElseIf (X < 1#) Then
' The following is an approximate empirical fit beyond Maynard's data
      RSRfPT = 1# - 0.02863 * ((1# - X) / 0.55) ^ 6
   Else
      RSRfPT = 1#
   End If
End Function


Private Sub AMLap(TminM As Double, TmaxA As Double, D As Double)
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Computes the overlap temperatures between Arp and McCarty equations
' as a function of density [kg/m3].
' TminM is the minimum temperature for McCarty (in compressed liquid).
' TmaxA is the maximum temperature for Arp
   Const TM As Double = 2.53
   Const TA As Double = 2.98
   Const Dmin As Double = 140#
   Const Dmax As Double = 190#
   Dim DT As Double
   If ((D < Dmin) Or (D > Dmax)) Then
      TminM = 0#
      TmaxA = 0.1
   Else
'     Note: -0.28 [K] / (Dmax-Dmin) = -0.0056 = slope of boundary line
      DT = -0.0056 * (D - Dmin)
      TminM = TM + DT
      TmaxA = TA + DT
      If (D > 180#) Then
         DT = -0.035 * (D - 180#)
         TminM = TminM + DT
         TmaxA = TmaxA + DT
      End If
   End If
End Sub


Private Function D2LfPT(xPsat As Double, Tsat As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Saturated liquid density, 0.8 to 5.1953 K, T76 scale.
' Output is the density such that Press(D,T) = xPsat; will differ slightly
' from the density DfSat, except when T > TNRC.
' Valid input of both xPsat and Tsat must be made by the calling program.
' Version Nov. 8, 1990
   Const P0 As Double = 5000#
   Dim DL As Double, D As Double, DPD As Double, DelP As Double
   Dim K As Integer
   DL = DFsat(Tsat)
   If (Tsat < TNRC) Then
      D = DL
      DPD = dPdD(D, Tsat)
      For K = 1 To 8
         DelP = Press(D, Tsat) - xPsat
         D = D - DelP / DPD
         If (Abs(DelP / (xPsat + P0)) < E(12)) Then Exit For
      Next
      If (Tsat > Tlow) Then
         DL = (D * (TNRC - Tsat) + DL * (Tsat - Tlow)) / DelT
      Else
         DL = D
      End If
   End If
   D2LfPT = DL
End Function


Private Function DFsat(T As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Density of saturated liquid [kg/m3]; 0.8 < T < 5.1953; T76 scale.
' This is the best independent estimate of saturation density;
' However, it is not exactly consistent with function D2LfPT.
' HeI range above TNRC:  theoretical form, fitted for continuity at TNRC
' HeI range above 3.2 K: R.D.McCarty equation May 4, 1988.
' HeI range below 3.2 K: Van Degrift' s equation, quoted by Barenghi
'     Lucas and Donnelly, J. Low Temp. Physics 44, 491 (1981),
'     shifted to the T76 temperature scale.
'
' HeII range: V.Arp equation Dec 16, 1987.
' Version Nov. 8, 1990
   Const TBB As Double = 3.1853
   Const DTPL As Double = 36.514
   Const Den0 As Double = 145.188
   Const TNRC As Double = 5.106
   Const DNRC As Double = 93.3275
   Dim Dcc As Double, D As Double, P As Double
   Dim X As Double, Y As Double, Z As Double
   Dim J As Integer
   Dim A(6) As Double, B(4) As Double, C(5) As Double
   A(1) = -34.34500882
   A(2) = 0.4501470666
   A(3) = -6.89546017
   A(4) = 54.42223002
   A(5) = -97.16273774
   A(6) = 17.85061367
   B(1) = -1.29782247
   B(2) = -4.911239491
   B(3) = 1.021276082
   B(4) = -2.626058542
   C(1) = -243.3560558
   C(2) = 1814.476908
   C(3) = 2425.346403
   C(4) = -460.4869802
   C(5) = 245.3561962
   Dcc = Dcrit / GMWT
   If (T >= TBB) Then
      If (T > TNRC) Then
         P = (Max(1E-20, ((Tcrit - T) / (Tcrit - TNRC)))) ^ One3RD
         D = Dcc + P * (DNRC / GMWT - Dcc)
      Else
         X = (T - Tcrit) / (TLamL - Tcrit)
         Y = X ^ One3RD
         Z = 1# / X
         P = A(1) * Log(X)
         For J = 2 To 6
            P = P + A(J) * (1# - Z)
            If (J = 4) Then
               Z = Y
            Else
               Z = Z * Y
            End If
         Next
         D = Dcc + Exp(P) * (DTPL - Dcc)
      End If
      DFsat = D * GMWT
   Else
      Z = T - TLamL
      If (Z >= 0#) Then
         DFsat = DLamLf + Z * (B(2) + Z * (B(3) + Z * B(4)))
         If (Z <> 0#) Then DFsat = DFsat + B(1) * Z * Log(Z)
      Else
         Y = Z * (Log(-Z) - 1#)
         P = T * T
         DFsat = ((P * P * (C(1) * Y + C(2) * Z + C(3) + C(4) * P) + _
                 C(5) * P) * 0.000001 + 1#) * Den0
      End If
   End If
End Function


Private Function D2VfPT(xPsat As Double, Tsat As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Saturated vapor density, 0.8 to 5.1953 K, T76 scale.
' Output is the density such that Press(D,T) = xPsat; may differ slightly
' from the density DGsat, except when T > TNRC.
' Valid input of both xPsat and Tsat must be made by the calling program.
' Version Nov. 6, 1990
   Const P0 As Double = 5000#
   Dim K As Integer
   Dim DV As Double, D As Double, DPD As Double, DelP As Double
   DV = DGsat(Tsat)
   If (Tsat < TNRC) Then
      D = DV
      DPD = dPdD(D, Tsat)
      For K = 1 To 8
         DelP = Press(D, Tsat) - xPsat
         D = D - DelP / DPD
         If (Abs(DelP / (xPsat + P0)) < E(12)) Then Exit For
      Next
      If (Tsat > Tlow) Then
         DV = (D * (TNRC - Tsat) + DV * (Tsat - Tlow)) / DelT
      Else
         DV = D
      End If
   End If
   D2VfPT = DV
End Function


Private Function DGsat(T As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Saturated vapor density, 0.8 < T < 5.1953.
' This is the best independent estimate of saturation density;
' Theoretical form above TNRC, fitted for continuity at TNRC.
' Version Nov. 8, 1990
   Const DTPV As Double = 0.294058864
   Const TNRC As Double = 5.106
   Const DNRC As Double = 46.5087
   Dim P As Double, D As Double, Dcc As Double
   Dim X As Double, Y As Double, Z As Double
   Dim J As Integer
   Dim A(8) As Double, C(8) As Double
   A(1) = -561.8978079
   A(2) = 3.300736785
   A(3) = -60.31200561
   A(4) = 612.990156
   A(5) = -2718.577178
   A(6) = 1285.18502
   A(7) = -440.6873907
   A(8) = 71.63145577
   C(1) = -7.41816
   C(2) = 5.42128
   C(3) = 9.903203
   C(4) = -9.617095
   C(5) = 6.804602
   C(6) = -3.0154606
   C(7) = 0.7461357
   C(8) = -0.0791791
   Dcc = Dcrit / GMWT
   If (T > TLamL) Then
      If (T > TNRC) Then
         P = (Max(1E-20, ((Tcrit - T) / (Tcrit - TNRC)))) ^ One3RD
         D = Dcc + P * (DNRC / GMWT - Dcc)
      Else
         X = (T - Tcrit) / (TLamL - Tcrit)
         Y = X ^ One3RD
         Z = 1# / X
         P = A(1) * Log(X)
         For J = 2 To 8
            P = P + A(J) * (1# - Z)
            If (J = 4) Then
               Z = Y
            Else
               Z = Z * Y
            End If
         Next
         D = Dcc + Exp(P) * (DTPV - Dcc)
      End If
      DGsat = D * GMWT
   ElseIf (T > 0.799) Then
      P = C(1) / T + C(2)
      Y = 1#
      For J = 3 To 8
         Y = Y * T
         P = P + C(J) * Y
      Next
      D = Exp(P) / (2077.2258 * T)
      Z = (0.0537 - 0.514 / T) / 4.0026
      DGsat = D / (1# + Z * D) + 0.00164 * D ^ 3
   Else
      DGsat = -1#
   End If
End Function


Private Sub Kierst(P As Double, xdPdT As Double, D As Double, _
                   dDdT As Double, T As Double)
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
'-----OUTPUT
' P    = Lambda-line pressure          [Pa]
'      = 5041.8 at T=2.1768
' xdPdT= dP/dT along the lambda line   [Pa/K]
' D    = Density at the lambda line    [kG/M3]
'      = 146.15 at T=2.1768
' dDdT = dD/dT along the lambda line   [kG/(M3-K)]
'-----INPUT
' T    = temperature  (T76 scale)      [K]
'        valid range is 2.1768 to 1.7673
'-----REFERENCE
' H.A. Kierstead, Phys Rev 162, 153 (1967); constants revised for T76 scale
' Note: Kierstead's density equation disagrees slightly with HEPAK,
' and is not used by HEPAK.
'-----
   Dim X As Double, Z As Double, Y1 As Double, Y2 As Double
   Dim A(7) As Double, B(7) As Double
   A(1) = 0.42774167
   A(2) = -94.820469
   A(3) = -85.817089
   A(4) = -102.39597
   A(5) = -76.73524
   A(6) = -0.37798315
   A(7) = 42.148155
   B(1) = 0.14841216
   B(2) = -0.15036724
   B(3) = -0.32811465
   B(4) = -0.52635312
   B(5) = -0.37937084
   B(6) = -0.00226216
   B(7) = 36.64523
   X = T - TLamL
   Z = Exp(A(7) * X)
   Y1 = A(3) + X * (A(4) + X * A(5))
   Y2 = B(3) + X * (B(4) + X * B(5))
   P = (A(1) + X * (A(2) + X * Y1) + A(6) * Z) * 101325#
   xdPdT = (A(2) + X * (2# * A(3) + X * (3# * A(4) + X * 4# * A(5))) + _
            A(6) * A(7) * Z) * 101325#
   Z = Exp(B(7) * X)
   D = (B(1) + X * (B(2) + X * Y2) + B(6) * Z) * 1000#
   dDdT = (B(2) + X * (2# * B(3) + X * (3# * B(4) + X * 4# * B(5))) + _
           B(6) * B(7) * Z) * 1000#
End Sub


Private Function D175K(Pascal As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' density [kg/m3] of HeII as a function of P [Pa] at T(58)=1.75
' accuracy probably about 0.2%.
   Dim C(6) As Double, P As Double
   Dim A1 As Double, A2 As Double
   C(1) = -0.8129582813
   C(2) = 18.5990126
   C(3) = -6.344068036
   C(4) = 2.844939685
   C(5) = -0.8216755033
   C(6) = 0.1044392629
   Const Den0 As Double = 146.081
   P = 0.000001 * Pascal
   A1 = C(4) + P * (C(5) + P * C(6))
   A2 = C(2) + P * (C(3) + P * A1)
   D175K = Den0 + C(1) + P * A2
End Function


Private Function PLfT(T As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' P = Lambda-line pressure = 5041.8 Pa    at T = 2.1768 K
'                          = 30.134E+5 Pa at T = 1.7673 K
' reference:
' H.A. Kierstead, Phys Rev 162, 153 (1967) (T58 scale)
' refitted to T76 scale.
   Dim B(7) As Double, X As Double, P As Double
   B(1) = 0.42774167
   B(2) = -94.820469
   B(3) = -85.817089
   B(4) = -102.39597
   B(5) = -76.73524
   B(6) = -0.37798315
   B(7) = 42.148155
   X = T - TLamL
   P = B(1) + (B(2) + (B(3) + (B(4) + B(5) * X) * X) * X) * X + _
       B(6) * Exp(B(7) * X)
   PLfT = P * 101325#
End Function


Private Function PMfT(TT As Double) As Double
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
'  melting pressure as a function of temperature, T76 scale below 5.1953
'  F1 to F3 have been refitted to T58 equations of Grilly and Mills.
'  F4 and F5 have been fitted to Grilly and Mills data.
   Const DelTT As Double = 0.005
   Dim T As Double, W As Double, X As Double
   T = TT
' Cut off at the upper limit of HEPAK
   If (T > 14.0104) Then
      PMfT = 101325000#
   ElseIf (T >= Tcrit + DelTT) Then
      PMfT = F1PMfT(T)
   ElseIf (T > Tcrit - DelTT) Then
      W = (T - Tcrit) / (2# * DelTT) + 0.5
      PMfT = W * F1PMfT(T) + (1# - W) * F2PMfT(T)
   ElseIf (T >= 2.0044 + DelTT) Then
      PMfT = F2PMfT(T)
   ElseIf (T > 2.0044 - DelTT) Then
      W = (T - 2.0044) / (2# * DelTT) + 0.5
      PMfT = W * F2PMfT(T) + (1# - W) * F3PMfT(T)
   ElseIf (T >= 1.766 + DelTT) Then
      PMfT = F3PMfT(T)
   ElseIf (T > 1.766 - DelTT) Then
      W = (T - 1.766) / (2# * DelTT) + 0.5
      PMfT = W * F3PMfT(T) + (1# - W) * F4PMfT(T)
   ElseIf (T >= 1.4676 + DelTT) Then
      PMfT = F4PMfT(T)
   ElseIf (T > 1.4676 - DelTT) Then
      W = (T - 1.4676) / (2# * DelTT) + 0.5
      PMfT = W * F4PMfT(T) + (1# - W) * F5PMfT(T)
   ElseIf (T >= Tmin) Then
      PMfT = F5PMfT(T)
   Else
      PMfT = -1#
   End If
End Function


Private Function F1PMfT(X As Double) As Double
   Const Conv As Double = 98066.5
   F1PMfT = (-17.8 + 17.31457 * X ^ 1.555414) * Conv
End Function


Private Function F2PMfT(X As Double) As Double
   Const Conv As Double = 98066.5
   Dim A1 As Double
   A1 = 32.26926 + X * (-4.91282 + X * 0.310795)
   F2PMfT = (34.2097 + X * (-45.31231 + X * A1)) * Conv
End Function


Private Function F3PMfT(X As Double) As Double
   Const Con2 As Double = 101325#
   F3PMfT = (17.8537 + X * (-15.49444 + X * 12.57562)) * Con2
End Function


Private Function F4PMfT(X As Double) As Double
   Const Con2 As Double = 101325#
   F4PMfT = (99.3328 + X * (-101.4497 + X * 35.1175)) * Con2
End Function


Private Function F5PMfT(X As Double) As Double
   Const Con2 As Double = 101325#
   F5PMfT = (24.997 + 3.3693 * (X - 0.725) ^ 4) * Con2
End Function


Private Sub PsatfT(P As Double, dPdTs As Double, T As Double)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
'  saturation pressure and dP/dT as a function of temperature
'  on the T76 scale, from 0.5 to 5.1953 K; Equation is given by
'  Durieux and Rusby, Metrologia 19, 67 (1983).
'
' Note: c(11,1) has been corrected in this version, August 19, 1992.
' Earlier versions of this routine had c(11,1)=14.5333d0 (missing a "3");
' the corresponding error in pressure is never more than 0.01 pct.
'
   Dim M As Integer, MX As Integer, J As Integer
   Dim Q0 As Double, Q1 As Double, TN As Double
   Dim X As Double, Y As Double
   Dim C(11, 2) As Double
   C(1, 1) = -30.93285
   C(2, 1) = 392.47361
   C(3, 1) = -2328.04587
   C(4, 1) = 8111.30347
   C(5, 1) = -17809.80901
   C(6, 1) = 25766.52747
   C(7, 1) = -24601.4
   C(8, 1) = 14944.65142
   C(9, 1) = -5240.36518
   C(10, 1) = 807.93168
   C(11, 1) = 14.53333
   C(1, 2) = -7.41816
   C(2, 2) = 5.42128
   C(3, 2) = 9.903203
   C(4, 2) = -9.617095
   C(5, 2) = 6.804602
   C(6, 2) = -3.0154606
   C(7, 2) = 0.7461357
   C(8, 2) = -0.0791791
   C(9, 2) = 0#
   C(10, 2) = 0#
   C(11, 2) = 0#
' at TL, P=5041.8 and dPdT=12407.9; fixed points on the T76 scale
' Note: Durieux and Rusby state Pc=227.46 kPa at 5.1953 K
' However, this (corrected) equation gives Pc=227.4623 kPa at 5.1953 K.
   If ((T > Tcrit + 0.00001) Or (T < Tmin)) Then
      P = -1#
      dPdTs = -1#
      Exit Sub
   ElseIf (T > TLamL) Then
      X = Min(T / Tcrit, 1#)
      M = 1
      MX = 10
   Else
      X = T
      M = 2
      MX = 8
   End If
   Q0 = C(1, M) / X + C(2, M)
   Q1 = -C(1, M) / (X * X)
   TN = 1#
   For J = 3 To MX
      Q1 = Q1 + C(J, M) * TN * CDbl(J - 2)
      TN = TN * X
      Q0 = Q0 + C(J, M) * TN
   Next
   X = 1# - X
   If ((M = 1) And (X > 0#)) Then
      Y = X ^ 0.9
      Q0 = Q0 + X * C(11, M) * Y
      Q1 = Q1 - 1.9 * C(11, M) * Y
   End If
   P = Exp(Q0)
   dPdTs = P * Q1
   If (M = 1) Then dPdTs = dPdTs / Tcrit

End Sub


Private Function Psat(T As Double) As Double
' Converted to Visual Basic by Jay Theilacker (1998)
   Dim TT As Double, P As Double, xdPdT As Double
   TT = T
   Call PsatfT(P, xdPdT, TT)
   Psat = P
End Function


Private Function DLfT(TT As Double) As Double
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
'
' Lambda line density as a function of temperature, T76 scale
' valid range is 2.1768 to 1.7673
   Dim D As Double, T As Double, JX As Integer
   T = TT
   Call RootAZ(D, 146.14, 147#, 155#, 182#, E(8), E(9), _
               E(11), E(13), 1, T, JX)
   If (JX <= 0) Then
      DLfT = -1#
   Else
      DLfT = D
   End If
End Function


Private Function DLfP(P As Double) As Double
'  lambda line density [kg/m3] as a function of pressure [Pa]
   Dim T As Double
   T = TLfP(P)
   DLfP = DLfT(T)
End Function


Private Function DMfP(P As Double) As Double
'-----melting density as a function of pressure
   Dim T As Double
   T = TMfP(P)
   DMfP = DMfT(T)
End Function


Private Function DMfT(T As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
'
' Liquid density at the melting line [kg/m3] as a function of T [K]
' Range 0.8 to 14.0 K; accuracy generally better than 0.3 kg/m3
' Version March 22, 1989.
   Dim Clow(3) As Double, Cmed(3) As Double, Chigh(4) As Double
   Dim Z As Double
   Const T1 As Double = 1.7673
   Const T2 As Double = 3.53
   Clow(1) = 167.9361577
   Clow(2) = 6.728283584
   Clow(3) = -10.18218341
   Cmed(1) = 146.9940814
   Cmed(2) = 18.58136626
   Cmed(3) = -1.497696476
   Chigh(1) = 154.2809083
   Chigh(2) = 18.90207406
   Chigh(3) = -0.8746341757
   Chigh(4) = 0.02235656147
   If (T > T2) Then
      DMfT = Chigh(1) + T * (Chigh(2) + T * (Chigh(3) + T * Chigh(4)))
   Else
      Z = T - T1
      If (Abs(Z) < 0.0001) Then
         DMfT = 0#
      Else
         DMfT = Z * Log(Abs(Z))
      End If
      If (Z >= 0#) Then
         DMfT = Cmed(3) * DMfT + Cmed(1) + T * Cmed(2)
      Else
         DMfT = Clow(3) * DMfT + Clow(1) + T * Clow(2)
      End If
   End If
End Function


Private Function TLfD(D As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' lambda line temperature as a function of density
   Dim Vcgs As Double
   Vcgs = 1000# / D
   Call LamDer(Vcgs)
   TLfD = TLam
End Function


Private Function TLfP(P As Double) As Double
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' lambda line temperature as a function of pressure [Pa]
   Dim PP As Double, TT As Double, JX As Integer
   PP = P
' T=2.17725 corresponds to P=0
   Call RootAZ(TT, 2.17725, 2.1768, 2#, 1.7672, _
               E(9), 0#, E(2), E(11), 2, PP, JX)
   TLfP = TT
End Function


Private Function TMfD(DD As Double) As Single
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
   Dim D As Double, X As Double, TA As Double, TB As Double
   Dim T As Double, JX As Integer
   D = DD
   X = (D - 171.2) / 136#
   TA = Max(1#, 0.8 + 13.2 * X)
   TB = 2# + Min(TA, 12#)
   TA = TB - 2#
   Call RootAZ(T, 0.8, TA, TB, 14.011, E(9), 0#, E(5), E(10), _
               3, D, JX)
   If (JX < 0) Then T = 0.8
   TMfD = T
End Function
'
'
'
Private Function TMfP(PP As Double) As Double
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' melting temperature as a function of pressure [Pa]
   Dim P As Double, TT As Double, JX As Integer
   P = PP
   Call RootAZ(TT, 0.8, 2#, 8#, 14.0104, _
               E(8), E(10), E(2), E(10), 4, P, JX)
   If (JX <= 0) Then
      TMfP = -1#
   Else
      TMfP = TT
   End If
End Function


Private Function TsatfP(PP As Double) As Double
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' saturation temperature as a function of pressure
   Dim P As Double, PL As Double, JX As Integer
   Dim T As Double, TB As Double, Tc As Double
   P = Max(PP, 1.475154)
   PL = Log(P)
   TB = 0.6985201046 + PL * (0.2654401757 + PL * _
        (-0.05434555937 + PL * 0.005045477179))
   TB = Max(Tmin, TB - 0.04)
   Tc = Min(Tcrit, TB + 0.08)
   Call RootAZ(T, Tmin, TB, Tc, Tcrit, _
               E(10), E(13), E(3), E(10), 5, P, JX)
   If (JX <= 0) Then
      TsatfP = -1#
   Else
      TsatfP = T
   End If
End Function


Private Function DenMax(T As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
'-----OUTPUT
' DenMax = Density [kG/m3] ; maximum density for which
'          the equations are valid (melting pressure to 14 K;
'          1000 atmos to 75 K, 20000 bars to 300 K, 1000 atmos to 1500 K)
'-----INPUT
' T  = Temperature [K] ; 0.8 < T < 1500.
' Accuracy of fit: better than 0.3 kg/m3 for all T
'-----
   Dim A(4) As Double, B(4) As Double, C(3) As Double
   Dim D As Double, X As Double
   A(1) = -16.70040531
   A(2) = 485.7565167
   A(3) = -258.5274317
   A(4) = 51.61716546
   B(1) = 292.4157352
   B(2) = 1202.25329
   B(3) = -1859.861872
   B(4) = 1060.045153
   C(1) = 0.9714767056
   C(2) = 474.5050407
   C(3) = -302.8968312
   If (T < 14.05) Then
      D = DMfT(T)
   Else
      X = 100# / (T + 50#)
      If (T < 74.99) Then
         D = A(1) + X * (A(2) + X * (A(3) + X * A(4)))
      ElseIf (T < 300.01) Then
         D = B(1) + X * (B(2) + X * (B(3) + X * B(4)))
      Else
         D = C(1) + X * (C(2) + X * C(3))
      End If
   End If
' Add a little for error tolerance.
   DenMax = D + 0.5
End Function


Private Function H3K(P As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
'-----OUTPUT
' H3K = Enthalpy [J/kG] at 3 K, for P > P(sat)
'-----INPUT
' P   = Pressure [Pa] ; fitted to 70 atmospheres
'-----ACCURACY
' Accuracy about +/- 4.  J/kG for pressures below 2 bars
'          about +/- 25. J/kG for pressures above 2 bars
'-----Version Sept 18, 1987
   Dim C(4) As Double, Z As Double
   C(1) = 4954.376389
   C(2) = -130.2167945
   C(3) = 657.1877607
   C(4) = -13.82310121
   If (P <= 0#) Then
      H3K = -1#
   Else
      Z = Sqr(0.00001 * P)
      H3K = C(1) + Z * (C(2) + Z * (C(3) + Z * C(4)))
   End If
End Function


Private Function S3K(P As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
'-----OUTPUT
' S3K = Entropy [J/(kG-K)] at 3 K, for P > P(sat)
'-----INPUT
' P   = Pressure [Pa] ; fitted to 70 atmospheres
'-----ACCURACY
' Accuracy about +/- 1. J/(kG-K) for pressures below 2 bars
'          about +/- 8. J/(kG-K) for pressures above 2 bars
'-----Version Sept 18, 1987
   Dim C(5) As Double, Z As Double, A1 As Double
   C(1) = 2365.990305
   C(2) = -30.12076857
   C(3) = -34.92403215
   C(4) = 5.141720002
   C(5) = -0.261267804
   If (P <= 0#) Then
      S3K = -1#
   Else
      Z = Sqr(0.00001 * P)
      A1 = C(3) + Z * (C(4) + Z * C(5))
      S3K = C(1) + Z * (C(2) + Z * A1)
   End If
End Function


Private Function U3K(P As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' This function written by V. Arp.
' Converted to Visual Basic by Jay Theilacker (1998)
'
'-----OUTPUT
' U3K = Internal Energy [J/kG] at 3 K, for P > P(sat)
'-----INPUT
' P   = Pressure [Pa] ; fitted to 70 atmospheres
'-----ACCURACY
' Accuracy about +/- 4.  J/kG for pressures below 2 bars
'          about +/- 10. J/kG for pressures above 2 bars
'-----Version Sept 18, 1987
   Dim C(5) As Double, Z As Double, A1 As Double
   C(1) = 4947.998498
   C(2) = -100.8309679
   C(3) = -106.6139445
   C(4) = 29.95380841
   C(5) = -1.53619894
   If (P <= 0#) Then
      U3K = -1#
   Else
      Z = Sqr(0.00001 * P)
      A1 = C(3) + Z * (C(4) + Z * C(5))
      U3K = C(1) + Z * (C(2) + Z * A1)
   End If
End Function


Private Function D3K(P As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
'-----OUTPUT
' D3K = Density [kG/m3] at 3 K, for P > P(sat)
'-----INPUT
' P   = Pressure [Pa] ; fitted to 74 bars
'-----ACCURACY
' Accuracy about +/- 0.1 kG/m3 for pressures below 15 bars
'          about +/- 0.3 kG/m3 for pressures above 15 bars
'-----Version Sept 18, 1987
   Dim C(5) As Double, Z As Double, A1 As Double
   C(1) = 139.8535389
   C(2) = 1.496220385
   C(3) = 2.334289084
   C(4) = -0.3419086874
   C(5) = 0.0171614572
   If (P <= 0#) Then
      D3K = -1#
   Else
      Z = Sqr(0.00001 * P)
      A1 = C(3) + Z * (C(4) + Z * C(5))
      D3K = C(1) + Z * (C(2) + Z * A1)
   End If
End Function


Private Sub SatfS(IDID As Integer, Qual As Double, PS As Double, _
                  xdPdT As Double, TS As Double, DL As Double, _
                  DV As Double, SS As Double)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Quality, saturated pressure, temperature, and densities as a function
' of (saturated) entropy.  SI units.  IDID=2 when successful.
'     Range 0.8 to 5.1953 K
'-----Version Oct. 17, 1992
   Dim S As Double, TX As Double, TA As Double, TB As Double
   Dim TT As Double, JX As Integer
   S = SS
   If (S >= Scrit) Then
      JR = 2
      TX = TsatSv(S)
   Else
      JR = 1
      TX = TsatSL(S)
   End If
   TA = Min(Max(TX, 0.81), 5.12)
   TB = TA + 0.06
   Call RootAZ(TT, Tmin, TA, TB, Tcrit, E(9), 0#, _
               E(5), E(8), 6, S, JX)
   If (JX < 0) Then
      IDID = -124
   Else
      IDID = 2
      TS = TT
      PS = SubRT(1)
      xdPdT = SubRT(2)
      If (JR = 2) Then
         DV = SubRT(3)
         DL = D2LfPT(PS, TS)
         Qual = 1#
      Else
         DL = SubRT(3)
         DV = D2VfPT(PS, TS)
         Qual = 0#
      End If
   End If
End Sub


Private Function SatS(T As Double) As Double
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' function used by SatfS
   Dim P As Double, xdPdT As Double, D As Double
   Call PsatfT(P, xdPdT, T)
   If (JR = 1) Then
      D = D2LfPT(P, T)
   Else
      D = D2VfPT(P, T)
   End If
   SatS = Entrop(D, T)
   SubRT(1) = P
   SubRT(2) = xdPdT
   SubRT(3) = D
End Function


Private Function TsatSL(SJKG As Double) As Single
' Converted to Visual Basic by Jay Theilacker (1998)
' Saturated liquid Temperature [K] as a function of entropy [J/kg-K]
' Valid for 4.7128 < SJKG < 5768.4 J/kg-K, corresponding to
'           0.8 < T < 5.1953 K
' Fitted to HEPAK data, Oct 17, 1992; rms accuracy (+-) 0.035 K
   Dim C(6) As Double, S As Double, S2 As Double, A1 As Double
   C(1) = 2.017414715
   C(2) = -1.192464994
   C(3) = 2.441002711
   C(4) = -1.804761829
   C(5) = 0.580036035
   C(6) = -0.06892288428
   Const Z As Double = 0.1736
   S = Max((SJKG * 0.001), 0.0047)
   S2 = Sqr(S)
   A1 = C(4) + S2 * (C(5) + S2 * C(6))
   TsatSL = C(1) * (S ^ Z) + S ^ 4 * (C(2) + S2 * (C(3) + _
            S2 * A1))
End Function


Private Function TsatSv(SJKG As Double) As Single
' Converted to Visual Basic by Jay Theilacker (1998)
' Saturated vapor Temperature [K] as a function of entropy [J/kg-K]
' Valid for 5768.4 < SJKG < 23936 J/kg-K, corresponding to
'           5.1953 > T > 0.8 K
' Fitted to HEPAK data, Oct 17, 1992; rms accuracy (+-) 0.035 K
   Dim C(6) As Double, S As Double, A1 As Double
   C(1) = -26.21235672
   C(2) = 24.31127551
   C(3) = -13.7770157
   C(4) = 3.157246354
   C(5) = -0.3331317936
   C(6) = 0.01342914271
   Const Z As Double = 0.6
   S = (SJKG * 0.001) ^ Z
   A1 = C(4) + S * (C(5) + S * C(6))
   TsatSv = C(1) + S * S * (C(2) + S * (C(3) + S * A1))
End Function


Private Sub SatfD(IDID As Integer, Qual As Double, xPsat As Double, _
                  xdPdT As Double, Tsat As Double, DL As Double, _
                  DV As Double, Rho As Double)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' This subroutine written by V. Arp.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Output: saturation properties: temperature, pressure, dP/dT,
'         liquid and vapor densities (one of these will be a repeat
'         of the input parameter, Rho, all SI units)
' Input:  Rho = saturated fluid density.
' NOTE: The maximum saturated liquid density of 146.1603 occurs at
'       T = 2.1852 K, according the HEPAK v3.2.  In the liquid, SatfD
'       always returns a temperature T > 2.1852; the saturated superfluid
'       branch is not found in the following iteration.
' Version August 29, 1992
   Const DelTB As Double = 0.0893
   Const DNRCL As Double = 93.3275
   Const DNRCV As Double = 46.5087
   Dim TB As Double, Tc As Double, Dref As Double, Fact As Double
   Dim Term As Double, T As Double, JX As Integer
   Dim TA1 As Double, TB1 As Double
'  DelTB = TB - TNRC;  TNRC is defined in DFSAT and DGSAT
   If ((Rho > 146.1603) Or (Rho < 0.000888)) Then
      IDID = -123
      Exit Sub
   End If
   If (Rho > DNRCL) Then
      JR = 0
      TB = 2.1852
      Tc = 4.5
   ElseIf (Rho > DNRCV) Then
      If (Rho >= Dcrit) Then
         Dref = DNRCL
'  Fact = (DNRCV - Dcrit) / (DNRCL - Dcrit)
         Fact = -0.9766194
         Qual = 0#
      Else
         Dref = DNRCV
'  Fact = (DNRCL - Dcrit) / (DREF - Dcrit)
         Fact = -1.0239403
         Qual = 1#
      End If
      Term = Max(Abs((Rho - Dcrit) / (Dref - Dcrit)), 0.00000001)
      Tsat = Tcrit - DelTB * (Term ^ 3)
      Call PsatfT(xPsat, xdPdT, Tsat)
      Term = Dcrit + Fact * (Rho - Dcrit)
      If (Qual < 0.5) Then
         DL = Rho
         DV = Term
      Else
         DL = Term
         DV = Rho
      End If
      IDID = 2
      Exit Sub
   Else
      JR = 1
      TB = Tmin + 0.06 * Rho
      Tc = 0.5 * (TB + Tcrit)
   End If
   Call RootAZ(T, Tmin, TB, Tc, Tcrit, 0.03, 0#, _
               0.00000001, 0.0003, 7, Rho, JX)
   If (JX <= 0) Then
      IDID = -123
   Else
      TA1 = Min(T, Tcrit - 0.005)
      TB1 = Min(T + 0.005, Tcrit)
      Call RootAZ(T, Tmin, TA1, TB1, Tcrit, E(8), E(10), _
                  E(15), E(9), 8, Rho, JX)
      If (JX < 0) Then
         IDID = -180
         Exit Sub
      End If
      IDID = 2
      xPsat = SubRT(1)
      xdPdT = SubRT(2)
      Tsat = T
      If (JR = 0) Then
         DL = Rho
         DV = D2VfPT(xPsat, Tsat)
         Qual = 0#
      Else
         DV = Rho
         DL = D2LfPT(xPsat, Tsat)
         Qual = 1#
      End If
   End If
End Sub


Private Function SatD(T As Double) As Double
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' function used by SatfD
   Call PsatfT(SubRT(1), SubRT(2), T)
   If (JR = 0) Then
      SatD = D2LfPT(SubRT(1), T)
   Else
      SatD = D2VfPT(SubRT(1), T)
   End If
End Function


Private Function SatD1(T As Double) As Double
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' function used by SatfD
   If (JR = 0) Then
      SatD1 = DFsat(T)
   Else
      SatD1 = DGsat(T)
   End If
End Function


Private Sub SatfY(IDID As Integer, Qual As Double, PS As Double, _
                  xdPdT As Double, TS As Double, DL As Double, _
                  DV As Double, YY As Double, Label As String)
' Converted to Visual Basic by Jay Theilacker (1998)
   Dim Y As Double, Ymin As Double, Ycrit As Double, Ymax As Double
   Dim TYmin As Double, T As Double
   Dim Iter As Integer
   Y = YY
   Select Case Label
      Case "S"
         Call SatfS(IDID, Qual, PS, xdPdT, TS, DL, DV, Y)
         Exit Sub
      Case "H"
         Ymin = 1.8766
         Ycrit = Hcrit
         Ymax = Hmax2P
         TYmin = 4.18
         JR = 9
         IDID = -125
      Case "U"
         Ymin = 1.8665
         Ycrit = Ucrit
         Ymax = Umax2P
         TYmin = 4.21
         JR = 11
         IDID = -127
      Case "G"
'     Following Gibbs energies are calculated from HEPAK v3.2
         Ymax = -8021.2
         Ymin = -1.8963
         Ycrit = Ymax
         JR = 12
         IDID = -128
      Case Else
         IDID = -299
         Exit Sub
   End Select
   If ((Y - Ycrit) * (Y - Ymin) <= 0#) Then
      Call RootAZ(T, Tmin, 2#, 4.5, Tcrit, E(8), 0#, _
                  E(11), E(9), 9, Y, Iter)
      If (Iter > 0) Then
         IDID = 2
         PS = SubRT(1)
         xdPdT = SubRT(6)
         TS = T
         DL = SubRT(3)
         DV = D2VfPT(SubRT(1), T)
         Qual = 0#
      End If
   ElseIf (((Y - Ycrit) * (Y - Ymax) <= 0#) And (JR <> 12)) Then
      Call RootAZ(T, TYmin, 4.5, 5#, Tcrit, E(8), 0#, _
                  E(11), E(9), 10, Y, Iter)
      If (Iter > 0) Then
         IDID = 2
         PS = SubRT(1)
         xdPdT = SubRT(6)
         TS = T
         DV = SubRT(3)
         DL = D2LfPT(SubRT(1), T)
         Qual = 1#
      End If
   End If
End Sub


Private Function SatLY(T As Double) As Double
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
'
' External function for SatFY; see it for definitions; v 1/5/93
   Call PsatfT(SubRT(1), SubRT(6), T)
   SubRT(3) = D2LfPT(SubRT(1), T)
   SubRT(2) = T
   Call SHAUG(SubRT(1), SubRT(2), SubRT(3), SubRT(4), SubRT(5), _
              SubRT(6), SubRT(7), SubRT(8), SubRT(9), SubRT(10), _
              SubRT(11), SubRT(12))
   SatLY = SubRT(JR)
End Function


Private Function SatVY(T As Double) As Double
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
'
' External function for SatFY; see it for definitions; v 1/5/93
   Call PsatfT(SubRT(1), SubRT(6), T)
   SubRT(3) = D2VfPT(SubRT(1), T)
   SubRT(2) = T
   Call SHAUG(SubRT(1), SubRT(2), SubRT(3), SubRT(4), SubRT(5), _
              SubRT(6), SubRT(7), SubRT(8), SubRT(9), SubRT(10), _
              SubRT(11), SubRT(12))
   SatVY = SubRT(JR)
End Function


Private Function T7658(T58 As Double) As Double
' Converted to Visual Basic by Jay Theilacker (1998)
'
' T76 as a function of T58; valid from 0.8 to Tc
' RMS fitting error = 66 microdegrees.
' constrained for 0 error at T58 = 2.172 and 5.190
' This function is not used by HEPAK.  It is listed here for reference only.
   Dim A(4) As Double, B(3) As Double, C(3) As Double
   Dim X As Double
   A(1) = -0.001891260993
   A(2) = 0.00850893784
   A(3) = -0.004836428596
   A(4) = 0.001076075419
   B(1) = 0.007530201452
   B(2) = -0.002681259464
   B(3) = 0.000978462867
   C(1) = 0.006115602751
   C(2) = 0.001146742367
   C(3) = -0.000276784046
   If (T58 <= 2.172) Then
      T7658 = T58 + A(1) + T58 * (A(2) + T58 * (A(3) + T58 * A(4)))
   ElseIf (T58 <= 4.082) Then
      X = 1# / (T58 - 1.7)
      T7658 = T58 + B(1) + X * X * (B(2) + X * B(3))
   ElseIf (T58 <= 4.619) Then
      T7658 = T58 + 0.00713
   ElseIf (T58 <= 5.19) Then
      X = 1# / (5.4 - T58)
      T7658 = T58 + C(1) + X * (C(2) + X * C(3))
   Else
      T7658 = -1#
   End If
End Function


Private Function RhoBB(P As Double, T As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
'  Beattie-Bridgeman equation used for approximate density
   Const A0 As Double = 0.0216
   Const ASM As Double = 0.05984
   Const B0 As Double = 0.014
   Const BSM As Double = 0#
   Const CSM As Double = 40#
   Const RBB As Double = 0.08206
   Dim PI As Double
   PI = RBB * T * 101325# / P
   RhoBB = GMWT / ((PI + B0 * (1# - BSM / PI)) * _
                  (1# - CSM / PI * T ^ (-3)) - _
             A0 * (1# - ASM / PI) / (RBB * T))
End Function


Private Sub Df2PT(IDID As Integer, D As Double, PP As Double, TT As Double)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
' This subroutine iterates for density D [Kg/m3], given P [Pa] and T [K].
' It is valid only in the compressed liquid near the lambda line, where
' it is more precise than subroutine DfPT
   Dim T As Double, P As Double, DelP As Double, DelD As Double
   Dim J As Integer
   T = TT
   P = PP
   D = D175K(T)
   IDID = 1
   For J = 1 To 5
      DelP = P - PressA(D, T)
      DelD = dPdDA(D, T)
      If (Abs(DelD) < 0.00000001) Then
         IDID = -207
         Exit Sub
      End If
      DelD = DelP / DelD
      D = D + DelD
      If (Abs(DelD) < 0.001) Then Exit Sub
   Next
   IDID = -207
End Sub


Private Sub DfPT(IDID As Integer, D As Double, X As Double, _
                 P As Double, T As Double)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
'  density given pressure and temperature [SI units].
   Const P20 As Double = 2000000#
   Dim Dmin As Double, Dmax As Double, DA As Double, DB As Double
   Dim Rho As Double, Tline As Double
   Dim JX As Integer
   SubRT(1) = T
   If (P > Pmax) Then
      IDID = -102
      Exit Sub
   ElseIf (P <= 0#) Then
      IDID = -101
      Exit Sub
   ElseIf (T < Tmin) Then
      IDID = -103
      Exit Sub
   End If
   Tline = 11# + P * 0.000001
   If (T > Tline) Then
' close to ideal gas
      If (T > Tmax) Then
         IDID = -104
         Exit Sub
      End If
      Rho = RhoBB(P, T)
      Dmax = DenMax(T)
      Dmin = 0.00000001
      DA = Min(Rho, Dmax)
      DB = DA * 0.9
' check saturation line
   ElseIf (T < Tcrit) Then
      If (P >= P20) Then
         Dmax = DMfT(T) + 0.3
         DA = 170#
         DB = 0.5 * (DA + Dmax)
         Dmin = 150#
      ElseIf (P < Psat(T)) Then
         DA = P / (2200# * T)
         DB = DGsat(T)
         Dmin = 0.00000001
         Dmax = DB * 1.03
      Else
         Dmin = 0.97 * DFsat(T)
         DB = D175K(P) + 1#
         DA = 0.5 * (DB + Dmin)
         Dmax = 180#
      End If
   ElseIf (P < Pcrit) Then
      DA = P / (HeGasConst * T)
      Dmax = 70#
      DB = Min(Dmax, (2# * DA))
      Dmin = 0.00000001
   Else
      Dmax = DenMax(T)
      Dmin = 5#
      If (P < P20) Then
         DA = 320# / T + P * 0.00002
         DB = 1.5 * DA
      Else
         DA = Min((80# + P * 0.0000004), (0.7 * Dmax))
         DB = 0.5 * (DA + Dmax)
      End If
   End If
   Call RootAZ(Rho, Dmin, DA, DB, Dmax, E(15), E(10), E(15), _
               E(10), 11, P, JX)
   If (JX <= 0) Then
      IDID = -201
      If (T < 14.05) Then IDID = -107
   Else
      D = Rho
      IDID = 1
      If (Rho > Dcrit) Then
         X = -1#
      Else
         X = 2#
      End If
   End If
End Sub


Private Function DPT(D As Double) As Double
' Function used by DfPT
   DPT = Press(D, SubRT(1))
End Function


Private Sub DfST(IDID As Integer, Rho As Double, X As Double, _
                 RF As Double, RG As Double, P As Double, _
                 dPdTs As Double, S As Double, T As Double)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
'  density given entropy and temperature [SI units].
   Dim Dmax As Double, Dmin As Double, DA As Double, DB As Double
   Dim Dum As Double, xPsat As Double, SF As Double, SG As Double
   Dim EYA As Double, JX As Integer
   If (S * (S - 64000#) >= 0#) Then
      IDID = -108
      Exit Sub
   End If
   If (T < Tmin) Then
      IDID = -103
      Exit Sub
   End If
   If (T > Tmax) Then
      IDID = -104
      Exit Sub
   End If
   X = 2#
   SubRT(1) = T
   Dmax = DenMax(T)
   If (T < Tcrit) Then
      Call PsatfT(xPsat, Dum, T)
      RG = D2VfPT(xPsat, T)
      SG = Entrop(RG, T)
      If (S >= SG) Then
         Dmax = 1.01 * RG
         DA = RG * Exp((SG - S) / HeGasConst)
         DB = Sqr(DA * Dmax)
         Dmin = 0.00000001
      Else
         RF = D2LfPT(xPsat, T)
         SF = Entrop(RF, T)
         If (S > SF) Then
            X = (S - SF) / (SG - SF)
            Rho = RF / (1# + X * (RF / RG - 1#))
            P = xPsat
            dPdTs = Dum
            IDID = 3
            Exit Sub
         Else
            Dmin = 0.99 * RF
            DA = 0.75 * Dmin + 0.25 * Dmax
            DB = 0.25 * Dmin + 0.75 * Dmax
         End If
      End If
   Else
      DA = 1.2 * DHI * Exp((SHI - S + CvHI * Log(T / THI)) / HeGasConst)
      DA = Min(DA, Dmax)
      Dmin = 0.00000001
      DB = Max((0.67 * DA), Dmin)
   End If
   EYA = CvHI * E(10)
   Call RootAZ(Rho, Dmin, DA, DB, Dmax, E(15), E(10), _
               EYA, E(10), 12, S, JX)
   If (JX <= 0) Then
      IDID = -108
   Else
      If (Rho > Dcrit) Then X = -1#
      IDID = 1
   End If
End Sub


Private Function DST(Y As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Function used by DfST
   DST = Entrop(Y, SubRT(1))
End Function


Private Sub DTfHS(IDID As Integer, D As Double, T As Double, _
                  X As Double, DL As Double, DV As Double, _
                  P As Double, xdPdT As Double, _
                  H As Double, S As Double)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Density, Temperature, and Quality as a function of Enthalpy and Entropy.
' In the two phase region (quality between 0 and 1), saturation densities,
' pressure, and dP/dT are also returned. SI units. Valid from 0.8 to 1500 K.
'
' Following are entropy and temperature when H(sat vap) is maximum
' Version August 29, 1992
   Const SHmax2 As Double = 8545#
   Const THmax2 As Double = 4.17
   
   Dim T1 As Double, D1 As Double, SL As Double, SV As Double
   Dim JX As Integer, I As Integer
' Following 2 elements of SubRT are not modified by the call to Root
   SubRT(4) = S
   SubRT(5) = H
   X = 2#
   If (S * (S - 64000#) >= 0#) Then
      IDID = -108
      Exit Sub
   End If
   If (H * (H - 8100000#) >= 0#) Then
      IDID = -109
      Exit Sub
   End If
   If (H < Hmax2P) Then
' If in the two-phase region, then we can find T such that H-TS = Gsat(T)
      Call RootAZ(T1, Tmin, 1.2, 3.6, Tcrit, E(8), 0#, _
                  E(9), E(10), 13, Zero, JX)
      If (JX > 0) Then
' Two phase
         T = T1
         Call PsatfT(P, xdPdT, T1)
         DL = D2LfPT(P, T1)
         DV = D2VfPT(P, T1)
         SL = Entrop(DL, T1)
         SV = Entrop(DV, T1)
         If (SV <> SL) Then
            X = (S - SL) / (SV - SL)
            D = 1# / ((1# / DL) * (1# - X) + (1# / DV) * X)
         Else
            X = 0.5
            D = Dcrit
         End If
         If (X < Zero) Then
            D1 = DL
         ElseIf (X > 1#) Then
            D1 = DV
         Else
            IDID = 3
            Exit Sub
         End If
      ElseIf (S <= Zero) Then
         IDID = -108
         Exit Sub
      ElseIf (S < Scrit) Then
' Compressed liquid
         T1 = Tmin + (S / Scrit) * (Tcrit - Tmin)
         D1 = 1.005 * DFsat(T1)
      ElseIf (S < SHmax2) Then
' Near saturated vapor, T > THMAX2
         T1 = (S - Scrit) / (SHmax2 - Scrit)
         T1 = Tcrit * (1# - T1) + THmax2 * T1
         D1 = 0.95 * DGsat(T1)
      ElseIf (H <= Zero) Then
         IDID = -109
         Exit Sub
      Else
' Vapor, T < Tcrit
' Latent heat at T = 0 is about 21000 J/kg
         T1 = Max(Tmin, ((H - 21000#) * THmax2 / Hmax2P))
         If (T1 >= Tcrit) Then T1 = Tcrit - 0.04
         D1 = DGsat(T1)
      End If
   Else
' use ideal gas approximation
      T1 = Max(THI + (H - HHI) / CpHI, Tcrit + 2#)
      T1 = Min(T1, Tmax)
      D1 = DHI * Exp((SHI - S - CvHI * Log(THI / T1)) / HeGasConst)
      D1 = Min(D1, DenMax(T1))
      D1 = Max(D1, 0.00001)
   End If
   Call DTILXY(D1, T1, I, "HS", H, S)
   If (I <= 0) Then
      IDID = -202
   Else
      D = D1
      T = T1
      DL = Zero
      DV = Zero
      If (D > Dcrit) Then X = -1#
      IDID = 1
   End If
End Sub


Private Function DTHS(T As Double) As Double
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' External function in DTfHS; Version Sept 9, 1989
   Dim Gsat As Double
   SubRT(2) = T
   Gsat = SubRT(5) - SubRT(2) * SubRT(4)
   Call PsatfT(SubRT(1), SubRT(6), SubRT(2))
   SubRT(3) = D2LfPT(SubRT(1), SubRT(2))
   Call SHAUG(SubRT(1), SubRT(2), SubRT(3), SubRT(4), SubRT(5), _
              SubRT(6), SubRT(7), SubRT(8), SubRT(9), SubRT(10), _
              SubRT(11), SubRT(12))
   DTHS = Gsat - SubRT(12)
End Function


Private Sub DTfPX(IDID As Integer, Rho As Double, TT As Double, _
                  Q As Double, DL As Double, DV As Double, _
                  dPdTs As Double, P As Double, _
                  X As Double, Label As String)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Output: Density (Rho), Temperature (TT), and Quality (Q)
'         If Q is between 0 and 1, the saturation pressure P, its
'         derivative dP/dT, and saturation densitites DL and DV
'         are also evaluated. SI units.
' Input:  Pressure (P) and parameter X, where X can be entropy (S),
'         enthalpy (H) or internal energy (U), depending on Label.
' Range:  valid over the full state equation, including 2-phase region
' Version 3.2; indices changed to match PROP array
   Dim FL(13) As Double, FV(13) As Double, Fder(11) As Double
   Dim FHIP1(8 To 11), FHIP2(8 To 11)
   Dim SpHt As Double, TA As Double, T As Double, TS As Double
   Dim D As Double, Dum As Double, DelTL As Double
   Dim X3 As Double, TM As Double, DM As Double
   Dim JF As Integer, I As Integer
   FHIP1(8) = 2900#
   FHIP1(9) = 712000#
   FHIP1(10) = 0#
   FHIP1(11) = 270000#
   FHIP2(8) = 17400#
   FHIP2(9) = 5450000#
   FHIP2(10) = 0#
   FHIP2(11) = 1510000#
   If (P <= 0#) Then
      IDID = -101
      Exit Sub
   ElseIf (P > Pmax) Then
      IDID = -102
      Exit Sub
   End If
   Select Case Label
      Case "PS"
         IDID = -108
         If (X * (X - 64000#) >= 0#) Then Exit Sub
         SpHt = CpHI
         JF = 8
         TA = THI * Exp((X - SHI + HeGasConst * Log(0.00001 * P)) / CpHI)
      Case "PH"
         IDID = -109
         If (X * (X - 8100000#) >= 0#) Then Exit Sub
         SpHt = CpHI
         JF = 9
         TA = THI + (X - HHI) / CpHI
      Case "PU"
         IDID = -110
         If (X * (X - 4730000#) >= 0#) Then Exit Sub
         SpHt = CvHI
         JF = 11
         TA = THI + (X - UHI) / CvHI
   End Select
   If (X < 1#) Then Exit Sub
   T = -1#
   If ((TA > 15#) And (P < 102000000#)) Then
' Check for approximate ideal gas
      D = RhoBB(P, TA)
      If ((P < 2000000#) Or (D < Dcrit)) Then T = TA
   End If
   If (T < 0#) Then
      If (P < Pcrit) Then
         If (P < 1.47) Then
'           Psat (0.8) = 1.4752 PA; low pressure vapor
            T = Max(Tmin, TA)
            D = P / (HeGasConst * T)
         Else
' Compare with saturated vapor
            TS = TsatfP(P)
            DV = D2VfPT(P, TS)
            If (JF = 8) Then
               FV(8) = Entrop(DV, TS)
               If (X > FV(8)) Then T = TS * Exp((X - FV(8)) / SpHt)
            Else
               FV(1) = P
               FV(2) = TS
               FV(3) = DV
               Call SHAUG(FV(1), FV(2), FV(3), FV(4), FV(5), _
                          FV(6), FV(7), FV(8), FV(9), FV(10), _
                          FV(11), FV(12))
               If (X > FV(JF)) Then T = TS + (X - FV(JF)) / SpHt
            End If
            If (T > 0) Then
' Single phase vapor
               D = DV * TS / T
            Else
' Not vapor; compare with saturated liquid
               DL = D2LfPT(P, TS)
               If (JF = 8) Then
                  FL(8) = Entrop(DL, TS)
               Else
                  FL(1) = P
                  FL(2) = TS
                  FL(3) = DL
                  Call SHAUG(FL(1), FL(2), FL(3), FL(4), FL(5), _
                             FL(6), FL(7), FL(8), FL(9), FL(10), _
                             FL(11), FL(12))
               End If
               If (X > FL(JF)) Then
' Two phase mixture
                  If (FL(JF) = FV(JF)) Then
                     Q = 0.5
                  Else
                     Q = Roun01((X - FL(JF)) / (FV(JF) - FL(JF)))
                  End If
                  Call PsatfT(Dum, dPdTs, TS)
                  Rho = DL / (1# + Q * (DL / DV - 1#))
                  TT = TS
                  IDID = 3
                  Exit Sub
               Else
' Compressed liquid
                  Call Deriv(Fder(1), Fder(2), Fder(3), Fder(4), Fder(5), _
                             Fder(6), Fder(7), Fder(8), Fder(9), Fder(10), _
                             Fder(11), DL, TS)
                  DelTL = (X - FL(JF)) / Fder(1)
                  If (JF = 8) Then DelTL = DelTL * T
                  T = TS + DelTL
                  If (T < Tmin) Then T = Tmin
                  D = Max(D3K(P), DL)
               End If
            End If
         End If
      ElseIf (P < 2540000#) Then
' Compressed liquid
         D = D3K(P)
         If (JF = 8) Then
            X3 = S3K(P)
            If (X > X3) Then
               T = 3# * Exp((X - X3) / SpHt)
               D = D * 3# / T
            Else
               T = Max(Tmin, (3# + (X - X3) / SpHt))
            End If
         Else
            If (JF = 9) Then
               X3 = H3K(P)
            Else
               X3 = U3K(P)
            End If
            T = Max(Tmin, (3# + (X - X3) / SpHt))
            If (X > X3) Then D = D * 3# / T
         End If
      ElseIf (P > 100000000#) Then
         If ((X - FHIP1(JF)) * (X - FHIP2(JF)) > 0#) Then Exit Sub
         D = 150# + 0.00000018 * P
         T = 180#
      Else
' Possible solid
         TM = TMfP(P)
         DM = DMfP(P)
         FL(1) = P
         FL(2) = TM
         FL(3) = DM
         Call SHAUG(FL(1), FL(2), FL(3), FL(4), FL(5), _
                    FL(6), FL(7), FL(8), FL(9), FL(10), _
                    FL(11), FL(12))
         If (X < FL(JF) * 0.999) Then
            IDID = -107
            Exit Sub
         End If
         T = TM + (X - FL(JF)) / SpHt
         If (JF = 8) Then T = TM * Exp((X - FL(JF)) / SpHt)
         D = DM * Sqr(TM / T)
      End If
   End If
   Call DTILXY(D, T, I, Label, P, X)
   If (I < 0) Then
      IDID = -206
   Else
      IDID = 1
      DL = 0#
      DV = 0#
      If (D > Dcrit) Then
         Q = -1#
      Else
         Q = 2#
      End If
      Rho = D
      TT = T
   End If
End Sub


Private Sub DTfPG(IDID As Integer, Dcalc As Double, Tcalc As Double, _
                  P As Double, G As Double)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Density and temperature as a function of pressure and Gibbs energy
' This subroutine valid only in the compressed liquid between 0.8 & 3 K
   Dim Prpty(13) As Double, T As Double, Pcalc As Double
   Dim G3K As Double, D As Double, Gcalc As Double
   Dim J As Integer
   If ((P < 1.46) Or (P > 3020000#)) Then
      IDID = -113
      Exit Sub
   End If
   Pcalc = P
   G3K = H3K(Pcalc) - 3# * S3K(Pcalc)
   D = D175K(Pcalc)
   T = 0.8
   Prpty(1) = Pcalc
   Prpty(2) = T
   Prpty(3) = D
   Call SHAUG(Prpty(1), Prpty(2), Prpty(3), Prpty(4), Prpty(5), _
              Prpty(6), Prpty(7), Prpty(8), Prpty(9), Prpty(10), _
              Prpty(11), Prpty(12))
   Gcalc = Prpty(12)
   If ((G - Gcalc) * (G - G3K) > 0#) Then
      IDID = -111
      Exit Sub
   End If
   T = 0.8 + (G - Gcalc) / (G3K - Gcalc) * 2.2
   Call DTILXY(D, T, J, "PG", P, G)
   If (J <= 0) Then
      IDID = -203
   Else
      IDID = 1
      Dcalc = D
      Tcalc = T
   End If
End Sub


Private Sub DTILXY(D As Double, T As Double, J As Integer, _
                   Label As String, X As Double, Y As Double)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
'-----OUTPUT
' D = Density     [kG/m3]
' T = Temperature [K]
' J = Iteration index, which should be a positive number
'     If J is negative, it signals convergence failure.
'-----INPUT
' D = an initial estimate for density
' T = an initial estimate for temperature
' LABEL = a character variable specifying the pair (X,Y)
' X = Thermodynamic variable identified by Label
' Y = Thermodynamic variable identified by Label
'----- Version August 29, 1992
   Dim CallDR As Boolean, I As Integer
   Dim F(11) As Double, Prpty(13) As Double
   Dim Fract As Double, Pcalc As Double, C2 As Double
   Dim DelP As Double, DelH As Double, DelD As Double
   Dim DelTT As Double, DelS As Double, DelU As Double
   Dim DelG As Double, CKT As Double, OMPKS As Double
   J = 0
   Fract = 0.3
   CallDR = True
   For I = 1 To 30
      Pcalc = Press(D, T)
      Prpty(1) = Pcalc
      Prpty(2) = T
      Prpty(3) = D
      Call SHAUG(Prpty(1), Prpty(2), Prpty(3), Prpty(4), Prpty(5), _
                 Prpty(6), Prpty(7), Prpty(8), Prpty(9), Prpty(10), _
                 Prpty(11), Prpty(12))
      If (CallDR) Then Call Deriv(F(1), F(2), F(3), F(4), F(5), F(6), _
                                  F(7), F(8), F(9), F(10), F(11), _
                                  D, T)
      C2 = F(7) * F(7)
      Select Case Label
         Case "PH"
'--------X is Pressure, Y is Enthalpy
            DelP = X - Pcalc
            DelH = Y - Prpty(9)
            If ((Abs(DelP) < E(9) * X + E(1)) And _
               (Abs(DelH) < E(9) * Y + E(1))) Then J = I
            DelD = ((1# + F(5)) * DelP - D * F(5) * DelH) / C2
            DelTT = F(8) * DelP + DelH / F(1)
         Case "PS"
'--------X is Pressure, Y is Entropy
            DelP = X - Pcalc
            DelS = Y - Prpty(8)
            If ((Abs(DelP / X) < E(9)) And _
               (Abs(T * DelS) < E(8) * F(1))) Then J = I
            DelD = (DelP - F(5) * D * T * DelS) / C2
            DelTT = (F(4) * DelP / D + T * DelS) / F(1)
         Case "PU"
'--------X is Pressure, Y is Energy
            DelP = X - Pcalc
            DelU = Y - Prpty(11)
            If ((Abs(DelP) < E(9) * X + E(1)) And _
               (Abs(DelU) < E(9) * Y + E(2))) Then J = I
            CKT = Pcalc * F(6)
            OMPKS = 1# - F(5) * CKT / F(3)
            DelD = (DelP - F(5) * D * DelU) / (C2 * OMPKS)
            DelTT = ((F(4) - CKT) * DelP / D + DelU) / (F(1) * OMPKS)
         Case "PG"
'--------X is Pressure, Y is Gibbs energy
            DelP = X - Pcalc
            DelG = Y - Prpty(12)
            If (((Abs(DelP / X) < E(10)) Or (Abs(DelP) < E(1))) _
               And (Abs(DelG / Prpty(8)) < E(6))) Then J = I
            DelD = (F(4) / (Prpty(8) * T)) * ((Prpty(8) / _
                   (F(5) * F(2)) - 1#) * DelP + D * DelG)
            DelTT = (DelP / D - DelG) / Prpty(8)
         Case "HS"
'--------X is Enthalpy, Y is Entropy
            DelH = X - Prpty(9)
            DelS = Y - Prpty(8)
            If ((Abs(DelH) < E(10) * Abs(X) + E(2)) And _
               (Abs(T * DelS) < E(10) * Abs(Y) + E(2))) Then J = I
            DelD = (DelH - (1# + F(5)) * T * DelS) * D / C2
            DelTT = (F(5) * DelH / C2 - D * F(8) * DelS) * T
      End Select
'-----iteration
      If ((Abs(DelD / D) + Abs(DelTT / T) < 0.0004) And (I > 3)) Then
         CallDR = False
         Fract = 0#
      End If
'-----Limit large steps
      If (Fract > 0#) Then
         If ((I > 7) And ((I / 4) * 4 = I)) Then
            DelD = DelD * 0.5
            DelTT = DelTT * 0.5
         End If
         DelD = Sign(Min(Abs(DelD), Abs(Fract * D)), DelD)
         DelTT = Sign(Min(Abs(DelTT), Abs(Fract * T)), DelTT)
      End If
      D = D + DelD
      T = T + DelTT
      If ((Abs(DelD) < E(12) * D) And (Abs(DelTT) < E(12) * T)) Then J = I
      If (J > 0) Then Exit Sub
   Next
   If (Fract < 0.01) Then
' Probable oscillation around the true answer
      D = D - 0.5 * DelD
      T = T - 0.5 * DelTT
      J = 61
      Exit Sub
   End If
   J = -200
End Sub


Private Sub TfDP(IDID As Integer, T As Double, X As Double, _
                 DL As Double, DV As Double, D As Double, P As Double)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Output: Temperature, quality, and saturation densities
' Input:  Density and pressure  (All SI units)
' GrunCv = Gruneisen parameter * Cv, and is reasonably constant.
' Used for integration w.r to P along isochores.
   Const GrunCv As Double = 2078#
   Const Pmx As Double = 101330000#
   Dim Tmx As Double, Tmn As Double, TA As Double, TB As Double
   Dim TT As Double, XX As Double, PS As Double, xdPdT As Double
   Dim TS As Double, PM As Double, JX As Integer
   X = -1#
   IDID = 1
   Tmx = 0#
   If (D <= 0#) Then
      IDID = -105
      Exit Sub
   ElseIf (P <= 0#) Then
      IDID = -101
      Exit Sub
   ElseIf (P > Pmx) Then
      IDID = -102
      Exit Sub
   End If
   If ((P < 1600000#) Or (D < Dcrit)) Then
      TA = Min((P / (HeGasConst * D)), 1400#)
      Tmn = Tmin
      If (TA > 30#) Then
         TB = 1.05 * TA
         Tmx = Tmax
      ElseIf ((TA > 6#) Or (P < 1.48)) Then
         TB = TA + 3#
         Tmx = Tmax
      End If
   End If
   If (Tmx = 0#) Then
      If (D < DLamLf) Then
         Call SatfD(JX, XX, PS, xdPdT, TS, DL, DV, D)
         If (JX < 0) Then
            IDID = -203
            Exit Sub
         End If
         If (P < PS) Then
            TT = TsatfP(P)
            DL = D2LfPT(P, TT)
            DV = D2VfPT(P, TT)
            If (DL = DV) Then
               X = 0.5
            Else
               X = (DL / D - 1#) / (DL / DV - 1#)
               X = Roun01(X)
            End If
            T = TT
            IDID = 3
            Exit Sub
         Else
            If (D < Dcrit) Then X = 2#
            Tmn = TS - 0.001
            Tmx = Tmax
            TA = Min((TS + 0.75 * (P - PS) / (D * GrunCv)), 1400#)
            TB = Min((TA * 1.25), Tmx)
         End If
      Else
' 179.84 = density at upper lambda point, where T = 1.7673
         If (D < 179.84) Then
            Tmn = TLfD(D) + 0.2
            JX = -1
         Else
            Tmn = TMfD(D)
            JX = 0
         End If
         PM = Press(D, Tmn)
         TA = 0.75 * (P - PM) / (D * GrunCv)
         If (TA < 0#) Then
            If (JX < 0) Then
               IDID = -180
            Else
               IDID = -107
            End If
            Exit Sub
         End If
         TA = TA + Tmn
         TB = TA * 1.25
         Tmx = 400#
      End If
   End If
   SubRT(1) = D
   Call RootAZ(TT, Tmn, TA, TB, Tmx, E(8), E(11), _
               E(2), E(11), 14, P, JX)
   If (JX < 0) Then
      IDID = -106
   Else
      T = TT
   End If
End Sub


Private Function TDP(Y As Double) As Double
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Function used by TfDP
   TDP = Press(SubRT(1), Y)
End Function


Private Sub TfDY(IDID As Integer, T As Double, Q As Double, _
                 PS As Double, xdPdT As Double, DL As Double, _
                 DV As Double, DD As Double, _
                 YY As Double, Label As String)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Output: Temperature (T) and Quality (Q)
'         If Q is between 0 and 1, the saturation pressure P, its
'         derivative dP/dT, and saturation densitites DL and DV
'         are also evaluated.  SI units.
' Input:  Density D and parameter X, where Y can be entropy (S),
'         enthalpy (H) or internal energy (U), depending on Label.
' Range:  valid over the full state equation, including 2-phase region
' Version Aug. 5, 1992; for HEPAK v3.2
   Dim Y As Double, D As Double, JX As Integer
   Dim TA As Double, TB As Double, Tc As Double
   Y = YY
   D = DD
   If (D * (D - 599#) >= 0#) Then
      IDID = -105
      Exit Sub
   End If
   If (D > 171.7) Then
      TA = TMfD(D)
   Else
      TA = Tmin
   End If
   Select Case Label
      Case "S"
         IDID = -108
         If (Y * (Y - 64000#) >= 0#) Then Exit Sub
         TB = Exp(((0.001 * Y - 10.037) / 2.07723 + Log(D)) / 1.5)
         TB = Max(TA, TB)
         JR = 8
      Case "H"
         IDID = -109
         If (Y * (Y - 8100000#) >= 0#) Then Exit Sub
         TB = Max((Y / CpHI), TA)
         JR = 9
      Case "U"
         IDID = -110
         If (Y * (Y - 4730000#) >= 0#) Then Exit Sub
         TB = Max((Y / CvHI), TA)
         JR = 11
      Case Else
         IDID = -299
         Exit Sub
   End Select
   TB = Max(TB, 1#)
   If (TB > 1200#) Then
      TB = 1200#
      Tc = 1400#
   Else
      Tc = 2# * TB
      Tc = Min(Tc, (TB + 50#))
   End If
   SubRT(3) = D
   Call RootAZ(T, TA, TB, Tc, Tmax, E(9), E(11), _
               E(3), E(10), 15, Y, JX)
   If (JX < 0) Then Exit Sub
   IDID = CInt(SubRT(7))
   If (IDID = 3) Then
      Q = SubRT(8)
      PS = SubRT(1)
      DL = SubRT(4)
      DV = SubRT(5)
      xdPdT = SubRT(6)
   ElseIf (D < Dcrit) Then
      Q = 2#
   Else
      Q = -1#
   End If
End Sub


Private Function YJfDT(T As Double) As Single
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' External function in TfDY.
' Output:
'     YJfDT = entropy, enthlpy, or energy [SI units]
'     (if IDID = 3) saturation properties P, dPdT, Q, DL, DV, via SubRT
' Input:  T = temperature (called from RootAZ)
'         D = density, from TfDY via SubRT
' Version Sept 13, 1992
   Dim F(13) As Double, D As Double, Q As Double, PS As Double
   Dim xdPdT As Double, DL As Double, DV As Double
   Dim IDID As Integer
   D = SubRT(3)
   Call FDT(IDID, Q, PS, xdPdT, DL, DV, D, T)
   SubRT(7) = CDbl(IDID)
   If (IDID < 0) Then IDID = 1
   If (IDID = 1) Then
      DL = D
      PS = Press(D, T)
   End If
   If (JR = 8) Then
      YJfDT = Entrop(DL, T)
   Else
      F(1) = PS
      F(2) = T
      F(3) = DL
      Call SHAUG(F(1), F(2), F(3), F(4), F(5), F(6), _
                 F(7), F(8), F(9), F(10), F(11), F(12))
      YJfDT = F(JR)
   End If
   If (IDID = 3) Then
      If (JR = 8) Then
         YJfDT = Q * Entrop(DV, T) + (1# - Q) * YJfDT
      Else
         F(3) = DV
         Call SHAUG(F(1), F(2), F(3), F(4), F(5), F(6), _
                    F(7), F(8), F(9), F(10), F(11), F(12))
         YJfDT = Q * F(JR) + (1# - Q) * YJfDT
      End If
      SubRT(8) = Q
      SubRT(1) = PS
      SubRT(4) = DL
      SubRT(5) = DV
      SubRT(6) = xdPdT
   End If
End Function


Private Sub DfTG(IDID As Integer, D As Double, T As Double, G As Double)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
' Density as a function of Temperature and Gibbs Energy (all SI units).
' This function operates only in the compressed liquid between 0.8 and
' 3 K (applicable to the fountain effect in superfluid).
   Dim Dmax As Double, Dmin As Double, D1 As Double, D2 As Double
   Dim JX As Integer
   If (T < 0.8) Then
      IDID = -103
      Exit Sub
   ElseIf (T > 3#) Then
      IDID = -112
      Exit Sub
   End If
   SubRT(2) = T
   Dmin = 0.995 * DFsat(T)
   Dmax = DMfT(T)
   D1 = 0.9 * Dmin + 0.1 * Dmax
   D2 = 0.5 * (Dmin + Dmax)
   Call RootAZ(D, Dmin, D1, D2, Dmax, E(7), E(11), _
               E(6), E(11), 16, G, JX)
   If (JX < 0) Then
      IDID = -111
   Else
      IDID = 1
   End If
End Sub


Private Function DTG(D As Double) As Double
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
'
' Gibbs energy as a function of density and (from SubRT) T.
' This is an external function to subroutine DfTG.
   SubRT(3) = D
   SubRT(1) = PressA(SubRT(3), SubRT(2))
   Call SHAUG(SubRT(1), SubRT(2), SubRT(3), SubRT(4), SubRT(5), _
              SubRT(6), SubRT(7), SubRT(8), SubRT(9), SubRT(10), _
              SubRT(11), SubRT(12))
   DTG = SubRT(12)
End Function


Private Function DielHe(Rho As Double) As Double
' Converted to Visual Basic by Jay Theilacker (1998)
'    (Dielectric constant - 1.0) as a function of density [Kg/m3].
'-----Version 3.23: Output changed to (Dielectric constant - 1.0)
'                   instead of just the dielectric constant.
'    Equation from D. Friend, NIST, 12/20/90.  Includes values from Gugan
'    and Michel,  Metrologia 1980, 1984 with 1986 CODATA constants;
'    temperature dependence of B is eliminated.
'
   Dim P As Double
   P = (0.123493 - 0.00000586 * Rho) * 1.046516792
   DielHe = (2# * P * Rho + 1000#) / (1000# - P * Rho) - 1#
End Function


Private Function Refr(Rho As Double, W As Double) As Double
' Converted to Visual Basic by Jay Theilacker (1998)
'     (Refractive index -1.0) as a function of density and wavelength [um].
'-----Version 3.23: Output changed to (refr. index - 1.0)
'                   instead of just the refractive index.
'     Uses new dielectric constant correlation for infinite wavelength.
'     Uses Sellmeier fit from Peck (App. Opt. 22, 2906, 1983) for dispersion
'     between 0.09 and 2 micrometer [= um] wavelength.
'     For larger wavelengths, linearly interpolates between 2 um value and
'     dielectric constant value.  The dispersion at 298.15 k, 0.101325 mpa
'     is assumed to be unchanged at other state points.
'     Equation from D. Friend (NIST), Jan 3 1991.
   Const RInf As Double = 0.000031704
   Const R2min As Double = 0.00003464
   Dim W2I As Double
   Dim RInInf As Double, RDisp As Double
   If (W <= 0#) Then
      W2I = 1#
   Else
      W2I = 1# / (W * W)
   End If
   RInInf = Sqr(1# + DielHe(Rho))
   If (W2I < 0.25) Then
      RDisp = W2I * (R2min - RInf) / 0.25
   Else
      RDisp = 0.0000073123 + 0.0092797 / (339.82 - W2I) - RInf
   End If
   Refr = RInInf + RDisp - 1#
End Function


Private Function STen(T As Double) As Double
' Converted to Visual Basic by Jay Theilacker (1998)
'  Surface tension [N/m] as a function of temperature [K].
'  Equation from D. Friend, NIST, Jan 8, 1991.  Fit of data from
'  Iino, Suzuki, Ikushima, JLTP 61, 155 (1985).
   Const A1 As Double = 0.000389057
   Const A2 As Double = 0.00052141
   Const A3 As Double = -0.000579737
   Dim TSTR As Double
   If (T >= Tcrit) Then
      STen = 0#
   Else
      TSTR = (Tcrit - T) / Tcrit
      STen = TSTR * (A1 + TSTR * (A2 + TSTR * A3))
   End If
End Function


Private Function TCon(Rho As Double, T As Double) As Double
' (C) Copyright (1990), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
' This function written by B.A. Hands.
'
'  thermal conductivity as a function of density and temperature
   Dim D As Double, D2 As Double, D3 As Double, DL As Double
   Dim T1 As Double, T2 As Double
   Dim Kt As Double, KT1 As Double, Kcrit As Double
   Dim DelD As Double, DelTc As Double
   Dim X1 As Double, XX2B As Double, XX2BE As Double
   Dim R2 As Double, W As Double, H As Double, A As Double
   Dim DHDX As Double, D2KT As Double
   Dim PDT As Double, TT As Double, K0 As Double
   Dim I As Integer
   Dim F(12) As Double, C(5) As Double
   F(1) = 3.726229668
   F(2) = 0.000186297053
   F(3) = -7.275964435E-07
   F(4) = -0.0001427549651
   F(5) = 0.00003290833592
   F(6) = -5.213335363E-08
   F(7) = 4.492659933E-08
   F(8) = -5.924416513E-09
   F(9) = 0.000007087321137
   F(10) = -0.000006013335678
   F(11) = 8.067145814E-07
   F(12) = 3.995125013E-07
   C(1) = 0.7034007057
   C(2) = 3.739232544
   C(3) = -26.20316969
   C(4) = 59.82252246
   C(5) = -49.26397634
   Const X0 As Double = 0.392
   Const E1 As Double = 2.8461
   Const E2 As Double = 0.27156
   Const Beta As Double = 0.3554
   Const Gamma As Double = 1.1743
   Const Delta As Double = 4.304
   Const Dc As Double = 69.58
   Const Tc As Double = 5.18992
   Const Pc As Double = 227460#
   Const Cons As Double = -5.882788298
   Const Con As Double = 3.4685233E-17
   D = Rho
   If (T < 3.5) Then
      If (D < Dc) Then
         TCon = ConLPT(Rho, T)
      Else
         TCon = 0#
      End If
      Exit Function
   End If
   T1 = T ^ (1# / 3#)
   T2 = T1 * T1
   D2 = Rho * Rho
   D3 = Rho * D2
   If (Rho < 0.00000001) Then
      DL = 0#
   Else
      DL = D2 * Log(Rho / 69.64)
   End If
' Calculate critical enhancement
   If (T > 12# Or T < 3.5) Then
      Kcrit = 0
   Else
' Calculate compressibility
      Kt = 1# / dPdD(Rho, T) / Rho
      DelD = Abs((Rho - Dc) / Dc)
      DelTc = Abs((T - Tc) / Tc)
      R2 = (DelTc / 0.2) ^ 2 + (DelD / 0.25) ^ 2
      If (R2 < 1) Then
         W = DelTc / (DelD ^ (1# / Beta))
         X1 = (W + X0) / X0
         XX2B = X1 ^ (2# * Beta)
         XX2BE = (1# + E2 * XX2B) ^ ((Gamma - 1#) / 2# / Beta)
         H = E1 * X1 * XX2BE
         DHDX = E1 * XX2BE / X0 + _
                E1 * E2 / X0 * XX2B * XX2BE / (1# + E2 * XX2B) * _
                (Gamma - 1#)
         D2KT = (Delta * H - W * DHDX / Beta) * (DelD ^ (Delta - 1#))
         KT1 = Dc * Dc / D2 / D2KT / Pc
         Kt = R2 * Kt + (1# - R2) * KT1
      End If
' Calculate excess conductivity
      PDT = dPdT(Rho, T)
      Kcrit = T * T * Sqr(Kt) / Rho / Viscos(Rho, T) * PDT * PDT * _
              Exp(-18.66 * DelTc ^ 2 - 4.25 * DelD ^ 4)
      Kcrit = Kcrit * Con
   End If
   A = 0
   TT = T
   For I = 2 To 5
      A = C(I) / TT + A
      TT = TT * T
   Next
   K0 = T ^ (C(1)) * Exp(A + Cons)
   DL = D2 * Log(Rho / 68#)
   TCon = K0 + F(1) * Kcrit + Rho * (F(2) + F(3) * T + F(4) * T1 + _
          F(5) * T2) + D3 * (F(6) + F(7) * T1 + F(8) * T2) + _
          DL * (F(9) + F(10) * T1 + F(11) * T2 + F(12) / T)
End Function


Private Function ConLPT(D As Double, T As Double) As Single
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
' Conductivity of vapor, T < 3.5 K; extrapolation of HEPAK v3.0 data.
' The accuracy of this extrapolation is unknown.  Existing data from
' 3.5 to 20 K suggests anomalous behavior in the low T limit.
' The assigned lower limit of 1.2 K is arbitrary.
'
' Units: ConLPT [W/m-K], Density [Kg/m3], Temperature [K]
' V. Arp; August, 1992
   Const A As Double = 0.001916377113
   Const B As Double = 0.00004654575892
   If (T > 1.1999999) Then
      ConLPT = A * T + B * D
   Else
      ConLPT = 0#
   End If
End Function


Private Function Viscos(DKgM3 As Double, TK As Double) As Double
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
' Output: Viscosity [SI units = Pa-s]
' Input:  Density [Kg/m3] (single phase region only!), and Temperature [K]
'-----Version 3.2; 7/2/92
'
   Dim D As Double, T As Double, VX As Double, W As Double
   Const TL As Double = 3.5
   Const TH As Double = 3.75
   Const DT As Double = TH - TL
   D = DKgM3
   T = TK
   If (T > TL) Then
      Viscos = ViscHI(D, T)
      If (T >= TH) Then Exit Function
   End If
   If (D > Dcrit) Then
      VX = VisLam(D, T)
   Else
      VX = VisLPT(D, T)
   End If
   If (T <= TL) Then
      Viscos = VX
      Exit Function
   End If
   W = (T - TL) / DT
   Viscos = W * Viscos + (1# - W) * VX
End Function


Private Function ViscHI(DKgM3 As Double, T As Double) As Double
' Converted to Visual Basic by Jay Theilacker (1998)
' Output: Viscosity [SI units = Pa-s]
' Input:  Density [Kg/m3], and Temperature [K] (valid only for T>3.5K)
'
' This is consistent with NIST's viscosity function, January, 1991,
' by D.G.Friend, V. Arp and R.D. McCarty.
'
   Dim TL As Double
   Dim ViscL As Double, Visc0 As Double, ViscH As Double
   Dim R As Double, B As Double, C As Double, D As Double
   Const T1 As Double = 100#
   Const T2 As Double = 110#
   Const DelTT As Double = T2 - T1
   Const TmaxV As Double = 300
   If (T <= TmaxV) Then
      TL = Log(T)
   Else
      TL = Log(Tmax)
   End If
   ViscL = Exp(-0.135311743 / TL + 1.00347841 + TL * (1.20654649 + _
           TL * (-0.149564551 + TL * 0.0125208416)))
   Visc0 = ViscL
   If (T > T1) Then
      ViscH = 196# * T ^ 0.71938 * _
              Exp((12.451 - 295.67 / T) / T - 4.1249)
      If (T < T2) Then
         Visc0 = ViscL + (ViscH - ViscL) * (T - T1) / DelTT
      Else
         Visc0 = ViscH
      End If
   End If
   R = 0.001 * DKgM3
   B = -47.5295259 / TL + 87.6799309 + TL * (-42.0741589 + _
       TL * (8.33128289 - 0.589252385 * TL))
   C = 547.309267 / TL - 904.870586 + TL * (431.404928 + _
       TL * (-81.4504854 + TL * 5.37008433))
   D = -1684.39324 / TL + 3331.0863 + TL * (-1632.19172 + _
       TL * (308.804413 - 20.2936367 * TL))
   ViscHI = 0.0000001 * (Visc0 + _
            ViscL * (Exp(R * (B + R * (C + R * D))) - 1#))
End Function


Private Function VisLam(DD As Double, TT As Double) As Single
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
' Liquid helium viscosity, 1.2 to 3.8 K
'     Fitted to HMeyer's data along the saturation line.
'     Extended to higher densities by using Hoffman's (Brewer and Edwards)
'     data, shifted by 6.2 pct; joined smoothly to Steward's data
'     at T = 3.6 to 3.8 K.
'     The calling program must make certain that the input density
'     is within the liquid range (nominally 132 to 180+ Kg/m3).
' Units: VisLam [Pa-s], DD [Kg/m3], TT [K]
   Dim Dens As Double, TL As Double, T As Double
   Dim X As Double, D As Double, V As Double, E As Double
   Dim J As Integer, Y1 As Double, Y2 As Double
   Dim A(3) As Double, E0(2) As Double
   Dim B(8) As Double, C(8) As Double
   E0(1) = 1.8
   E0(2) = -4.79
   A(1) = 2.505885162
   A(2) = 0.5230553382
   A(3) = 0.5607799718
   B(1) = -112.7424846
   B(2) = 209.5894826
   B(3) = -128.6503418
   B(4) = -79.51538104
   B(5) = 201.5521019
   B(6) = -123.1199069
   B(7) = -64.60724357
   B(8) = 54.24829902
   C(1) = 877.2148954
   C(2) = -2515.234338
   C(3) = 2679.676294
   C(4) = 346.9587682
   C(5) = -1509.946785
   C(6) = 2056.348276
   C(7) = 288.6724375
   C(8) = -383.1832082
   Dens = DD
   TL = TLfD(Dens)
   T = TT
   If (T < 1.1999999) Then
      VisLam = 0#
      Exit Function
   End If
   X = Abs(1# - T / TL)
   D = 10# * (Dens / DLamLf - 1#)
   If (T <= TL) Then
      J = 2
      Y1 = C(4) + X * (C(5) + X * C(6)) + D * X * (C(7) + X * C(8))
      V = X * X * (C(1) + X * (C(2) + X * C(3)) + D * Y1)
   Else
      J = 1
      Y2 = B(4) + X * (B(5) + X * B(6)) + D * (B(7) + X * B(8))
      V = X * X * (B(1) + X * (B(2) + X * B(3)) + D * Y2)
   End If
   E = 10# * (1# + E0(J) * ((X + 0.00000001) ^ 0.84))
   V = V + E * (A(1) + D * (A(2) + D * A(3)))
   VisLam = V * 0.0000001
End Function


Private Function VisLPT(D As Double, T As Double) As Single
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
' Viscosity of vapor, T < 3.5 K; extrapolation of HEPAK v3.0 data.
' The accuracy of this extrapolation is unknown.  Existing data from
' 3.5 to 20 K suggests anomalous behavior in the low T limit.
' The assigned lower limit of 1.2 K is arbitrary.
'
' Units: VisLPt [Pa-s], Density [Kg/m3], Temperature [K]
   Const A As Double = 0.2539133848
   Const B As Double = 0.01024277783
   If (T > 1.1999999) Then
      VisLPT = (A * T + B * D) * 0.000001
   Else
      VisLPT = 0#
   End If
End Function


Private Sub RootAZ(R As Double, A As Double, B As Double, _
                   C As Double, Z As Double, _
                   EXA As Double, EXR As Double, _
                   EYA As Double, EYR As Double, _
                   I As Integer, F0 As Double, JX As Integer)
' (C) Copyright (1992), Cryodata Inc.; all rights reserved.
' Converted to Visual Basic by Jay Theilacker (1998)
'
'-----DESCRIPTION
' This subroutine finds the root, R, such that F(R) = F0,
' where F is an external function, and F0 is a specified constant.
' The function is accessed through the function RootFunc.
'
'-----EXTERNAL VARIABLES
' This version of RootAZ assumes that all external floating-point
' variables are Double
' OUTPUT VARIABLES
'   R           = the desired root,     if JX is positive;
'               = an approximate root,  if JX is zero;
'               = an invalid root,      if JX is negative.
'   JX          = number of calls to F, if RootAZ was successful.
'                 If JX=0, input error tolerances may be invalid.
'                 If JX = Jmax+4 = 54, iteration has been terminated without
'                 satisfying any convergence criterion.
'               = -1 if input values of A, B, C, and Z are invalid
'               = -2 if input values of EXA, EXR, EYA, or EYR are invalid.
'               = -3 if extrapolation beyond [B, C] is indeterminant.
'               = -4 if extrapolation beyond A or Z was attempted.
'               = -5 if extrapolation of a non-monotonic function
'                     was attempted.
' INPUT VARIABLES
'   A and Z     = limiting values of the independent variable which must
'                 bracket the root R.  A must not equal Z.  RootAZ will
'                 not extrapolate beyond these limits.
'   B and C     = Initial guesses for the root R.  B must not equal C,
'                 and both must be within or at the limits of the
'                 range [A, Z].
'   EXA and EXR = absolute and relative tolerances on the independent
'                 variable.  R is accepted as the root if it is
'                 bracketed within the range EXA + EXR*R
'                 EXA must be > zero.  EXR must be >= zero.
'   EYA and EYR = absolute and relative tolerances on the dependent
'                 variable.  R is accepted as the root if the absolute
'                 value of F(R) - F0 is less than EYA + EYR*F0
'                 EYA must be > zero.  EYR must be >= zero.
'   F0          = the specified value of the function
'
' Note:  RootAZ does not modify any of these input variables
'
'   I           = Reference number for the external function. See RootFunc
'                 for the list of function names.  It can be non-monotonic
'                 if B and C bracket R, but must be monotonic if
'                 extrapolation is required..
'-----INTERNAL VARIABLES
' All internal floating point variables are Double
' X(1)...X(4)     are successive approximations to R.
' Y(1)...Y(4)     are successive values of F(X(i)) - F0.
' Esys            specifies the machine precision of the external variables.
'                 It must be larger than the smallest number EPS such that
'                 1 + EPS is not equal to 1 in double precision arithmetic.
' Jmax            should be chosen such that 2**Jmax > (1/Esys)**2
'                 in order to assure convergence for any F.
' Other variables are used in intermediate steps.
'
'-----VERSION Sept. 13, 1992;  V. Arp
   Const Esys As Double = 0.00000001
   Const OME As Double = 1# - Esys
   Const Jmax As Integer = 50
   Dim Limit As Boolean
   Dim X(4) As Double, Y(4) As Double
   Dim EY As Double, R21 As Double, R32 As Double, W As Double
   Dim Step As Double, Xlim As Double, XI As Double, XE As Double
   Dim XJD As Double, XJS As Double, XJO As Double, YJD As Double
   Dim JD As Integer, JF As Integer, JG As Integer
   Dim JJ As Integer, JO As Integer, JS As Integer
'
' Check for valid parameters
'
   If ((A = Z) Or (B = C) Or ((B - A) * (B - Z) > 0#) Or _
      ((C - A) * (C - Z) > 0#)) Then
      JX = -1
      Exit Sub
   End If
   If ((EXA <= 0#) Or (EXR < 0#) Or (EYA <= 0#) _
      Or (EYR < 0#)) Then
      JX = -2
      Exit Sub
   End If
   Limit = False
   EY = EYA + (Esys + EYR) * Abs(F0)
'
'  First guess
'
   Y(1) = RootFunc(I, B) - F0
   If (Abs(Y(1)) < EY) Then
      JX = 1
      R = B
      Exit Sub
   End If
'
'  Second guess
'
   Y(2) = RootFunc(I, C) - F0
   If (Abs(Y(2)) < EY) Then
      JX = 2
      R = C
      Exit Sub
   End If
   R21 = Y(2) / Y(1)
'
'  Order such that Y(2) is smaller than Y(1).  This is required for
'  extrapolation, where monotonic F(x) is assumed.
'  If extrapolation is not required, F(x) can be non-monotonic.
'
   If (Abs(R21) > 1# + Esys) Then
      R32 = Y(2)
      Y(2) = Y(1)
      Y(1) = R32
      X(1) = C
      X(2) = B
   ElseIf (R21 < OME) Then
      X(1) = B
      X(2) = C
   Else
'  Failure: Constant function F; indeterminant extrapolation.
      JX = -3
      Exit Sub
   End If
'
'  Third try: linear interpolation or extrapolation
'
   Step = X(2) - X(1)
   X(3) = Y(2) * Step / (Y(1) - Y(2)) + X(2)
   If (R21 > 0#) Then
'
'  Extrapolation; check limits.
'  First, find which limit applies.
'
      If (((A - X(2)) / Step) >= 0#) Then
         Xlim = A
      Else
         Xlim = Z
      End If
      If (Abs(X(2) - Xlim) < _
         0.5 * (EXA + (EXR + Esys) * Abs(X(2)) + _
         (EXR + Esys) * Abs(X(2)))) Then
'  Failure: A=B or A=C or Z=B or Z=C; extrapolation is not permitted.
         JX = -4
         Exit Sub
      End If
'  Do not extrapolate beyond XLIM.
      If ((X(3) - Xlim) * (X(3) - X(1)) >= 0#) Then
         X(3) = Xlim
         Limit = True
      End If
   ElseIf (Abs(Step) <= _
          (EXA + (EXR + Esys) * Abs(X(1)) + _
          (EXR + Esys) * Abs(X(2)))) Then
'  B and C straddle the root.
'  Return if their separation is is less than the minimum acceptable
'  separation specified by EXA, EXR, and ESYS.   User should check for
'  input error.
      JX = 0
      R = X(3)
      Exit Sub
   End If
   Y(3) = RootFunc(I, X(3)) - F0
   If (Abs(Y(3)) < EY) Then
      R = X(3)
      JX = 3
      Exit Sub
   End If
   JF = 0
   JG = 0
'
'-----Do loop
'
   For JJ = 0 To Jmax
      R21 = Y(2) / Y(1)
      R32 = Y(3) / Y(2)
      If ((R32 > 0#) And (R21 > 0#)) Then
'
'  Root not bracketed; extrapolation still required.
'
         If (Limit) Then
'     Failure; further extrpolation is prohibited.
            JX = -4
            Exit Sub
         End If
         If (R32 >= OME) Then
'     Extrapolation failure: constant or non-monotonic function.
            JX = -5
            Exit Sub
         End If
         XI = Y(3) * (X(3) - X(2)) / (Y(2) - Y(3))
         If (R21 > R32) Then
'     Use inverse quadratic extrapolation if abs(slope) is decreasing.
            XE = Y(3) * (X(3) - X(1)) / (Y(1) - Y(3))
            W = Y(1) / (Y(1) - Y(2))
            XI = W * XI + (1# - W) * XE
         End If
         X(4) = XI + X(3)
         If (R32 > 0.1) Then
'     If extrapolation convergence is slow, give an artificial boost.
            If (JJ >= 2) Then
               X(4) = X(4) + CDbl(JJ - 1) * Step
            Else
               X(4) = X(4) + CDbl(JJ + 1) * (X(3) - X(2))
            End If
         End If
'     But do not extrapolate beyond the limit.
         If ((X(4) - Xlim) * (X(4) - X(1)) > 0#) Then
            X(4) = Xlim
            Limit = True
         End If
         JF = 0
         JS = 2
      Else
'
'  Root has been bracketed.
'    Define the geometry.
         If (R32 < 0#) Then
            JO = 1
            JF = 0
            If (R21 < 0#) Then
               JS = 2
               JD = 3
               JG = JG + 1
            Else
               JS = 3
               JD = 2
               JG = 0
            End If
         Else
            JS = 1
            JD = 3
            JO = 2
            JF = JF + 1
            JG = 0
         End If
         XJD = X(JD)
         XJS = X(JS)
         XJO = X(JO)
         YJD = Y(JD)
'    Exit if the bracket width is less than EX
         If (Abs(XJS - XJD) <= _
            (EXA + (EXR + Esys) * Abs(XJS) + _
            (EXR + Esys) * Abs(XJD))) Then
            JX = JJ + 3
            R = X(3)
            Exit Sub
         End If
         If ((Abs(R32) >= OME) Or _
            ((JG >= 2) And (R32 < -0.7))) Then
'    Not converging: bisect.
            X(4) = 0.5 * (XJS + XJD)
            JF = 0
         Else
'
'    The following approximates quadratic interpolation when converging
'    nicely.  Otherwise it can accelerate the search when large changes
'    in slope and/or slow movement toward a distant root are found.
'
            XI = YJD * (XJD - XJS) / (Y(JS) - YJD) + XJD
            If (Abs(Y(JO) - YJD) < EY) Then
               XE = XJS
            Else
               XE = YJD * (XJD - XJO) / (Y(JO) - YJD) + XJD
               If ((XE - XJS) * (XE - XJD) > 0#) Then XE = XJS
            End If
            W = (XJD - XJO) / (XJS - XJO)
            If ((JF >= 2) And (W < 0.3) And _
               (R32 > 0.25)) Then
               X(4) = X(3) + (X(1) - X(3)) * Sqr(W)
            Else
               X(4) = W * XI + (1# - W) * XE
            End If
         End If
      End If
'  X(4), Y(4) are the new estimate.
      Y(4) = RootFunc(I, X(4)) - F0
      If (Abs(Y(4)) < EY) Then
         JX = JJ + 4
         R = X(4)
         Exit Sub
      End If
'  Update indices.
      If (JS >= 2) Then
         X(1) = X(2)
         Y(1) = Y(2)
      End If
      X(2) = X(3)
      X(3) = X(4)
      Y(2) = Y(3)
      Y(3) = Y(4)
   Next
'  This statement should never be reached
   R = X(4)
   JX = Jmax + 4
End Sub