Homec4science
Diffusion Gtools (master) Browse

gtools/Gtools/gscripts/gatime2Myr
3f5b7a0e0005master

Last Change
Raw File

gatime2Myr
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
# DlogA = da / a
#
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parse = pt.add_units_options(parser)
parser.add_option("--a1",
action="store",
dest="a1",
type="float",
default = 0.1,
help="a1",
metavar=" FLOAT")
parser.add_option("--a2",
action="store",
dest="a2",
type="float",
default = 1.,
help="a2",
metavar=" FLOAT")
parser.add_option("--a",
action="store",
dest="a",
type="float",
default = None,
help="a",
metavar=" FLOAT")
parser.add_option("--da",
action="store",
dest="da",
type="float",
default = None,
help="da",
metavar=" FLOAT")
parser.add_option("--Dloga",
action="store",
dest="Dloga",
type="float",
default = None,
help="Dloga",
metavar=" FLOAT")
(options, args) = parser.parse_args()
#pt.check_files_number(args)
files = args
return files,options
#######################################
# MakePlot
#######################################
def MakePlot(files,opt):
# define local units
unit_params = pt.do_units_options(opt)
system_of_units = units.Set_SystemUnits_From_Params(unit_params)
# convert a to Myrs
out_units = units.UnitSystem('local',[units.Unit_cm,units.Unit_Msol,units.Unit_Myr,units.Unit_K])
hubbleparam = unit_params["HubbleParam"]
omegalambda = unit_params["OmegaLambda"]
omega0 = unit_params["Omega0"]
Hubble = ctes.HUBBLE.into(system_of_units)
pars = {"Hubble":Hubble,"HubbleParam":hubbleparam,"OmegaLambda":omegalambda,"Omega0":omega0}
funit = system_of_units.convertionFactorTo(out_units.UnitTime)
dt = None
if opt.a != None:
a1 = opt.a
if opt.da != None:
da = opt.da
a2 = a1 + da
dt = cosmo.dt_da(da,a1,pars)/ hubbleparam * funit
elif opt.Dloga != None:
da = opt.Dloga*a1
a2 = a1 + da
dt = cosmo.dt_da(da,a1,pars)/ hubbleparam * funit
else:
print "with --a you need to define either --da or --Dloga"
sys.exit()
else:
a1 = opt.a1
a2 = opt.a2
T1 = a1
T2 = a2
T1Myr = cosmo.CosmicTime_a(T1,pars) / hubbleparam * funit
T2Myr = cosmo.CosmicTime_a(T2,pars) / hubbleparam * funit
print "a1 = %8.6f = %10.4f [Myr]"%(T1,T1Myr)
print "a2 = %8.6f = %10.4f [Myr]"%(T2,T2Myr)
if dt==None:
dt = T2Myr-T1Myr
print " Delta t = %g [Myr]"%(dt)
########################################################################
# MAIN
########################################################################
if __name__ == '__main__':
files,opt = parse_options()
MakePlot(files,opt)
c4science · Help