From d58c617e633b928009ae047142ec00ff87ca9502 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@iap.fr>
Date: Mon, 29 Sep 2014 10:02:03 +0200
Subject: [PATCH] Do 3d bubble for galaxies in ksz

---
 python_sample/ksz/gal_prof.py | 40 +++++++++++++++++------------------
 1 file changed, 19 insertions(+), 21 deletions(-)

diff --git a/python_sample/ksz/gal_prof.py b/python_sample/ksz/gal_prof.py
index ee7c0ab..eb44610 100644
--- a/python_sample/ksz/gal_prof.py
+++ b/python_sample/ksz/gal_prof.py
@@ -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]