Homec4science
Diffusion Gtools (master) Browse

gtools/Gtools/vanderwaals.py
c2cb5b12a7e2master

Last Change
Raw File

vanderwaals.py
View Options

from numpy import *
from pNbody.ctes import *
####################################################################################################################################
#
# VAN DER WAALS THERMODYNAMIC RELATIONS
#
####################################################################################################################################
# deflauts parameters in cgs
defaultpars = {"k":BOLTZMANN,"mh":PROTONMASS,"mu":2,"gamma":5/3.,"a":AV,"b":BV}
###################
def Prt(rho,T,pars=defaultpars):
###################
'''
P(rho,T)
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return k*T/(mumh/rho-b) - a*(rho/mumh)**2
###################
def Trp(rho,P,pars=defaultpars):
###################
'''
T(rho,P)
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return (mumh/rho-b)/k * ( P + a*(rho/mumh)**2 )
###################
def Art(rho,T,pars=defaultpars):
###################
'''
A(rho,T)
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return k/mumh * (rho*mumh/(mumh-b*rho))**(1.-gamma) * T
###################
def Tra(rho,A,pars=defaultpars):
###################
'''
T(rho,A)
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return mumh/k * (rho*mumh/(mumh-b*rho))**(gamma-1.) * A
###################
def Urt(rho,T,pars=defaultpars):
###################
'''
U(rho,T)
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return 1./(gamma-1.)* k*T/mumh - a*rho/mumh**2
###################
def Tru(rho,U,pars=defaultpars):
###################
'''
T(rho,U)
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return (gamma-1.)* mumh/k * ( U + a*rho/mumh**2 )
###################
def Pra(rho,A,pars=defaultpars):
###################
'''
P(rho,A)
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return ((rho*mumh)/(mumh-b*rho))**gamma * A - a*(rho/mumh)**2
###################
def Arp(rho,P,pars=defaultpars):
###################
'''
A(rho,P)
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return ((mumh-b*rho)/(rho*mumh))**gamma * ( P + a*(rho/mumh)**2 )
###################
def Pru(rho,U,pars=defaultpars):
###################
'''
P(rho,U)
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return (gamma-1.) * ((rho*mumh)/(mumh-b*rho)) * ( U + a*rho/mumh**2 ) - a*(rho/mumh)**2
###################
def Urp(rho,P,pars=defaultpars):
###################
'''
U(rho,P)
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return 1./(gamma-1.) * ((mumh-b*rho)/(rho*mumh)) * ( P + a*(rho/mumh)**2 ) - a*rho/mumh**2
###################
def Ura(rho,A,pars=defaultpars):
###################
'''
U(rho,A)
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return 1./(gamma-1.) * ((rho*mumh)/(mumh-b*rho))**(gamma-1.) * A - a*rho/mumh**2
###################
def Aru(rho,U,pars=defaultpars):
###################
'''
A(rho,U)
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return (gamma-1.) * ((mumh-b*rho)/(rho*mumh))**(gamma-1.) * ( U + a*rho/mumh**2 )
####################################################################################################################################
# SOME CRITICAL VALUES
####################################################################################################################################
###################
def Rl(pars=defaultpars):
###################
'''
Limit density
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
if b!=0:
return mu*mh/b
else:
return 0.
###################
def Rc(pars=defaultpars):
###################
'''
Critical density
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
if b!=0:
return mu*mh/(3*b)
else:
return 0.
###################
def Pc(pars=defaultpars):
###################
'''
Critical pressure
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
if b!=0:
return a/(27*b**2)
else:
return 0.
###################
def Tc(pars=defaultpars):
###################
'''
Critical tempertature
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
if b!=0:
return 8*a/(27*k*b)
else:
return 0.
###################
def Aneg(rho,pars=defaultpars):
###################
'''
A(rho) below wich, the pressure is negative
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return a*(rho/mumh)**2 * ((mumh-b*rho)/(rho*mumh))**gamma
###################
def Uneg(rho,pars=defaultpars):
###################
'''
U(rho) below wich, the pressure is negative
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return a*(rho/mumh)**2 * 1./(gamma-1.) * ((mumh-b*rho)/(rho*mumh)) - a*rho/mumh**2
###################
def Tneg(rho,pars=defaultpars):
###################
'''
T(rho) below wich, the pressure is negative
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return a*(rho/mumh)**2 * (mumh/rho - b)/k
###################
def Tmin(rho,pars=defaultpars):
###################
'''
T(rho) below wich, the specific energy must be negative
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return (gamma-1.)* mumh/k * (a*rho/mumh**2 )
###################
def Amin(rho,pars=defaultpars):
###################
'''
A(rho) below wich, the specific energy must be negative
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
mumh = mu*mh
return (gamma-1.) * ((mumh-b*rho)/(rho*mumh))**(gamma-1.) * (a*rho/mumh**2 )
###################
def Extremas(T,pars=defaultpars):
###################
'''
Give the 3 local extremas (if they exists) of P(rho,T) as a function of T
'''
k = pars['k']
mh = pars['mh']
mu = pars['mu']
gamma = pars['gamma']
a = pars['a']
b = pars['b']
if b!=0:
Tc = 8*a/(27*k*b)
else:
return zeros(len(T)).real, zeros(len(T)).real, zeros(len(T)).real, T!=T
p = -2*mu*mh/b
q = (mu*mh/b)**2
r = (-k*T*(mu*mh)**3)/(2*a*b**2)
aa = (1/3. )*(3*q -p**2)
bb = (1/27.)*(2*p**3 -9*p*q +27*r)
C = sqrt( (bb**2)/4. + (aa**3)/27. +0j)
A = (-bb/2.+C)**(1/3.)
B = (-bb/2.-C)**(1/3.)
# intermediate solution
x1 = A+B
x2 = -(A+B)/2 + (A-B)/2.* sqrt(-3+0j)
x3 = -(A+B)/2 - (A-B)/2.* sqrt(-3+0j)
# final solution
rho1 = x1 - p/3.
rho2 = x2 - p/3.
rho3 = x3 - p/3.
# real part
return rho2.real, rho3.real, rho1.real, (T<Tc)
c4science · Help