Homec4science
Diffusion Gtools (master) Browse

gtools/Gtools/__init__.py
c2cb5b12a7e2master

Last Change
Raw File

__init__.py
View Options

'''
Some usefull function used to deal with Gadget-type files
(or sometimes other types)
'''
from numpy import *
from pNbody import ctes
from pNbody import units
from io import *
from pNbody import myNumeric
try:
import SM
except:
pass
####################################################################################################################################
#
# UNITS CONVERSION
#
####################################################################################################################################
#################################
def Length2cgs(lenght,units):
#################################
'''
Convert lenght (in user units) into cgs units
'''
UnitLength_in_cm = units[0]
UnitMass_in_g = units[1]
UnitVelocity_in_cm_per_s = units[2]
return lenght*UnitLength_in_cm
#################################
def Mass2cgs(mass,units):
#################################
'''
Convert mass (in user units) into cgs units
'''
UnitLength_in_cm = units[0]
UnitMass_in_g = units[1]
UnitVelocity_in_cm_per_s = units[2]
return mass*UnitMass_in_g
#################################
def Velocity2cgs(velocity,units):
#################################
'''
Convert velocity (in user units) into cgs units
'''
UnitLength_in_cm = units[0]
UnitMass_in_g = units[1]
UnitVelocity_in_cm_per_s = units[2]
return velocity*UnitVelocity_in_cm_per_s
def Density2cgs(density,units):
'''
Convert density (in user units) into cgs units
'''
UnitLength_in_cm = units[0]
UnitMass_in_g = units[1]
UnitVelocity_in_cm_per_s = units[2]
return density*UnitMass_in_g/UnitLength_in_cm**3
#################################
def Energy2cgs(energy,units):
#################################
'''
Convert energy (in user units) into cgs units
'''
UnitLength_in_cm = units[0]
UnitMass_in_g = units[1]
UnitVelocity_in_cm_per_s = units[2]
return energy*UnitMass_in_g*UnitVelocity_in_cm_per_s**2
#################################
def EnergySpec2cgs(energy,units):
#################################
'''
Convert energy sepc (in user units) into cgs units
'''
UnitLength_in_cm = units[0]
UnitMass_in_g = units[1]
UnitVelocity_in_cm_per_s = units[2]
return energy*UnitVelocity_in_cm_per_s**2
#################################
def Time2cgs(energy,units):
#################################
'''
Convert time (in user units) into cgs units
'''
UnitLength_in_cm = units[0]
UnitMass_in_g = units[1]
UnitVelocity_in_cm_per_s = units[2]
return time*UnitLength_in_cm/UnitVelocity_in_cm_per_s
#################################
def Pressure2cgs(pressure,units):
#################################
'''
Convert pressure (in user units) into cgs units
'''
UnitLength_in_cm = units[0]
UnitMass_in_g = units[1]
UnitVelocity_in_cm_per_s = units[2]
UnitTime_in_s = UnitLength_in_cm / UnitVelocity_in_cm_per_s
return pressure*UnitMass_in_g / UnitLength_in_cm / UnitTime_in_s**2
#################################
def Set_SystemUnits_From_Params(params):
#################################
'''
return a system of units from given parameters
params is a dictionary that must constains at least
params['UnitVelocity_in_cm_per_s']
params['UnitMass_in_g']
params['UnitLength_in_cm']
'''
UnitVelocity_in_cm_per_s = params['UnitVelocity_in_cm_per_s']
UnitMass_in_g = params['UnitMass_in_g']
UnitLength_in_cm = params['UnitLength_in_cm']
UnitTime_in_s = UnitLength_in_cm / UnitVelocity_in_cm_per_s
# now from the params, define a system of units
Unit_lenght = units.Unit_cm * UnitLength_in_cm
Unit_mass = units.Unit_g * UnitMass_in_g
Unit_time = units.Unit_s * UnitTime_in_s
localsystem = units.UnitSystem('local',[Unit_lenght, Unit_mass, Unit_time , units.Unit_K, units.Unit_mol, units.Unit_C])
return localsystem
####################################################################################################################################
#
# OPTIONS FACILITIES
#
####################################################################################################################################
#################################
def add_limits_options(parser,xmin=None,xmax=None,ymin=None,ymax=None,zmin=None,zmax=None):
#################################
'''
This function allow to add limits options to a parser object
'''
parser.add_option("--xmin",
action="store",
dest="xmin",
type="float",
default=xmin,
help="min value in x")
parser.add_option("--xmax",
action="store",
dest="xmax",
type="float",
default=xmax,
help="max value in x")
parser.add_option("--ymin",
action="store",
dest="ymin",
type="float",
default=ymin,
help="min value in y")
parser.add_option("--ymax",
action="store",
dest="ymax",
type="float",
default=ymax,
help="max value in y")
parser.add_option("--zmin",
action="store",
dest="zmin",
type="float",
default=zmin,
help="min value in z")
parser.add_option("--zmax",
action="store",
dest="zmax",
type="float",
default=zmax,
help="max value in z")
return parser
#################################
def add_color_options(parser):
#################################
'''
This function allow to add color options to a parser object
'''
parser.add_option("-c",
action="store",
dest="colors",
type="string",
default = None,
help="colors",
metavar=" 0,64,192")
return parser
#################################
def add_labels_options(parser):
#################################
'''
This function allow to add labels options to a parser object
'''
parser.add_option("--labels",
action="store",
dest="labels",
type="string",
default = None,
help="labels",
metavar=''' "['name1','name2',...]" ''')
parser.add_option("--labx",
action="store",
dest="labx",
type="float",
default = 0.7,
help="box labels position x",
metavar="FLOAT")
parser.add_option("--laby",
action="store",
dest="laby",
type="float",
default = 0.7,
help="box labels position y",
metavar="FLOAT")
parser.add_option("--labdx",
action="store",
dest="labdx",
type="float",
default = 0.85,
help="box labels width : width = (1-labx)*labdx",
metavar="FLOAT")
parser.add_option("--labdy",
action="store",
dest="labdy",
type="float",
default = 0.85,
help="box labels height : height = (1-laby)*labdy",
metavar="FLOAT")
parser.add_option("--labex",
action="store",
dest="labex",
type="float",
default = 0.1,
help="labels offset in x with respect to the label box dxlab = dxbox*labex",
metavar="FLOAT")
parser.add_option("--labbox",
action="store",
dest="labbox",
type="int",
default = 1,
help="draw labels box or not",
metavar="INT")
return parser
#################################
def add_postscript_options(parser,xmin=None,xmax=None,ymin=None,ymax=None,zmin=None,zmax=None):
#################################
'''
This function allow to add postscript options to a parser object
'''
parser.add_option("-p",
action="store",
dest="ps",
type="string",
default = None,
help="postscript filename",
metavar=" FILE")
return parser
#################################
def add_units_options(parser):
#################################
'''
This function allow to add postscript options to a parser object
'''
parser.add_option("--UnitLength_in_cm",
action="store",
dest="UnitLength_in_cm",
type="float",
default=None,
help="UnitLength in cm")
parser.add_option("--UnitMass_in_g",
action="store",
dest="UnitMass_in_g",
type="float",
default=None,
help="UnitMass in g")
parser.add_option("--UnitVelocity_in_cm_per_s",
action="store",
dest="UnitVelocity_in_cm_per_s",
type="float",
default=None,
help="UnitVelocity in cm per s")
parser.add_option("--param",
action="store",
dest="GadgetParameterFile",
type="string",
default = None,
help="Gadget parameter file",
metavar=" FILE")
return parser
#################################
def add_log_options(parser):
#################################
'''
This function allow to add postscript options to a parser object
'''
parser.add_option("--log",
action="store",
dest="log",
type="string",
default=None,
help="Use log to plot values (x,y,xy)")
return parser
#################################
def add_reduc_options(parser):
#################################
'''
This function allow to reduc the number of particles
'''
parser.add_option("--reduc",
action="store",
dest="reduc",
type="int",
default=None,
help="reduc from a factor n")
return parser
#################################
def add_histogram_options(parser):
#################################
'''
This function allow to add histogram options to a parser object
'''
parser.add_option("--histogram",
action="store",
dest="histogram",
type="string",
default = 'none',
help="add histogram (none,add,only)")
return parser
#################################
def add_gas_options(parser):
#################################
'''
This function allow to give gas properties
'''
parser.add_option("--gamma",
action="store",
dest="gamma",
type="float",
default=5/3.,
help="adiabatic index")
parser.add_option("--mu",
action="store",
dest="mu",
type="float",
default=2,
help="mean molecular mass")
parser.add_option("--av",
action="store",
dest="av",
type="float",
default=0,
help="Av consante (in cgs)")
parser.add_option("--bv",
action="store",
dest="bv",
type="float",
default=0,
help="Bv consante (in cgs)")
return parser
#################################
def add_center_options(parser):
#################################
'''
This function allow to center the model
'''
parser.add_option("--center",
action="store",
dest="center",
type="string",
default=None,
help="center the model (histocenter,hdcenter)")
return parser
#################################
def add_select_options(parser):
#################################
'''
This function allow to select particles from the model
'''
parser.add_option("--select",
action="store",
dest="select",
type="string",
default='gas',
help="select particles from the model ('gas','sph','sticky',...)")
return parser
#################################
def add_cmd_options(parser):
#################################
'''
This function allow to execute a command on the model
'''
parser.add_option("--cmd",
action="store",
dest="cmd",
type="string",
default='None',
help="python command 'nb = nb.selectc((nb.T()>10))'")
return parser
#################################
def add_display_options(parser):
#################################
'''
This function allow to display the model
'''
parser.add_option("--display",
action="store",
dest="display",
type="string",
default='None',
help="display the model")
return parser
#################################
def add_ptype_options(parser):
#################################
'''
This function allow to choose point type
'''
parser.add_option("--ptype",
action="store",
dest="ptype",
type="string",
default=None,
help="ptype option (10,3)")
return parser
#################################
def get_units_options(options):
#################################
'''
This function allow to extract units options from option object
'''
try:
UnitLength_in_cm = options.UnitLength_in_cm
except:
UnitLength_in_cm = None
try:
UnitMass_in_g = options.UnitMass_in_g
except:
UnitMass_in_g = None
try:
UnitVelocity_in_cm_per_s = options.UnitVelocity_in_cm_per_s
except:
UnitVelocity_in_cm_per_s = None
try:
GadgetParameterFile = options.GadgetParameterFile
except:
GadgetParameterFile = None
if GadgetParameterFile != None:
params = io.read_params(GadgetParameterFile)
UnitLength_in_cm = params['UnitLength_in_cm']
UnitMass_in_g = params['UnitMass_in_g']
UnitVelocity_in_cm_per_s = params['UnitVelocity_in_cm_per_s']
if UnitLength_in_cm==None or UnitMass_in_g==None or UnitVelocity_in_cm_per_s==None:
UnitLength_in_cm=1.0
UnitMass_in_g=1.0
UnitVelocity_in_cm_per_s=1.0
print "One of the basic units is not set. Assuming cgs units."
return (UnitLength_in_cm,UnitMass_in_g,UnitVelocity_in_cm_per_s)
#################################
def Set_SystemUnits_From_Options(options):
#################################
'''
This function allow to define a system of units from units options
'''
UnitLength_in_cm,UnitMass_in_g,UnitVelocity_in_cm_per_s = get_units_options(options)
UnitTime_in_s = UnitLength_in_cm / UnitVelocity_in_cm_per_s
# now from the params, define a system of units
Unit_lenght = units.Unit_cm * UnitLength_in_cm
Unit_mass = units.Unit_g * UnitMass_in_g
Unit_time = units.Unit_s * UnitTime_in_s
localsystem = units.UnitSystem('local',[Unit_lenght, Unit_mass, Unit_time , units.Unit_K, units.Unit_mol, units.Unit_C])
return localsystem
####################################################################################################################################
#
# GRAPHIC FACILITIES
#
####################################################################################################################################
#################################
def Graph_Init(ps=None):
#################################
if ps==None:
g = SM.plot("x11 -bg white -fg black ")
else:
g = SM.plot("postencap %s"%ps)
return g
#################################
def Graph_SetColorsForFiles(files,cols=None):
#################################
# set colors
colors = {}
i = 0
for file in files:
if cols!=None:
colors[file]=cols[i]
else:
colors[file] = i*255/len(files)
i = i + 1
return colors
#################################
def Graph_SetLabelsForFiles(files,labels=None,ctypes=None,ltypes=None,labx=0.75,laby=0.7,labdx=0.85,labdy=0.85,labex=0.1,fdxline=0.3,labbox=1):
#################################
lab_x = labx
lab_y = laby
lab_dx = (1-lab_x)*labdx
lab_dy = (1-lab_y)*labdy
lab_box = labbox
lab_ex = lab_dx*labex
nf = len(files)
lab_oy0 = float(lab_dy)/(nf+1) * arange(nf,0,-1)
i=0
lab_oy = {}
lab_txt = {}
lab_ltype = {}
lab_ctype = {}
if labels != None:
exec("labels=%s"%(labels))
lab_show = 1
for file in files:
lab_oy[file] = lab_oy0[i]
lab_txt[file] = labels[i]
if ltypes != None:
lab_ltype[file]=ltypes[i]
else:
lab_ltype[file]=None
if ctypes != None:
lab_ctype[file]=ctypes[i]
else:
lab_ctype[file]=None
i = i+1
else:
lab_show = 0
params = {}
params['lab_show'] = lab_show
params['lab_txt'] = lab_txt
params['lab_x'] = lab_x
params['lab_y'] = lab_y
params['lab_dx'] = lab_dx
params['lab_dy'] = lab_dy
params['lab_ex'] = lab_ex
params['lab_oy'] = lab_oy
params['lab_box'] = lab_box
params['lab_ltype'] = lab_ltype
params['lab_dxline'] = lab_dx*fdxline
params['lab_ctype'] = lab_ctype
return params
#################################
def Graph_DrawLabelsForFiles(g,params,file):
#################################
lab_show = params['lab_show']
lab_txt = params['lab_txt']
lab_x = params['lab_x']
lab_y = params['lab_y']
lab_dx = params['lab_dx']
lab_dy = params['lab_dy']
lab_ex = params['lab_ex']
lab_oy = params['lab_oy']
lab_box = params['lab_box']
lab_ltype= params['lab_ltype']
lab_dxline= params['lab_dxline']
lab_ctype= params['lab_ctype']
flab_ex = 1
if lab_show != 0:
g.limits(0,1,0,1)
if lab_ctype[file]!=None:
g.ctype(lab_ctype[file])
if g.ltype(lab_ltype[file])!=None:
g.relocate(lab_x+lab_ex,lab_y+lab_dy-lab_oy[file])
g.ltype(lab_ltype[file])
g.draw(lab_x+lab_ex+lab_dxline,lab_y+lab_dy-lab_oy[file])
g.ltype(0)
flab_ex = 2
g.relocate(lab_x+lab_ex*flab_ex+lab_dxline,lab_y+lab_dy-lab_oy[file])
g.putlabel(6,lab_txt[file])
g.ctype(0)
if lab_box:
g.ctype(0)
g.relocate(lab_x,lab_y)
g.draw(lab_x,lab_y+lab_dy)
g.relocate(lab_x,lab_y)
g.draw(lab_x+lab_dx,lab_y)
g.relocate(lab_x+lab_dx,lab_y+lab_dy)
g.draw(lab_x,lab_y+lab_dy)
g.relocate(lab_x+lab_dx,lab_y+lab_dy)
g.draw(lab_x+lab_dx,lab_y)
#################################
def Graph_SetDefaultsGraphSettings(g):
#################################
g.palette('bgyrw')
g.expand(0.999)
g.setvariable('TeX_strings', '1')
g.ptype(0,0)
g.location(3500, 31000, 3500, 31000)
#################################
def Graph_SetLimits(g,xmin,xmax,ymin,ymax,x,y):
#################################
if xmin == None:
xmin = min(x)
if xmax == None:
xmax = max(x)
if xmin==xmax:
xmin = xmin-0.05*xmin
xmax = xmax+0.05*xmax
else:
xmin = xmin-0.05*(xmax-xmin)
xmax = xmax+0.05*(xmax-xmin)
# cut y values based on x
c = (x>=xmin)*(x<=xmax)
y = compress(c,y)
if ymin == None:
ymin = min(y)
if ymax == None:
ymax = max(y)
if ymin==ymax:
ymin = ymin-0.05*ymin
ymax = ymax+0.05*ymax
else:
ymin = ymin-0.05*(ymax-ymin)
ymax = ymax+0.05*(ymax-ymin)
g.limits(xmin,xmax,ymin,ymax)
return xmin,xmax,ymin,ymax
#################################
def Graph_UseLog(x,y,log):
#################################
if log == 'x':
x = log10(x)
elif log == 'y':
y = log10(y)
elif log == 'xy' or log =='xy':
x = log10(x)
y = log10(y)
else:
pass
c = isinf(x)+isinf(y)
x = compress(logical_not(c),x)
y = compress(logical_not(c),y)
return x,y
#################################
def Graph_DrawBox(g,xmin,xmax,ymin,ymax,log):
#################################
dx= xmax-xmin
dy= ymax-ymin
xs=xb=ys=yb=0
if log == 'x':
if log10(dx)<0:
xb = -1
else:
xb = 0
xs = -1
elif log == 'y':
if log10(dy)<0:
yb = -1
else:
yb = 0
ys = -1
elif log == 'xy' or log =='xy':
if log10(dx)<0:
xb = -1
else:
xb = 0
xs = -1
if log10(dy)<0:
yb = -1
else:
yb = 0
ys = -1
g.ticksize(xs,xb,ys,yb)
g.box()
#################################
def Graph_MakeHistrogram(g,xmin,xmax,ymin,ymax,x,y,nb=50):
#################################
xb = arange(xmin,xmax,(xmax-xmin)/nb)
yb = myNumeric.whistogram(x.astype(float),y.astype(float),xb.astype(float))
h = histogram(x.astype(float),xb)
h2 = where((h>0),h,1)
yb = where((h>0),yb/h2,0)
return xb[1:-2],yb[1:-2]
#################################
def Graph_End(g,ps=None):
#################################
if ps==None:
g.show()
else:
g.write()
g.clean()
####################################################################################################################################
#
# USEFULL FUNCTION
#
####################################################################################################################################
#################################
def histogram(a,bins):
#################################
'''
Return the histogram (n x 1 Float array) of the
n x 1 array "a".
"bins" (m x 1 array) specify the bins of the histogram.
'''
n = searchsorted(sort(a),bins)
n = concatenate([n,[len(a)]])
return n[1:]-n[:-1]
c4science · Help