Page MenuHomec4science
Files F105541473

gcosmo2giso
No OneTemporary
Actions

File Metadata

Created
Mon, Mar 17, 23:38

gcosmo2giso
View Options

#!/usr/bin/env python
import Ptools as pt
from optparse import OptionParser
from pNbody import *
import sys
from scipy import optimize
from numpy import *
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("--a0",
action="store",
dest="a0",
type="float",
default = 1.,
help="initial scaling factor",
metavar=" FLOAT")
parser.add_option("--Omega0",
action="store",
dest="Omega0",
type="float",
default = 0.24,
help="Omega0",
metavar=" FLOAT")
parser.add_option("--OmegaLambda",
action="store",
dest="OmegaLambda",
type="float",
default = 0.76,
help="OmegaLambda",
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()
# read file
file = files[0]
nb = Nbody(file,ftype=opt.ftype)
a = nb.atime
# unit options
unit_params1 = io.read_params(opt.GadgetParameterFile1)
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'])
# set some parameters
Hubble = ctes.HUBBLE.into(nb.localsystem_of_units)
cosmo.Omega0 = unit_params1['Omega0']
cosmo.OmegaLambda= unit_params1['OmegaLambda']
cosmo.HubbleParam= unit_params1["HubbleParam"]
cosmo.Hubble = Hubble
cosmo.setdefault() # this set defaultpars
# compute some factors
Ha = cosmo.Hubble_a(a,pars=cosmo.defaultpars)
pos = nb.pos
vel = nb.vel
for ar in nb.get_list_of_array():
if ar == "pos":
#print "pos"
nb.pos = pos * a
elif ar == "vel":
#print "vel",
nb.vel = vel*sqrt(a) + pos*a*Ha
elif ar == "mass":
#print "mass"
pass
elif ar == "u":
#print "u"
pass
elif ar == "rho":
#print "rho"
nb.rho = nb.rho / a**3
elif ar == "rsp":
#print "rsp"
nb.rsp = nb.rsp * a
elif ar == "tpe":
pass
elif ar == "num":
pass
elif ar == "minit":
#print "minit"
pass
else:
print "skipping %s"%ar
# compute the time
hubbleparam = unit_params1["HubbleParam"]
omegalambda = unit_params1["OmegaLambda"]
omega0 = unit_params1["Omega0"]
units = local_units_2.UnitTime
funit = nb.localsystem_of_units.convertionFactorTo(units)
Hubble = ctes.HUBBLE.into(nb.localsystem_of_units)
pars = {"Hubble":Hubble,"HubbleParam":hubbleparam,"OmegaLambda":omegalambda,"Omega0":omega0}
age = cosmo.CosmicTime_a(nb.atime,pars) / hubbleparam * funit
print "--------------------------------------"
print "a= %g -> t=%g "%(nb.atime,age)
print "--------------------------------------"
nb.atime = age
if opt.outputfile!=None:
# kill comobile flag
nb.comovingintegration= False
nb.omega0=0
nb.omegalambda=0
nb.rename(opt.outputfile)
nb.nzero=None
nb.massarr=None
nb.write()

Event Timeline

c4science · Help