Commit 08f2f6bd authored by Zolnierczuk, Piotr's avatar Zolnierczuk, Piotr
Browse files

added water density to water_model

parent d4d72fde

pysen/physics/water_model.py

+24 −9
Original line number Diff line number Diff line
@@ -11,6 +11,8 @@ from numpy import ( asarray, arange, ones_like, newaxis, sqrt, exp, log, 
                    interp, polyval )

from numpy import sum as npsum



T0C = 273.15 # 0C in Kelvin



# Model parameters, see Table

ns   = 1000.0       # unit of time is ps in this model

kB   = 0.08617      # meV/K
@@ -327,15 +329,15 @@ RA = 1.9872036 # cal/mol/K 
# TABLE I in [3], sorted in increasing temp

Teixeira_Data = asarray([

    #T[K]	     t0[ps]	L[A]

    273.15-20.0, 22.7, 	2.39,

    273.15-17.0, 10.8, 	1.80,

    273.15-15.0, 8.90, 	1.73,

    273.15-12.0, 7.63, 	1.70,

    273.15-10.0, 6.47, 	1.65,

    273.15 -5.0, 4.66,	1.54,

    273.15 +5.0, 2.33,	1.32,

    273.15+12.0, 1.66,	1.25,

    273.15+20.0, 1.25,	1.29,

    T0C-20.0, 22.7, 	2.39,

    T0C-17.0, 10.8, 	1.80,

    T0C-15.0, 8.90, 	1.73,

    T0C-12.0, 7.63, 	1.70,

    T0C-10.0, 6.47, 	1.65,

    T0C -5.0, 4.66,	1.54,

    T0C +5.0, 2.33,	1.32,

    T0C+12.0, 1.66,	1.25,

    T0C+20.0, 1.25,	1.29,

    ]).reshape((-1,3))



def tau1_Teixeira(T):
@@ -424,3 +426,16 @@ def eta_d2o(T): 
    """

    return _water_viscosity(T-231.832, 885.60402e-3, 1.55255,

                            (2.9067e-8, -1.6342e-5, 2.7990e-3))



# simple fit

def rho_h2o(T):

    """returns water density in g/cm^3 T in Kelvin """

    temp = T-T0C  # convert to Celsius

    assert 5.0<=temp<=40.0, "temperature must be between 5C and 40C"

    return 1.00029 - 8.57143e-6*temp - 4.85714e-6*temp**2



def rho_d2o(T):

    """returns heavy water density in g/cm^3 T in Kelvin """

    temp = T-T0C  # convert to Celsius

    assert 5.0<=temp<=40.0, "temperature must be between 5C and 40C"

    return 1.10517 + 1.42619e-4*temp - 6.80952e-6*temp**2

pysen/revision.py

+2 −2
Original line number Diff line number Diff line
@@ -3,8 +3,8 @@ PySEN revision module 
"""

import sys

__version__  = "2.0"

__release__  = ".2b3"

__date__     = "July 3, 2025"

__release__  = ".2b4"

__date__     = "July 16, 2025"



def version(full=False):

    "get pysen version number"

test/test_physics_water_model.py

+13 −1
Original line number Diff line number Diff line
@@ -59,7 +59,7 @@ class WaterModelTestCase(unittest.TestCase): 


    def test_viscosity(self):

        "test viscosity"

        C = 273.15

        C = wm.T0C

        mili=1e-3

        self.assertAlmostEqual(wm.eta_h2o( 5.0+C)/mili, 1.51893, 5)

        self.assertAlmostEqual(wm.eta_h2o(20.0+C)/mili, 1.00203, 5)
@@ -68,6 +68,18 @@ class WaterModelTestCase(unittest.TestCase): 
        self.assertAlmostEqual(wm.eta_d2o(20.0+C)/mili, 1.25136, 5)

        self.assertAlmostEqual(wm.eta_d2o(70.0+C)/mili, 0.47435, 5)



    def test_density(self):

        "test density"

        C = wm.T0C

        exp_rho_h2o = { 10.0: 0.9997, 15.0: 0.9991, 20.0: 0.9982, 

                        25.0: 0.9970, 30.0: 0.9957, 35.0: 0.9940 }

        exp_rho_d2o = { 10.0: 1.1059, 15.0: 1.1058, 20.0: 1.1053,

                        25.0: 1.1045, 30.0: 1.1033} 

        for t in exp_rho_h2o:

            self.assertAlmostEqual(wm.rho_h2o(t+C), exp_rho_h2o[t], 4, f"rho_h2o at {t}C failed")

        for t in exp_rho_d2o:

            self.assertAlmostEqual(wm.rho_d2o(t+C), exp_rho_d2o[t], 4, f"rho_d2o at {t}C failed")



    def test_teixeira(self):

        "test Teixeira model"

        #def uev2ps(uev):