def water_specific_weight(temperature=25*unit('degC'),pressure=1*unit('atm')):
    '''    
    Return the specific weight of water in N/m3 at the specified temperature and pressure.
    
    Parameters
    ----------
    temperature : Quantity, optional
                  The temperature. Defaults to 25 degC if omitted.
    pressure    : Quantity, optional
                  The ambient pressure of the solution. 
                  Defaults to atmospheric pressure (1 atm) if omitted.
                  
    Returns
    -------
    Quantity
            The specific weight of water in N/m3.  
    
    Examples
    --------
    >>> water_specific_weight() #doctest: +ELLIPSIS
    <Quantity(9777.637025975, 'newton / meter ** 3')>
            
    See Also
    --------
    water_density
    
    '''
    spweight = water_density(temperature,pressure) * unit.g_n
    logger.info('Computed specific weight of water as %s at T=%s and P = %s' % (spweight,temperature,pressure))
    return spweight.to('N/m ** 3')

def water_viscosity_kinematic(temperature=25*unit('degC'),pressure=1*unit('atm')):
    '''
    Return the kinematic viscosity of water in m2/s = Stokes
    at the specified temperature.
    
    Parameters
    ----------
    temperature : Quantity, optional
                  The temperature. Defaults to 25 degC if omitted.
    pressure    : Quantity, optional
                  The ambient pressure of the solution. 
                  Defaults to atmospheric pressure (1 atm) if omitted.
                  
    Returns
    -------
    Quantity
            The kinematic viscosity of water in Stokes (m2/s)
    
    Examples
    --------
    >>> water_viscosity_kinematic()  #doctest: +ELLIPSIS
    <Quantity(8.899146003595295e-07, 'meter ** 2 / second')>
            
    See Also
    --------
    water_viscosity_dynamic
    water_density
    
    '''
    kviscosity = water_viscosity_dynamic(temperature,pressure) / water_density(temperature,pressure)
    logger.info('Computed kinematic viscosity of water as %s at T=%s and P = %s ' % (kviscosity,temperature,pressure)) 
    return kviscosity.to('m**2 / s')

def _debye_parameter_activity(temperature="25 degC"):
    """
    Return the constant A for use in the Debye-Huckel limiting law (base 10)
    
    Parameters
    ----------
    temperature : str Quantity, optional
                  String representing the temperature of the solution. Defaults to '25 degC' if not specified.
    
    Returns
    -------
    Quantity          The parameter A for use in the Debye-Huckel limiting law (base e)
    
    Notes
    -----
     
    The parameter A is equal to: [#]_
     
    ..  math::    
        A^{\\gamma} = {e^3 ( 2 \\pi N_A {\\rho})^{0.5} \\over (4 \\pi \\epsilon_o \\epsilon_r k T)^{1.5}}
    
    Note that this equation returns the parameter value that can be used to calculate
    the natural logarithm of the activity coefficient. For base 10, divide the
    value returned by 2.303. The value is often given in base 10 terms (0.509 at
    25 degC) in older textbooks.
    
    References
    ----------
    .. [#] Archer, Donald G. and Wang, Peiming. "The Dielectric Constant of Water \
    and Debye-Huckel Limiting Law Slopes." /J. Phys. Chem. Ref. Data/ 19(2), 1990.
        
    Examples
    --------
    >>> _debye_parameter_activity() #doctest: +ELLIPSIS
    1.17499...
    
    See Also
    --------
    _debye_parameter_osmotic
    
    """

    debyeparam = (
        unit.elementary_charge ** 3
        * (2 * math.pi * unit.avogadro_number * h2o.water_density(unit(temperature))) ** 0.5
        / (
            4
            * math.pi
            * unit.epsilon_0
            * h2o.water_dielectric_constant(unit(temperature))
            * unit.boltzmann_constant
            * unit(temperature)
        )
        ** 1.5
    )

    logger.info("Computed Debye-Huckel Limiting Law Constant A^{\\gamma} = %s at %s" % (debyeparam, temperature))
    return debyeparam.to("kg ** 0.5 / mol ** 0.5")

def water_dielectric_constant(temperature=25*unit('degC')):
    '''    
    Return the dielectric constant of water at the specified temperature.
    
    Parameters
    ----------
    temperature : Quantity, optional
                  The temperature. Defaults to 25 degC if omitted.
                  
    Returns
    -------
    float
            The dielectric constant (or permittivity) of water relative to the
            permittivity of a vacuum. Dimensionless.
    
    Notes
    -----
    This function implements a quadratic fit of measured permittivity data as
    reported in the CRC Handbook [#]_. The parameters given are valid over the
    range 273 K to 372 K. Permittivity should not be extrapolated beyond this
    range.
    
    .. math:: \\epsilon(T) = a + b T + c T^2
    
    References
    ----------
    .. [#] "Permittivity (Dielectric Constant) of Liquids." *CRC Handbook of 
            Chemistry and Physics*, 92nd ed, pp 6-187 - 6-208.
    
    Examples
    --------
    >>> water_dielectric_constant(unit('20 degC')) #doctest: +ELLIPSIS
    80.15060...
    
    Display an error if 'temperature' is outside the valid range
    
    >>> water_dielectric_constant(-5*unit('degC'))
    
     
    '''
    # do not return anything if 'temperature' is outside the range for which
    # this fit applies
    if temperature < 273 * unit('K') or temperature > 372 * unit('K'):
        logger.error('Specified temperature (%s) exceeds valid range of data. Cannot extrapolate.' % temperature.to('K'))
        return None
    
    # otherwise, calculate the dielectric constant using the quadratic fit    
    a = 0.24921e3
    b = -0.79069e0
    c = 0.72997e-3
    dielectric = a + b * temperature.to('K').magnitude + c * temperature.to('K').magnitude ** 2
    
    logger.info('Computed dielectric constant of water as %s at %s' % (dielectric,temperature))
    
    logger.debug('Computed dielectric constant of water using empirical equation given in "Permittivity (Dielectric Constant) of Liquids." CRC Handbook of Chemistry and Physics, 92nd ed, pp 6-187 - 6-208.')
    
    return dielectric

 def test_molar_conductivity_chloride(self):
     # Cl- - 76.31 x 10 ** -4 m ** 2 S / mol
     self.s1 = pyEQL.Solution([['Na+','0.001 mol/L'],['Cl-','0.001 mol/L']],temperature='25 degC')
     result = self.s1.get_molar_conductivity('Cl-').to('m**2*S/mol').magnitude
     expected = pyEQL.unit('76.31e-4 m**2 * S / mol').magnitude
     
     self.assertWithinExperimentalError(result,expected,self.tol)

 def test_molar_conductivity_magnesium(self):
     # Mg+2 - 106 x 10 ** -4 m ** 2 S / mol
     self.s1 = pyEQL.Solution([['Mg+2','0.001 mol/L'],['Cl-','0.002 mol/L']],temperature='25 degC')
     result = self.s1.get_molar_conductivity('Mg+2').to('m**2*S/mol').magnitude
     expected = pyEQL.unit('106e-4 m**2 * S / mol').magnitude
     
     self.assertWithinExperimentalError(result,expected,self.tol)

 def test_molar_conductivity_potassium(self):
     # K+ - 73.48 x 10 ** -4 m ** 2 S / mol
     self.s1 = pyEQL.Solution([['K+','0.001 mol/L'],['Cl-','0.001 mol/L']],temperature='25 degC')
     result = self.s1.get_molar_conductivity('K+').to('m**2*S/mol').magnitude
     expected = pyEQL.unit('73.48e-4 m**2 * S / mol').magnitude
     
     self.assertWithinExperimentalError(result,expected,self.tol)

 def test_molar_conductivity_hydrogen(self):
     # H+ - 349.65 x 10 ** -4 m ** 2 S / mol
     self.s1 = pyEQL.Solution(temperature='25 degC')
     result = self.s1.get_molar_conductivity('H+').to('m**2*S/mol').magnitude
     expected = pyEQL.unit('349.65e-4 m**2 * S / mol').magnitude
     
     self.assertWithinExperimentalError(result,expected,self.tol)

 def test_molar_conductivity_hydroxide(self):
     # OH- - 198 x 10 ** -4 m ** 2 S / mol
     self.s1 = pyEQL.Solution(temperature='25 degC')
     result = self.s1.get_molar_conductivity('OH-').to('m**2*S/mol').magnitude
     expected = pyEQL.unit('198e-4 m**2 * S / mol').magnitude
     
     self.assertWithinExperimentalError(result,expected,self.tol)

 def test_molar_conductivity_sulfate(self):
     # SO4-2 - 160 x 10 ** -4 m ** 2 S / mol
     self.s1 = pyEQL.Solution([['Na+','0.002 mol/L'],['SO4-2','0.001 mol/L']],temperature='25 degC')
     result = self.s1.get_molar_conductivity('SO4-2').to('m**2*S/mol').magnitude
     expected = pyEQL.unit('160.0e-4 m**2 * S / mol').magnitude
     
     self.assertWithinExperimentalError(result,expected,self.tol)

 def test_effective_pitzer_mgcl2_activity(self):
     # test the activity coefficient of MgCl2
     # corresponds to 0.515m, 1.03m, 2.58m, and 4.1m
     multiple = [1,2,5,8]
     expected=[0.5,0.5,0.67,1.15]
     
     # import the parameters database
     from pyEQL import paramsDB as db
         
     for item in range(len(multiple)):
         s1 = self.mock_seawater(multiple[item])
         Salt = pyEQL.salt_ion_match.Salt('Mg+2','Cl-')
         db.search_parameters(Salt.formula)
         param = db.get_parameter(Salt.formula,'pitzer_parameters_activity')
         alpha1 = 2
         alpha2 = 0
         molality = Salt.get_effective_molality(s1.get_ionic_strength())
         temperature = str(s1.get_temperature())
         
         activity_coefficient=pyEQL.activity_correction.get_activity_coefficient_pitzer(s1.get_ionic_strength(), \
         molality,alpha1,alpha2,param.get_value()[0],param.get_value()[1],param.get_value()[2],param.get_value()[3], \
         Salt.z_cation,Salt.z_anion,Salt.nu_cation,Salt.nu_anion,temperature)
         
         # convert the result to a rational activity coefficient
         result = activity_coefficient * (1+pyEQL.unit('0.018 kg/mol')*s1.get_total_moles_solute()/s1.get_solvent_mass())
         #print(result,expected[item])
         self.assertWithinExperimentalError(result,expected[item],self.tol) 

def _debye_parameter_B(temperature='25 degC'):
    '''
    Return the constant B used in the extended Debye-Huckel equation
    
    Parameters
    ----------
    temperature : str Quantity, optional
                  String representing the temperature of the solution. Defaults to '25 degC' if not specified.
    
    Notes
    -----
    The parameter B is equal to: [#]_

    .. math:: B = ( {8 \\pi N_A e^2 \\over 1000 \\epsilon k T} ) ^ {1 \\over 2}
    
    .. [#] Bockris and Reddy. /Modern Electrochemistry/, vol 1. Plenum/Rosetta, 1977, p.210.    
    
    Examples
    --------
    >>> _debye_parameter_B() #doctest: +ELLIPSIS
    0.3291...
    
    '''
    # TODO - fix this and resolve units
    param_B = ( 8 * math.pi * unit.avogadro_number * unit.elementary_charge ** 2 
    / (h2o.water_density(unit(temperature)) * unit.epsilon_0 * h2o.water_dielectric_constant(unit(temperature)) * unit.boltzmann_constant * unit(temperature)) )** 0.5
    return param_B.to_base_units()

def water_density(temperature=25*unit('degC'),pressure=1*unit('atm')):
    # TODO add pressure??
    # TODO more up to date equation??
    '''    
    Return the density of water in kg/m3 at the specified temperature and pressure.
    
    Parameters
    ----------
    temperature : float or int, optional
                  The temperature in Celsius. Defaults to 25 degrees if not specified.
    pressure    : float or int, optional
                  The ambient pressure of the solution in Pascals (N/m2). 
                  Defaults to atmospheric pressure (101325 Pa) if not specified.
    
    Returns
    -------
    float
            The density of water in kg/m3.
    
    Notes
    -----
    Based on the following empirical equation reported in [#]_
    
    
    .. math:: \\rho_W = 999.65 + 0.20438 T - 6.1744e-2 T ^ {1.5}
    
    Where :math:`T` is the temperature in Celsius.
    
    
    .. [#] Sohnel, O and Novotny, P. *Densities of Aqueous Solutions of Inorganic Substances.* Elsevier Science, Amsterdam, 1985.
    
    Examples
    --------
    >>> water_density(25*unit('degC')) #doctest: +ELLIPSIS
    <Quantity(997.0415, 'kilogram / meter ** 3')>
    
    '''
    # calculate the magnitude
    density = 999.65 + 0.20438 * temperature.to('degC').magnitude - 6.1744e-2 * temperature.to('degC').magnitude ** 1.5
    # assign the proper units
    density = density  * unit('kg/m**3')
    logger.info('Computed density of water as %s at T= %s and P = %s' % (density,temperature,pressure))
    logger.debug('Computed density of water using empirical relation in Sohnel and Novotny, "Densities of Aqueous Solutions of Inorganic Substances," 1985' )
    return density.to('kg/m**3')

def _debye_parameter_volume(temperature='25 degC'):
    '''
    Return the constant A_V, the Debye-Huckel limiting slope for apparent
    molar volume.
    
    Parameters
    ----------
    temperature : str Quantity, optional
                  String representing the temperature of the solution. Defaults to '25 degC' if not specified.
    
    Notes
    -----
    Takes the value 1.8305 cm ** 3 * kg ** 0.5 /  mol ** 1.5 at 25 C.
    This constant is calculated according to: [#]_

     .. math:: A_V = -2 A_{\\phi} R T [ {3 \\over \\epsilon} {{\\partial \\epsilon \\over \\partial p} \
     } - {{1 \\over \\rho}{\\partial \\rho \\over \\partial p} }]
     
    NOTE: at this time, the term in brackets (containing the partial derivatives) is approximate.
    These approximations give the correct value of the slope at 25 degC and 
    produce estimates with less than 10% error between 0 and 60 degC.
     
    The derivative of epsilon with respect to pressure is assumed constant (for atmospheric pressure)
    at -0.01275 1/MPa. Note that the negative sign does not make sense in light
    of real data, but is required to give the correct result.
     
    The second term is equivalent to the inverse of the bulk modulus of water, which
    is taken to be 2.2 GPa. [#]_
    
    References
    ----------
    .. [#] Archer, Donald G. and Wang, Peiming. "The Dielectric Constant of Water \
    and Debye-Huckel Limiting Law Slopes." /J. Phys. Chem. Ref. Data/ 19(2), 1990.
        
    .. [#] http://hyperphysics.phy-astr.gsu.edu/hbase/permot3.html
    
    Examples
    --------
    TODO
  
    
    See Also
    --------
    _debye_parameter_osmotic
    
    '''
    
    # TODO - add partial derivatives to calculation
    epsilon = h2o.water_dielectric_constant(unit(temperature))
    dedp = unit('-0.01275 1/MPa')
    result = -2 * _debye_parameter_osmotic(temperature) * unit.R * unit(temperature) * \
    (3 / epsilon * dedp - 1/unit('2.2 GPa'))
    #result = unit('1.898 cm ** 3 * kg ** 0.5 /  mol ** 1.5')
    
    if unit(temperature) != unit('25 degC'):
        logger.warning('Debye-Huckel limiting slope for volume is approximate when T is not equal to 25 degC')
    
    logger.info('Computed Debye-Huckel Limiting Slope for volume A^V = %s at %s' % (result,temperature))
    
    return result.to('cm ** 3 * kg ** 0.5 /  mol ** 1.5')

def get_activity_coefficient_davies(ionic_strength, formal_charge=1, temperature="25 degC"):
    """
    Return the activity coefficient of solute in the parent solution according to the Davies equation.
    
    Parameters
    ----------
    formal_charge : int, optional      
                    The charge on the solute, including sign. Defaults to +1 if not specified.
    ionic_strength : Quantity
                     The ionic strength of the parent solution, mol/kg
    temperature : str Quantity, optional
                     String representing the temperature of the solution. Defaults to '25 degC' if not specified.
                  
    Returns
    -------
    Quantity
         The mean molal (mol/kg) scale ionic activity coefficient of solute, dimensionless.

    See Also
    --------
    _debye_parameter_activity
    
    Notes
    -----
    Activity coefficient is calculated according to: [#]_

    .. math:: \\ln \\gamma = A^{\\gamma} z_i^2 ({\sqrt I \\over (1 + \sqrt I)} + 0.2 I)
    
    Valid for 0.1 < I < 0.5
    
    References
    ----------
    .. [#] Stumm, Werner and Morgan, James J. Aquatic Chemistry, 3rd ed, 
           pp 103. Wiley Interscience, 1996.
    
    """
    # check if this method is valid for the given ionic strength
    if not ionic_strength.magnitude <= 0.5 and ionic_strength.magnitude >= 0.1:
        logger.warning("Ionic strength exceeds valid range of the Davies equation")

    # the units in this empirical equation don't work out, so we must use magnitudes
    log_f = (
        -_debye_parameter_activity(temperature).magnitude
        * formal_charge ** 2
        * (ionic_strength.magnitude ** 0.5 / (1 + ionic_strength.magnitude ** 0.5) - 0.2 * ionic_strength.magnitude)
    )

    return math.exp(log_f) * unit("1 dimensionless")

def get_activity_coefficient_guntelberg(ionic_strength,formal_charge=1,temperature='25 degC'):
    '''
    Return the activity coefficient of solute in the parent solution according to the Guntelberg approximation.
    
    Parameters
    ----------
    formal_charge : int, optional      
                    The charge on the solute, including sign. Defaults to +1 if not specified.
    ionic_strength : Quantity
                     The ionic strength of the parent solution, mol/kg
    temperature : str Quantity, optional
                     String representing the temperature of the solution. Defaults to '25 degC' if not specified.
                  
    Returns
    -------
    Quantity
         The mean molal (mol/kg) scale ionic activity coefficient of solute, dimensionless.
         
    See Also
    --------
    _debye_parameter_activity
    
    Notes
    ------
    Activity coefficient is calculated according to: [#]_ 

    .. math:: \\ln \\gamma = A^{\\gamma} z_i^2 {\sqrt I \\over (1 + \sqrt I)}
    
    Valid for I < 0.1
    
    References
    ----------
    .. [#] Stumm, Werner and Morgan, James J. Aquatic Chemistry, 3rd ed, 
           pp 103. Wiley Interscience, 1996.
    
    '''
    # check if this method is valid for the given ionic strength
    if not ionic_strength.magnitude <= 0.1:
        logger.warning('Ionic strength exceeds valid range of the Guntelberg approximation')
    
    log_f = - _debye_parameter_activity(temperature) * formal_charge ** 2 * ionic_strength ** 0.5 / (1+ionic_strength.magnitude ** 0.5)

    return math.exp(log_f) * unit('1 dimensionless')

def get_activity_coefficient_debyehuckel(ionic_strength, formal_charge=1, temperature="25 degC"):
    """
    Return the activity coefficient of solute in the parent solution according to the Debye-Huckel limiting law.
    
    Parameters
    ----------
    formal_charge : int, optional      
                    The charge on the solute, including sign. Defaults to +1 if not specified.
    ionic_strength : Quantity
                     The ionic strength of the parent solution, mol/kg
    temperature : str Quantity, optional
                     String representing the temperature of the solution. Defaults to '25 degC' if not specified.
                  
    Returns
    -------
    Quantity
         The mean molal (mol/kg) scale ionic activity coefficient of solute, dimensionless.

    See Also
    --------
    _debye_parameter_activity
    
    Notes
    -----
    Activity coefficient is calculated according to: [#]_ 

    .. math:: \\ln \\gamma = A^{\\gamma} z_i^2 \sqrt I
    
    Valid only for I < 0.005
    
    References
    ----------
    .. [#] Stumm, Werner and Morgan, James J. Aquatic Chemistry, 3rd ed, 
           pp 103. Wiley Interscience, 1996.
    
    """
    # check if this method is valid for the given ionic strength
    if not ionic_strength.magnitude <= 0.005:
        logger.warning("Ionic strength exceeds valid range of the Debye-Huckel limiting law")

    log_f = -_debye_parameter_activity(temperature) * formal_charge ** 2 * ionic_strength ** 0.5

    return math.exp(log_f) * unit("1 dimensionless")

def adjust_temp_pitzer(c1,c2,c3,c4,c5,temp,temp_ref=unit('298.15 K')):
    '''
    Calculate a parameter for th e Pitzer model based on temperature-dependent
    coefficients c1,c2,c3,c4,and c5.
    
    Parameters
    ----------
    c1, c2, c3, c4, c5: float
                Temperature-dependent coefficients for the pitzer parameter of 
                interest.
    temp: Quantity
                The temperature at which the Pitzer parameter is to be calculated
    temp_ref: Quantity, optional
                The reference temperature on which the parameters are based.
                298.15 K if omitted.
    
    As described in the PHREEQC documentation
    
    
    '''
    pitzer_param = c1 + c2 * (1/temp + 1/temp_ref) + c2 * math.log(temp/temp_ref) \
    + c3 * (temp - temp_ref) + c4 * (temp ** 2 - temp_ref ** 2) + c5 * (temp ** -2 - temp_ref ** -2)
    
    return pitzer_param

def water_viscosity_dynamic(temperature=25*unit('degC'),pressure=1*unit('atm')):
    '''
    Return the dynamic (absolute) viscosity of water in N-s/m2 = Pa-s = kg/m-s
    at the specified temperature.
    
    Parameters
    ----------
    temperature : Quantity, optional
                  The temperature. Defaults to 25 degC if omitted.
    pressure    : Quantity, optional
                  The ambient pressure of the solution. 
                  Defaults to atmospheric pressure (1 atm) if omitted.
    
    Returns
    -------
    Quantity 
                The dynamic (absolute) viscosity of water in N-s/m2 = Pa-s = kg/m-s
                  
    Notes
    -----
    Implements the international equation for viscosity of water as specified by NIST [#]_
    
    Valid for 273 < temperature < 1073 K and 0 < pressure < 100,000,000 Pa
    
    References
    ----------
    .. [#] Sengers, J.V. "Representative Equations for the Viscosity of Water Substance." 
        *J. Phys. Chem. Ref. Data* 13(1), 1984.http://www.nist.gov/data/PDFfiles/jpcrd243.pdf
    
    Examples
    --------
    >>> water_viscosity_dynamic(20*unit('degC')) #doctest: +ELLIPSIS
    <Quantity(0.000998588610804179, 'kilogram / meter / second')>
    >>> water_viscosity_dynamic(unit('100 degC'),unit('25 MPa')) #doctest: +ELLIPSIS
    <Quantity(0.00028165034364318573, 'kilogram / meter / second')>
    >>> water_viscosity_dynamic(25*unit('degC'),0.1*unit('MPa')) #doctest: +ELLIPSIS
    <Quantity(0.0008872817880143659, 'kilogram / meter / second')>
    
    #TODO - check these again after I implement pressure-dependent density function
    
    '''
    # generate warnings if temp or pressure are outside valid range of equation
    if temperature < 273 * unit('K') or temperature > 1073 * unit('K'):
        logger.error('Specified temperature (%s) exceeds valid range of NIST equation for viscosity of water. Cannot extrapolate.' % temperature)
        return None
        
    if pressure < 0 * unit('Pa') or pressure > 100000000 * unit ('Pa'):
        logger.error('Specified pressure (%s) exceeds valid range of NIST equation for viscosity of water. Cannot extrapolate.' % pressure)
        return None
    
    # calculate dimensionless temperature and pressure
    T_star = 647.27 #K
    P_star = 22115000 #Pa
    rho_star = 317.763 #kg/m3
    
    T_bar = temperature.to('K').magnitude / T_star
    P_bar = pressure.to('Pa').magnitude / P_star
    rho_bar = water_density(temperature,pressure).magnitude / rho_star
    
    # calculate the first function, mu_o
    mu_star = 1e-6 #Pa-s
    a = [0.0181583,0.0177624,0.0105287,-0.0036477]
    sum_o = 0
    mu_temp = 0
    for index in range(len(a)):
        sum_o += a[index] * T_bar ** -index
    
    mu_o = mu_star * math.sqrt(T_bar) / sum_o
    
    # calculate the second fucntion, mu_1
    b=[[0.501938,0.235622,-0.274637,0.145831,-0.0270448],[0.162888,0.789393,-0.743539,0.263129,-0.0253093],[-0.130356,0.673665,-0.959456,0.347247,-0.0267758],[0.907919,1.207552,-0.687343,0.213486,-0.0822904],[-0.551119,0.0670665,-0.497089,0.100754,0.0602253],[0.146543,-0.0843370,0.195286,-0.032932,-0.0202595]]
    mu_1 = 0
    
    for i in range(len(b)):
        for j in range(len(b[i])):
            mu_temp += rho_bar * b[i][j] * (1/T_bar -1 ) ** i * (rho_bar -1) ** j
    
    mu_1 = math.exp(mu_temp)
    # multiply the functions to return the viscosity
    viscosity = mu_o * mu_1 *unit('kg/m/s')
    
    logger.info('Computed dynamic (absolute) viscosity of water as %s at T=%s and P = %s'  % (viscosity,temperature,pressure)) 
    
    logger.debug('Computed dynamic (absolute) viscosity of water using empirical NIST equation described in Sengers, J.V. "Representative Equations for the Viscosity of Water Substance." J. Phys. Chem. Ref. Data 13(1), 1984.')
    
    return viscosity.to('kg/m/s')