Homec4science
Diffusion Gtools (master) Browse

gtools/Gtools/pscripts/pJeansTemperature
8f16877662e4master

Last Change
Raw File

pJeansTemperature
View Options

#!/usr/bin/env python
from optparse import OptionParser
import Ptools as pt
from pNbody import *
from pNbody import units
from pNbody import ctes
from pNbody import thermodyn
import string
from scipy import optimize
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = pt.add_postscript_options(parser)
parser = pt.add_limits_options(parser)
parser.add_option("--Nsph",
action="store",
dest="Nsph",
type="float",
default = 1.,
help="Number of neighbouring particles for SPH",
metavar=" FLOAT")
parser.add_option("--NJeans",
action="store",
dest="NJeans",
type="float",
default = 10.,
help="Number of particles sampling a Jeans Mass",
metavar=" FLOAT")
parser.add_option("--ParticleMass",
action="store",
dest="ParticleMass",
type="float",
default = 1,
help="Particle mass in Msol",
metavar=" FLOAT")
(options, args) = parser.parse_args()
files = args
return files,options
#######################################
# MakePlot
#######################################
def MakePlot(files,opt):
# some inits
colors = pt.Colors(n=len(files))
datas = []
# some constants
gamma = 5/3.
G = ctes.GRAVITY.value
PROTONMASS = ctes.PROTONMASS.value
BOLTZMANN = ctes.BOLTZMANN.value
# input units
Unit_atom = ctes.PROTONMASS.into(units.cgs)*units.Unit_g
Unit_atom.set_symbol('atom')
in_units = units.UnitSystem('local',[units.Unit_cm,Unit_atom,units.Unit_s,units.Unit_K])
in_units_mass = units.UnitSystem('local',[units.Unit_cm,units.Unit_Msol,units.Unit_s,units.Unit_K])
# cgs
cgs_units = units.UnitSystem('local',[units.Unit_cm,units.Unit_g,units.Unit_s,units.Unit_K])
# Density : assume units in a/cc
dx = (opt.xmax - opt.xmin)/1000.
Density = 10**arange(opt.xmin,opt.xmax,dx)
DensityCGS = Density * in_units.convertionFactorTo(cgs_units.UnitDensity)
# Mass
ParticleMassCGS = opt.ParticleMass * in_units_mass.convertionFactorTo(cgs_units.UnitMass)
# deflauts parameters in cgs
pars = {"k":BOLTZMANN,"mh":PROTONMASS,"mu":2,"gamma":gamma,"G":G}
EnergyInt = (36./pi**5)**(1/3.) * (float(opt.NJeans)/opt.Nsph)**(2/3.) * G * gamma**(-1) * (gamma-1)**(-1) * ParticleMassCGS**(2/3.) * DensityCGS**(1/3.)
Temperature = thermodyn.Tru(Density,EnergyInt,pars)
opt.log=None
data = pt.DataPoints(log10(Density),log10(Temperature))
datas.append(data)
for d in datas:
pt.plot(d.x,d.y,color=d.color)
# set limits and draw axis
xmin,xmax,ymin,ymax = pt.SetLimitsFromDataPoints(opt.xmin,opt.xmax,opt.ymin,opt.ymax,datas,opt.log)
xlabel = r'$\rm{log\,Density}\,\left[ \rm{atom/cm^3} \right]$'
ylabel = r'$\rm{log\,Jeans\,Temperature}\,\left[ \rm{K} \right]$'
# plot axis
pt.SetAxis(xmin,xmax,ymin,ymax,log=opt.log)
pt.xlabel(xlabel,fontsize=pt.labelfont)
pt.ylabel(ylabel,fontsize=pt.labelfont)
pt.grid(False)
########################################################################
# MAIN
########################################################################
if __name__ == '__main__':
files,opt = parse_options()
pt.InitPlot(files,opt)
pt.pcolors
MakePlot(files,opt)
pt.EndPlot(files,opt)
c4science · Help