Skip to content

Commit

Permalink
Cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
@NolanSmyth
NolanSmyth committed Aug 24, 2023
1 parent b91d840 commit 385548d
Show file tree
Hide file tree
Showing 15 changed files with 380 additions and 8,264 deletions.
15 changes: 0 additions & 15 deletions LensCalcPy/_modidx.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,32 +20,17 @@
'LensCalcPy.ffp.Ffp.differential_rate': ('ffp.html#ffp.differential_rate', 'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.differential_rate_integrand': ( 'ffp.html#ffp.differential_rate_integrand',
'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.differential_rate_integrand_m31': ( 'ffp.html#ffp.differential_rate_integrand_m31',
'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.differential_rate_integrand_m31_mass': ( 'ffp.html#ffp.differential_rate_integrand_m31_mass',
'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.differential_rate_integrand_mw': ( 'ffp.html#ffp.differential_rate_integrand_mw',
'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.differential_rate_integrand_mw_mass': ( 'ffp.html#ffp.differential_rate_integrand_mw_mass',
'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.differential_rate_m31': ('ffp.html#ffp.differential_rate_m31', 'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.differential_rate_m31_mass': ( 'ffp.html#ffp.differential_rate_m31_mass',
'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.differential_rate_m31_monochromatic': ( 'ffp.html#ffp.differential_rate_m31_monochromatic',
'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.differential_rate_mass': ('ffp.html#ffp.differential_rate_mass', 'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.differential_rate_monochromatic': ( 'ffp.html#ffp.differential_rate_monochromatic',
'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.differential_rate_mw': ('ffp.html#ffp.differential_rate_mw', 'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.differential_rate_mw_mass': ( 'ffp.html#ffp.differential_rate_mw_mass',
'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.differential_rate_mw_monochromatic': ( 'ffp.html#ffp.differential_rate_mw_monochromatic',
'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.differential_rate_total': ('ffp.html#ffp.differential_rate_total', 'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.f_m': ('ffp.html#ffp.f_m', 'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.mass_func': ('ffp.html#ffp.mass_func', 'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.pl_norm': ('ffp.html#ffp.pl_norm', 'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.pl_norm_new': ('ffp.html#ffp.pl_norm_new', 'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.sticking_point': ('ffp.html#ffp.sticking_point', 'LensCalcPy/ffp.py'),
'LensCalcPy.ffp.Ffp.umin_upper_bound': ('ffp.html#ffp.umin_upper_bound', 'LensCalcPy/ffp.py')},
'LensCalcPy.galaxy': { 'LensCalcPy.galaxy.Galaxy': ('galaxy.html#galaxy', 'LensCalcPy/galaxy.py'),
Expand Down
146 changes: 31 additions & 115 deletions LensCalcPy/ffp.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,8 +52,7 @@ def __init__(self,
self.m_min = m_min
self.m_max = m_max
self.M_norm = 1 #solar mass
# self.Z = self.pl_norm(self.p)
self.Z = self.pl_norm_new()
self.Z = self.pl_norm()
self.Zprime = self.dN_dM_norm()

# Instantiate or use existing MilkyWayModel and M31Model
Expand Down Expand Up @@ -93,35 +92,29 @@ def dN_dlogM(self, A, log10M, M_norm, p):
def dN_dlogM_wrapper(self, M):
return self.dN_dlogM(1, M, self.M_norm, self.p)

def pl_norm_new(self):
def pl_norm(self):
return 1/abs(nquad(self.dN_dlogM_wrapper,[[np.log10(self.m_min), np.log10(self.m_max)]], opts={'points': [np.log10(self.m_min), np.log10(self.m_min*1e3), np.log10(self.m_max*1e3)]})[0])

def mass_func(self, log10m):
#M_norm = 1 solar mass for now. This is dN/dlogM
m = 10**log10m
return self.Z * (m/1)**-self.p

def pl_norm(self, p):
N_ffp = 1 # Number of FFPs per star
return N_ffp/abs(nquad(self.mass_func,[[self.m_min, self.m_max]], opts={'points': [self.m_min, self.m_min*1e3]})[0])

# def differential_rate_integrand(self, umin, d, mf, t, dist_func, density_func, v_disp_func, finite=False):
def differential_rate_integrand(self, umin, d, t, mf, dist_func, density_func, v_disp_func, finite=False):
r = dist_func(d, self.l, self.b)
# ut = self.umin_upper_bound(d, mf) if (self.ut_interp and finite) else self.u_t

def differential_rate_integrand(self, umin, d, t, mf, model, finite=False, v_disp=None):
r = model.dist_center(d, self.l, self.b)
ut = self.umin_upper_bound(d, mf) if finite else self.u_t
if ut <= umin:
return 0
v_rad = velocity_radial(d, mf, umin, t * htosec, ut)
v_disp = v_disp_func(r)
v_rad = velocity_radial(d, mf, umin, t * htosec, ut)
if v_disp is None:
v_disp = model.velocity_dispersion_stars(r)
return 2 * (1 / (ut**2 - umin**2)**0.5 *
#For FFP number density, use stellar density for 1 solar mass stars
density_func(d) / (1 * v_disp**2) *
model.density_stars(d) / (1 * v_disp**2) *
v_rad**4 * (htosec / kpctokm)**2 *
np.exp(-(v_rad**2 / v_disp**2)) *
1)


def differential_rate(self, t, integrand_func, finite=False):
#rewrite using tplquad?
num = 40 # number of discretization points, empirically, result levels off for >~ 40
Expand All @@ -142,85 +135,41 @@ def differential_rate(self, t, integrand_func, finite=False):
0, self.d_upper_bound(10**mf),
lambda d: 0,
lambda d: self.umin_upper_bound(d, 10**mf),
# args=(mf, t),
args=(10**mf, t),
# epsabs=0,
# epsrel=1e-1,
)
else:
single_result, error = dblquad(integrand_func,
#Without finite size effects, integral blows up at M31 center
0, 10,
0, self.ds,
lambda d: 0,
lambda d: self.u_t,
args=(10**mf, t),
# epsabs=0,
# epsrel=1e-2,
)
# if single_result != 0 and error/abs(single_result) >=1:
# print("Warning: error in differential rate integration is large: {}".format(error/abs(single_result)))

result += single_result * ((10**mf/1) ** -self.p) * dm # multiply by mass function and by dlogm. This is for dN/dlogM

result *= self.Z # normalization
return result


def differential_rate_monochromatic(self, t, integrand_func, finite=False, m=1e-10):
if finite:
result, error = dblquad(integrand_func,
0, self.d_upper_bound(m),
lambda d: 0,
lambda d: self.umin_upper_bound(d, m),
# args=(m, t),
args=(m, t, {'points':[0, self.ds]}),
epsabs=0,
epsrel=1e-1,
)

else:
result, error = dblquad(integrand_func,
#Without finite size effects, integral blows up at M31 center
0, self.ds*0.99,
lambda d: 0,
lambda d: self.u_t,
args=(m, t),
)
return result


def differential_rate_integrand_mw(self, umin, d, mf, t, finite=True, vel_func=None):
if vel_func is None:
vel_func = self.mw_model.velocity_dispersion_stars
return self.differential_rate_integrand(umin, d, t, mf, self.mw_model.dist_center, self.mw_model.density_stars, vel_func, finite=finite)

def differential_rate_mw(self, t, finite=True, v_disp=None):
if v_disp:
vel_func = lambda r: v_disp
f = functools.partial(self.differential_rate_integrand_mw, vel_func=vel_func, finite=finite)
return self.differential_rate(t, f, finite=finite)
# return self.differential_rate(t, self.differential_rate_integrand_mw, finite=finite)
f = functools.partial(self.differential_rate_integrand_mw, finite=finite)
return self.differential_rate(t, f, finite=finite)

def differential_rate_mw_monochromatic(self, t, finite=False, m=1e-10):
return self.differential_rate_monochromatic(t, self.differential_rate_integrand_mw, finite=finite, m=m)

def differential_rate_integrand_m31(self, umin, d, mf, t, finite=True, vel_func=None):
if vel_func is None:
vel_func = self.m31_model.velocity_dispersion_stars
return self.differential_rate_integrand(umin, d, t, mf, self.m31_model.dist_center, self.m31_model.density_stars, vel_func, finite=finite)

def differential_rate_m31(self, t, finite=False, v_disp=None):
if v_disp:
vel_func = lambda r: v_disp
f = functools.partial(self.differential_rate_integrand_m31, vel_func=vel_func, finite=finite)
return self.differential_rate(t, f, finite=finite)
f = functools.partial(self.differential_rate_integrand_m31, finite=finite)
return self.differential_rate(t, self.differential_rate_integrand_m31, finite=finite)

def differential_rate_m31_monochromatic(self, t, finite=False, m=1e-10):
return self.differential_rate_monochromatic(t, self.differential_rate_integrand_m31, finite=finite, m=m)
def integrand_func(umin, d, mf, t):
return self.differential_rate_integrand(umin, d, t, mf, self.mw_model, finite=finite, v_disp=v_disp)
return self.differential_rate(t, integrand_func, finite=finite)

def differential_rate_m31(self, t, finite=True, v_disp=None):
def integrand_func(umin, d, mf, t):
return self.differential_rate_integrand(umin, d, t, mf, self.m31_model, finite=finite, v_disp=v_disp)
return self.differential_rate(t, integrand_func, finite=finite)

def differential_rate_mw_mass(self, m, finite=True, v_disp=None, tcad=0.07, tobs=3, epsabs = 1.49e-08, epsrel = 1.49e-08, efficiency=None, monochromatic=False):
def integrand_func(umin, d, t, mf):
return self.differential_rate_integrand(umin, d, t, mf, self.mw_model, finite=finite)
return self.differential_rate_mass(m, integrand_func, finite=finite, tcad=tcad, tobs=tobs, epsabs = epsabs, epsrel = epsrel, efficiency=efficiency, monochromatic=monochromatic)

def differential_rate_m31_mass(self, m, finite=True, v_disp=None, tcad=0.07, tobs=3, epsabs = 1.49e-08, epsrel = 1.49e-08, efficiency=None, monochromatic=False):
def integrand_func(umin, d, t, mf):
return self.differential_rate_integrand(umin, d, t, mf, self.m31_model, finite=finite)
return self.differential_rate_mass(m, integrand_func, finite=finite, tcad=tcad, tobs=tobs, epsabs = epsabs, epsrel = epsrel, efficiency=efficiency, monochromatic=monochromatic)

def umin_upper_bound(self, d, m):
rho = rho_func(m, d, self.ds)
Expand All @@ -247,18 +196,8 @@ def differential_rate_total(self, t, finite=False):

def compute_differential_rate(self, ts, finite=False):
return [self.differential_rate_total(t, finite=finite) for t in ts]

def differential_rate_integrand_mw_mass(self, umin, d, t, mf, finite=True, vel_func=None):
if vel_func is None:
vel_func = self.mw_model.velocity_dispersion_stars
return self.differential_rate_integrand(umin, d, t, mf, self.mw_model.dist_center, self.mw_model.density_stars, vel_func, finite=finite)

def differential_rate_integrand_m31_mass(self, umin, d, t, mf, finite=True, vel_func=None):
if vel_func is None:
vel_func = self.m31_model.velocity_dispersion_stars
return self.differential_rate_integrand(umin, d, t, mf, self.m31_model.dist_center, self.m31_model.density_stars, vel_func, finite=finite)

def differential_rate_mass(self, m, integrand_func, finite=True, tcad=0.07, tobs=3, d_lower = 0, d_upper=np.inf, epsabs = 1.49e-08, epsrel = 1.49e-08, efficiency=None, monochromatic=False):
def differential_rate_mass(self, m, integrand_func, finite=True, tcad=0.07, tobs=3, epsabs = 1.49e-08, epsrel = 1.49e-08, efficiency=None, monochromatic=False):

if efficiency is None:
def efficiency(t):
Expand All @@ -277,38 +216,15 @@ def integrand(t, d, m, finite):
bounds_t = [tcad, tobs]

if finite:
bounds_d = [d_lower, min(self.d_upper_bound(m), d_upper)]
bounds_d = [0, self.d_upper_bound(m)]
else:
bounds_d = [d_lower, min(self.ds, d_upper)]
bounds_d = [0, self.ds]

opts = {"epsabs": epsabs, "epsrel": epsrel, "points":[point, self.ds]}

result, error = nquad(integrand, [bounds_t, bounds_d], args=(m, finite), opts=opts)

# if result > 0 and error/result > 0.1: #redo with higher precision
# print('Error greater than 10%, recalculating integral with higher precision')
# opts = {"epsabs": 0, "epsrel": 1e-1, "points":[point, self.ds]}
# result, error = nquad(integrand, [bounds_t, bounds_d], args=(m, finite), opts=opts)

if monochromatic:
return result
return result*self.f_m(m)



def differential_rate_mw_mass(self, m, finite=True, v_disp=None, tcad=0.07, tobs=3, d_lower=0, d_upper=np.inf, epsabs = 1.49e-08, epsrel = 1.49e-08, efficiency=None, monochromatic=False):
if v_disp:
vel_func = lambda r: v_disp
f = functools.partial(self.differential_rate_integrand_mw_mass, vel_func=vel_func, finite=finite)
else:
f = functools.partial(self.differential_rate_integrand_mw_mass, finite=finite)
return self.differential_rate_mass(m, f, finite=finite, tcad=tcad, tobs=tobs, d_lower=d_lower, d_upper=d_upper, epsabs = epsabs, epsrel = epsrel, efficiency=efficiency, monochromatic=monochromatic)

def differential_rate_m31_mass(self, m, finite=True, v_disp=None, tcad=0.07, tobs=3, d_lower=0, d_upper=np.inf, epsabs = 1.49e-08, epsrel = 1.49e-08, efficiency=None, monochromatic=False):
if v_disp:
vel_func = lambda r: v_disp
f = functools.partial(self.differential_rate_integrand_m31_mass, vel_func=vel_func, finite=finite)
else:
f = functools.partial(self.differential_rate_integrand_m31_mass, finite=finite)
return self.differential_rate_mass(m, f, finite=finite, tcad=tcad, tobs=tobs, d_lower = self.ds*0.5, d_upper=d_upper, epsabs = epsabs, epsrel = epsrel, efficiency=efficiency, monochromatic=monochromatic)

8 changes: 6 additions & 2 deletions LensCalcPy/pbh.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -136,8 +136,12 @@ def rate_mw(self, finite=False):
return self.rate_total(self.differential_rate_integrand_mw, finite=finite)

def rate_m31(self, finite=False):
return self.rate_total(self.differential_rate_integrand_m31, finite=finite)

result = self.rate_total(self.differential_rate_integrand_m31, finite=finite)
if not np.isnan(result):
return result
else:
return 0

def rate_tot(self, finite=False):
return self.rate_mw(finite=finite) + self.rate_m31(finite=finite)

Expand Down
6 changes: 3 additions & 3 deletions LensCalcPy/utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,7 +124,7 @@ def scientific_format(x, pos):
b = int(b)
return r'${} \times 10^{{{}}}$'.format(a, b)

# %% ../nbs/04_utils.ipynb 9
# %% ../nbs/04_utils.ipynb 8
# Add finite size calculation following https://arxiv.org/pdf/1905.06066.pdf

# Compute 'w' parameter given the mass of the primordial black hole and the wavelength
Expand Down Expand Up @@ -239,7 +239,7 @@ def calc_ut_arr(m):
ut_interp = interp2d(d_arr, m_arr, ut_values)
return ut_interp

# %% ../nbs/04_utils.ipynb 24
# %% ../nbs/04_utils.ipynb 13
#Ratio of angular extent of source and lens in plane of lens
# rho == theta_s/theta_l
rho_min = 0.1
Expand All @@ -255,7 +255,7 @@ def calc_ut_arr(m):
# with open('../LensCalcPy/interpolations/ut_interp_rho.pkl', 'wb') as f:
# pickle.dump(ut_interp_rho, f)

# %% ../nbs/04_utils.ipynb 26
# %% ../nbs/04_utils.ipynb 15
def ut_func_new(rho, A_thresh):
# rho = np.asarray(rho)
rho = np.atleast_1d(rho) # Ensure rho is at least 1-dimensional
Expand Down
Loading

0 comments on commit 385548d

Please sign in to comment.