Homec4science
Diffusion Gtools (master) Browse

gtools/Gtools/gscripts/ggetZiforCosmoSims
6c451efb15ebmaster

Last Change
Raw File

ggetZiforCosmoSims
View Options

#!/usr/bin/env python
from optparse import OptionParser
import Ptools as pt
import sys
from numpy import *
import scipy.integrate as si
from scipy import optimize
def parse_options():
usage = "usage: %prog [options] file"
parser = OptionParser(usage=usage)
parser.add_option("--N1",
action="store",
dest="N1",
type="float",
default = 512,
help="number of particles for the higest level along one dimension",
metavar=" FLOAT")
parser.add_option("--N2",
action="store",
dest="N2",
type="float",
default = None,
help="number of particles for the higest level along one dimension",
metavar=" FLOAT")
parser.add_option("-L",
action="store",
dest="L",
type="float",
default = 100,
help="box size in Mpc/h",
metavar=" FLOAT")
parser.add_option("--param",
action="store",
dest="param",
type="string",
default = None,
help="Gadget parameter file",
metavar=" FILE")
parser.add_option("--Omega_b",
action="store",
dest="Omega_b",
type="float",
default = None,
help="Omega_b",
metavar=" FLOAT")
parser.add_option("--Omega_m",
action="store",
dest="Omega_m",
type="float",
default = None,
help="Omega_m",
metavar=" FLOAT")
parser.add_option("--Omega_de",
action="store",
dest="Omega_de",
type="float",
default = None,
help="Omega_de",
metavar=" FLOAT")
parser.add_option("--sigma8",
action="store",
dest="sigma8",
type="float",
default = None,
help="sigma8",
metavar=" FLOAT")
parser.add_option("--h",
action="store",
dest="h",
type="float",
default = None,
help="h",
metavar=" FLOAT")
parser.add_option("--n",
action="store",
dest="n",
type="float",
default = None,
help="n",
metavar=" FLOAT")
(options, args) = parser.parse_args()
#pt.check_files_number(args)
files = args
return files,options
########################################################################
# MAIN
########################################################################
try:
from cosmicpy import *
except ImportError:
print "%s requires cosmicpy"%sys.argv[0]
sys.exit()
files,opt = parse_options()
if opt.param != None:
print "You want to use a parameter file to change the cosmology, this is not implemented"
cosmo = cosmology()
if opt.h!=None:
cosmo.h = opt.h
if opt.Omega_b!=None:
cosmo.Omega_b = opt.Omega_b
if opt.Omega_m!=None:
cosmo.Omega_m = opt.Omega_m
if opt.Omega_de!=None:
cosmo.Omega_de = opt.Omega_de
if opt.sigma8!=None:
cosmo.sigma8 = opt.sigma8
if opt.n!=None:
cosmo.n = opt.n
def _sigmasq_integrand_log_l(logk, z):
"""Integrand used internally by the sigma_r function.
"""
k = exp(logk)
# The 1e-10 factor in the integrand is added to avoid roundoff
# error warnings. It is divided out later.
#print k,power_spectrum(k, z, **cosmology),z
a = z2a(z)
return (k *
(1.e-10 / (2. * math.pi**2.)) * k**2. *
cosmo.pk_lin(k,a,type='eisenhu')
)
def compute_Sigma(kmin,kmax,z):
# Integrate over logk from -infinity to infinity.
integral, error = si.quad(_sigmasq_integrand_log_l,
log(kmin),
log(kmax),
args=(z),
limit=10000)#, epsabs=1e-9, epsrel=1e-9)
sigma = sqrt(1.e10 * integral)
return sigma
def getZ(z,kmin,kmax,sigma):
return compute_Sigma(kmin,kmax,z) - sigma
#######################
# first level
#######################
kmin = 2*pi/opt.L
kmax = pi*opt.N1/opt.L
zi1 = optimize.bisect(getZ, a=5, b=500, args = (kmin,kmax,0.1), xtol = 1e-3, maxiter = 500)
zi2 = optimize.bisect(getZ, a=5, b=500, args = (kmin,kmax,0.2), xtol = 1e-3, maxiter = 500)
print "N=%5d zinit(max)=%4.1f zinit(min)=%4.1f"%(opt.N1,zi1,zi2)
#######################
# second level
#######################
if opt.N2!=None:
kmin = 2*pi/opt.L
kmax = pi*opt.N2/opt.L
zi1b = optimize.bisect(getZ, a=5, b=500, args = (kmin,kmax,0.1), xtol = 1e-3, maxiter = 500)
zi2b = optimize.bisect(getZ, a=5, b=500, args = (kmin,kmax,0.2), xtol = 1e-3, maxiter = 500)
print "N=%5d zinit(max)=%4.1f zinit(min)=%4.1f"%(opt.N2,zi1b,zi2b)
# find the intersection
print
if opt.N1 < opt.N2:
print "zinit should be in [%4.1f,%4.1f]"%(zi2b,zi1)
else:
print "zinit should be in [%4.1f,%4.1f]"%(zi2,zi1b)
print
c4science · Help