Page MenuHomec4science
Files F110311443

gchangeunits
No OneTemporary
Actions

File Metadata

Created
Fri, Apr 25, 16:52

gchangeunits
View Options

#!/usr/bin/env python
from optparse import OptionParser
from pNbody import *
import sys
from scipy import optimize
from numpy import *
import Ptools as pt
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser = pt.add_ftype_options(parser)
parser.add_option("-o",
action="store",
dest="outputfile",
type="string",
default = None,
help="output file name",
metavar=" STRING")
parser.add_option("--UnitLength_in_cm_in",
action="store",
dest="UnitLength_in_cm_in",
type="float",
default = 3.085e+21,
help="UnitLength_in_cm_in",
metavar=" FLOAT")
parser.add_option("--UnitMass_in_g_in",
action="store",
dest="UnitMass_in_g_in",
type="float",
default = 1.989e+43,
help="UnitMass_in_g_in",
metavar=" FLOAT")
parser.add_option("--UnitVelocity_in_cm_per_s_in",
action="store",
dest="UnitVelocity_in_cm_per_s_in",
type="float",
default = 20725573.785998672,
help="UnitVelocity_in_cm_per_s_in",
metavar=" FLOAT")
parser.add_option("--UnitLength_in_cm_out",
action="store",
dest="UnitLength_in_cm_out",
type="float",
default = 3.085678e+21,
help="UnitLength_in_cm_out",
metavar=" FLOAT")
parser.add_option("--UnitMass_in_g_out",
action="store",
dest="UnitMass_in_g_out",
type="float",
default = 1.989e+43,
help="UnitMass_in_g_out",
metavar=" FLOAT")
parser.add_option("--UnitVelocity_in_cm_per_s_out",
action="store",
dest="UnitVelocity_in_cm_per_s_out",
type="float",
default = 100000.0,
help="UnitVelocity_in_cm_per_s_out",
metavar=" FLOAT")
parser.add_option("--param1",
action="store",
dest="GadgetParameterFile1",
type="string",
default = None,
help="Gadget parameter file 1",
metavar=" FILE")
parser.add_option("--param2",
action="store",
dest="GadgetParameterFile2",
type="string",
default = None,
help="Gadget parameter file 2",
metavar=" FILE")
(options, args) = parser.parse_args()
files = args
return files,options
###################################################
# main
###################################################
files,opt = parse_options()
file = files[0]
# define local units
cosmo_corrections = []
if opt.GadgetParameterFile1!=None:
unit_params1 = io.read_params(opt.GadgetParameterFile1)
local_units_1 = units.Set_SystemUnits_From_Params(unit_params1)
if unit_params1.has_key('HubbleParam'):
local_units_1.CorrectFromHubbleParameter(unit_params1['HubbleParam'])
ccorrect=False
if unit_params1.has_key('Omega0'):
if unit_params1['Omega0'] !=0:
ccorrect=unit_params1
cosmo_corrections.append(ccorrect)
if opt.GadgetParameterFile2!=None:
unit_params2 = io.read_params(opt.GadgetParameterFile2)
local_units_2 = units.Set_SystemUnits_From_Params(unit_params2)
if unit_params2.has_key('HubbleParam'):
local_units_2.CorrectFromHubbleParameter(unit_params2['HubbleParam'])
ccorrect=False
if unit_params2.has_key('Omega0'):
if unit_params2['Omega0'] !=0:
ccorrect=unit_params2
cosmo_corrections.append(ccorrect)
UnitLengthRatio = local_units_1.convertionFactorTo(local_units_2.UnitLength)
UnitVelocityRatio = local_units_1.convertionFactorTo(local_units_2.UnitVelocity)
UnitMassRatio = local_units_1.convertionFactorTo(local_units_2.UnitMass)
UnitEnergySpecRatio = local_units_1.convertionFactorTo(local_units_2.UnitEnergy)/local_units_1.convertionFactorTo(local_units_2.UnitMass)
UnitDensityRatio = local_units_1.convertionFactorTo(local_units_2.UnitDensity)
UnitTimeRatio = local_units_1.convertionFactorTo(local_units_2.UnitTime)
print local_units_1.convertionFactorTo(local_units_2.UnitEnergy)
print local_units_1.convertionFactorTo(local_units_2.UnitMass)
print "UnitTimeRatio=%g"%UnitTimeRatio
nb = Nbody(file,ftype=opt.ftype)
# check
if unit_params1['HubbleParam'] != nb.hubbleparam:
print "unit_params1['HubbleParam'] = %f"%(unit_params1['HubbleParam'])
print "is different from"
print "nb.hubbleparam = %f"%(nb.hubbleparam)
sys.exit()
for a in nb.get_list_of_array():
if a == "pos":
#print "pos",UnitLengthRatio
nb.pos = nb.pos * UnitLengthRatio
elif a == "vel":
#print "vel",
nb.vel = nb.vel * UnitVelocityRatio
elif a == "mass":
#print "mass",UnitMassRatio
nb.mass = nb.mass * UnitMassRatio
elif a == "u":
#print "u",UnitEnergySpecRatio
nb.u = nb.u* UnitEnergySpecRatio
elif a == "rho":
#print "rho",UnitDensityRatio
nb.rho = nb.rho * UnitDensityRatio
elif a == "rsp":
#print "rsp",UnitLengthRatio
nb.rsp = nb.rsp * UnitLengthRatio
elif a == "tpe":
pass
elif a == "num":
pass
elif a == "minit":
#print "minit",UnitMassRatio
nb.minit = nb.minit * UnitMassRatio
elif a == "tstar":
tpe = 1
idx = compress(nb.tpe==tpe,arange(nb.nbody))
nbs = nb.select(tpe)
ct = nbs.CosmicTime(local_units_2.UnitTime)
put(nb.tstar,idx,ct)
elif a== "thtsnii":
tpe = 0
idx = compress(nb.tpe==tpe,arange(nb.nbody))
nbg = nb.select(tpe)
ct = nbg.CosmicTime(local_units_2.UnitTime,age=nbg.thtsnii)
print ct
put(nb.thtsnii,idx,ct)
elif a== "thtsnia":
tpe = 0
idx = compress(nb.tpe==tpe,arange(nb.nbody))
nbg = nb.select(tpe)
ct = nbg.CosmicTime(local_units_2.UnitTime,age=nbg.thtsnia)
put(nb.thtsnia,idx,ct)
else:
print "skipping %s"%a
# print some info concernint the units
units = local_units_2.UnitTime
funit = nb.localsystem_of_units.convertionFactorTo(units)
Hubble = ctes.HUBBLE.into(nb.localsystem_of_units)
pars = {"Hubble":Hubble,"HubbleParam":nb.hubbleparam,"OmegaLambda":nb.omegalambda,"Omega0":nb.omega0}
age = cosmo.CosmicTime_a(nb.atime,pars) / nb.hubbleparam * funit
print "--------------------------------------"
print "a= %g -> t=%g "%(nb.atime,age)
print "--------------------------------------"
if opt.outputfile!=None:
if unit_params2.has_key('HubbleParam'):
nb.hubbleparam = unit_params2['HubbleParam']
nb.rename(opt.outputfile)
nb.nzero=None
nb.massarr=None
nb.write()

Event Timeline

c4science · Help