guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Do 3d bubble for galaxies in ksz

Browse Source
This commit is contained in:
Guilhem Lavaux 2014-09-29 10:02:03 +02:00
parent adfc76a96d
commit d58c617e63
1 changed files with 19 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

40
python_sample/ksz/gal_prof.py
View File

@ -1,4 +1,5 @@
import numpy as np
import numexpr as ne
from .constants import *
# -----------------------------------------------------------------------------
@ -6,7 +7,7 @@ from .constants import *
# -----------------------------------------------------------------------------
class KSZ_Profile(object):
R_star= 0.015 # 15 kpc/h
R_star= 0.0 # 15 kpc/h
L_gal0 = 10**(0.4*(tmpp_cat['Msun']-tmpp_cat['Mstar']))
def __init__(self):
@ -58,31 +59,28 @@ class KSZ_Isothermal(KSZ_Profile):
self.rho0 = self.sigma_FP**2/(2*np.pi*G) # * (Lgal/L_gal0)**(2./3)
self.rGalaxy = self.R_gal*(Lgal/self.L_gal0)**(1./3)
self.rInnerGalaxy = self.R_innergal*(Lgal/self.L_gal0)**(1./3)
# self._prepare()
def _prepare(self, x_min, x_max):
from scipy.integrate import quad
from numpy import sqrt, log10
lmin = log10(x_min)
lmax = log10(x_max)
x = 10**(np.arange(100)*(lmax-lmin)/100.+lmin)
profile = np.empty(x.size)
self._prepare()
nu_tilde = lambda u: (1/(u**2*(1+u)))
for i in range(x.size):
profile[i] = 2*quad(lambda y: (nu_tilde(sqrt(x[i]**2+y**2))), 0, args.x)[0]
self._table = x,profile
def _prepare(self):
pass
def evaluate_profile(self,r):
rho0, rGalaxy, rInner = self.rho0, self.rGalaxy, self.rInnerGalaxy
Q=rho0*2/r*np.arctan(np.sqrt((rGalaxy/r)**2 - 1))/Mpc
# Q[r<rInner] = 0
D = {'rho0':rho0, 'rGalaxy':rGalaxy, 'rInner': rInner, 'Mpc':Mpc }
Q = np.zeros(r.size)
cond = r <= rInner
D['r'] = r[cond]
ne.evaluate('rho0*2/(Mpc*r) * arctan(sqrt( (rGalaxy/r)**2 -1 ) - arctan(sqrt( (rInner/r)**2 - 1 ))',
local_dict=D, out=Q[cond])
cond = (r > rInner)*(r <= rGalaxy)
D['r'] = r[cond]
ne.evaluate('rho0*2/(Mpc*r) * arctan(sqrt( (rGalaxy/r)**2 -1 ))',
local_dict=D, out=Q[cond])
return Q,np.where(r<rInner)[0]