Skip to content

Commit aa7da86

Browse files
martinkilbingermartinkilbinger
authored
Merge pull request #9 from martinkilbinger/cosmo
added cosmo (for critical surface mass density)
2 parents 671181f + 84d22cd commit aa7da86

12 files changed

Lines changed: 510 additions & 28 deletions

File tree

‎cs_util/__init__.py‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,3 @@
1-
# -*- coding: utf-8 -*-
2-
31
"""cs_util PACKAGE.
42
53
Provide a basic description of what your package contains.

‎cs_util/cosmo.py‎

Lines changed: 223 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,223 @@
1+
"""COSMO.
2+
3+
:Name: cosmo.py
4+
5+
:Description: This file contains methods for cosmological
6+
quantities.
7+
8+
:Author: Martin Kilbinger <martin.kilbinger@cea.fr>
9+
10+
"""
11+
12+
import numpy as np
13+
14+
from astropy import constants
15+
from astropy import units
16+
17+
18+
19+
def sigma_crit(z_lens, z_source, cosmo, d_lens=None, d_source=None):
20+
"""Critical surface mass density.
21+
22+
Parameters
23+
----------
24+
z_lens : float
25+
lens redshift
26+
z_source : float
27+
source redshift
28+
cosmo : pyccl.core.Cosmology
29+
cosmological parameters
30+
d_lens : astropy.units.Quantity, optional
31+
precomputed anguar diameter distance to lens, computed from z_lens
32+
if ``None`` (default)
33+
d_source : astropy.units.Quantity, optional
34+
precomputed anguar diameter distance to sourcce, computed from z_source
35+
if ``None`` (default)
36+
37+
Returns
38+
-------
39+
astropy.units.Quantity
40+
critical surface mass density with units of M_sol / pc^2
41+
42+
"""
43+
unit_return = units.Msun / units.pc**2
44+
45+
# Return 0 if lens behind source
46+
if z_lens >= z_source:
47+
return 0.0 * unit_return
48+
49+
a_lens = 1 / (1 + z_lens)
50+
a_source = 1 / (1 + z_source)
51+
if not d_lens:
52+
d_lens = cosmo.angular_diameter_distance(a_lens) * units.Mpc
53+
if not d_source:
54+
d_source = cosmo.angular_diameter_distance(a_source) * units.Mpc
55+
56+
d_lens_source = cosmo.angular_diameter_distance(
57+
a_lens,
58+
a_source
59+
) * units.Mpc
60+
61+
frac = d_source / (d_lens_source * d_lens)
62+
pref = constants.c**2 / (4 * np.pi * constants.G)
63+
64+
sigma_cr = (pref * frac).to(unit_return)
65+
66+
return sigma_cr
67+
68+
69+
def sigma_crit_eff(
70+
z_lens,
71+
z_source_arr,
72+
nz_source_arr,
73+
cosmo,
74+
d_lens=None,
75+
d_source_arr=None,
76+
):
77+
"""Effective critical surface mass density, which
78+
is sigma_crit(z_lens, z_source) weighted by nz_source.
79+
80+
Parameters
81+
----------
82+
z_lens : float
83+
lens redshift
84+
z_source_arr : list
85+
source redshifts
86+
nz_source_arr : list
87+
number of galaxies at z_source
88+
cosmo : pyccl.core.Cosmology
89+
cosmological parameters
90+
d_lens : astropy.units.Quantity, optional
91+
precomputed anguar diameter distance to lens;
92+
computed from z_lens if ``None`` (default)
93+
d_source_arr : list, optional
94+
precompuated angular diameter distances to sources;
95+
computed from z_source_arr if ``None`` (default);
96+
needs to be list of astropy.units.Quantity
97+
98+
Raises
99+
------
100+
IndexError
101+
If lists ``z_source_arr``, ``nz_source_arr``, and d_source_arr
102+
do not match
103+
104+
Returns
105+
-------
106+
astropy.units.Quantity
107+
effective critical surface mass density with units of M_sol / pc^2
108+
109+
"""
110+
n_source = len(z_source_arr)
111+
112+
if d_source_arr is None:
113+
d_source_arr = [None] * n_source
114+
115+
if (len(nz_source_arr) != n_source) or (len(d_source_arr) != n_source):
116+
raise IndexError(
117+
'Lists for source z, n(z), and/or d_ang have different lenghts'
118+
)
119+
120+
sigma_cr_arr = []
121+
for idx in range(n_source):
122+
sigma_cr = sigma_crit(
123+
z_lens,
124+
z_source_arr[idx],
125+
cosmo,
126+
d_lens=d_lens,
127+
d_source=d_source_arr[idx]
128+
)
129+
130+
# Get unit
131+
if len(sigma_cr_arr) == 0:
132+
unit = sigma_cr.unit
133+
134+
sigma_cr_arr.append(sigma_cr.value)
135+
136+
# Mean sigma_cr weighted by source redshifts.
137+
# np.average can only deal with unitless quantities.
138+
sigma_cr_eff = np.average(sigma_cr_arr, weights=nz_source_arr)
139+
140+
return sigma_cr_eff * unit
141+
142+
143+
def sigma_crit_m1_eff(
144+
z_lens,
145+
z_source_arr,
146+
nz_source_arr,
147+
cosmo,
148+
d_lens=None,
149+
d_source_arr=None,
150+
):
151+
"""Effective inverse critical surface mass density, which
152+
is sigma_crit^{-1}(z_lens, z_source) weighted by nz_source.
153+
See Eq. (17) in :cite:`2004AJ....127.2544S`.
154+
155+
Parameters
156+
----------
157+
z_lens : float
158+
lens redshift
159+
z_source_arr : list
160+
source redshifts
161+
nz_source_arr : list
162+
number of galaxies at z_source
163+
cosmo : pyccl.core.Cosmology
164+
cosmological parameters
165+
d_lens : astropy.units.Quantity, optional
166+
precomputed anguar diameter distance to lens;
167+
computed from z_lens if ``None`` (default)
168+
d_source_arr : float, optional
169+
precomputed anguar diameter distance to sources;
170+
computed from z_source_arr if ``None`` (default);
171+
needs to be list of astropy.units.Quantity
172+
173+
Raises
174+
------
175+
IndexError
176+
If lists ``z_source_arr``, ``nz_source_arr``, and ``d_source_arr``
177+
do not match
178+
179+
Returns
180+
-------
181+
astropy.units.Quantity
182+
effective inverse critical surface mass density with units of
183+
M_sol / pc^2
184+
185+
"""
186+
n_source = len(z_source_arr)
187+
188+
if d_source_arr is None:
189+
d_source_arr = [None] * n_source
190+
191+
if (len(nz_source_arr) != n_source) or (len(d_source_arr) != n_source):
192+
raise IndexError(
193+
'Lists for source z, n(z), and/or d_ang have different lenghts'
194+
)
195+
196+
sigma_cr_m1_arr = []
197+
weights = []
198+
199+
for idx in range(n_source):
200+
sigma_cr = sigma_crit(
201+
z_lens,
202+
z_source_arr[idx],
203+
cosmo,
204+
d_lens=d_lens,
205+
d_source=d_source_arr[idx]
206+
)
207+
208+
# Get unit
209+
if len(sigma_cr_m1_arr) == 0:
210+
unit = 1 / sigma_cr.unit
211+
212+
# If lens behind source: continue
213+
if sigma_cr == 0:
214+
continue
215+
216+
sigma_cr_m1 = 1 / sigma_cr
217+
218+
sigma_cr_m1_arr.append(sigma_cr_m1.value)
219+
weights.append(nz_source_arr[idx])
220+
221+
sigma_cr_m1_eff = np.average(sigma_cr_m1_arr, weights=weights)
222+
223+
return sigma_cr_m1_eff * unit

‎cs_util/logging.py‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
"""LOGGING.
22
3-
:Description: This script contains utility methods for job execution and progress logging.
3+
:Description: This script contains utility methods for job execution and
4+
progress logging.
45
56
:Author: Martin Kilbinger <martin.kilblinger@cea.fr>
67

‎cs_util/tests/__init__.py‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,3 @@
1-
# -*- coding: utf-8 -*-
2-
31
"""UNIT TESTS.
42
53
Unit testing framework for the package.

‎cs_util/tests/test_calc.py‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,3 @@
1-
# -*- coding: utf-8 -*-
2-
31
"""UNIT TESTS FOR CALC SUBPACKAGE.
42
53
This module contains unit tests for the calc subpackage.

‎cs_util/tests/test_cat.py‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,3 @@
1-
# -*- coding: utf-8 -*-
2-
31
"""UNIT TESTS FOR CAT SUBPACKAGE.
42
53
This module contains unit tests for the cat subpackage.

0 commit comments

Comments
 (0)