Skip to content

Commit

Permalink
Weak lensing profile (#824)
Browse files Browse the repository at this point in the history 
* Created weak_lensing_profile.py file: WL calculator.
Implemented function: Sigma_{crit}.
Work in progress.

* Created module to compute  weak lensing functions: convergence, shear, reduced shear...

Issue #791.


Co-authored-by: Sandro Dias Pinto Vitenti <sandro@isoftware.com.br>
Co-authored-by: David Alonso <dam.phys@gmail.com>
Co-authored-by: Nick Koukoufilippas <nikfilippas@gmail.com>
  • Loading branch information
@pennalima @damonge @nikfilippas
4 people committed Jul 14, 2021
1 parent e287373 commit c3d04fc
Show file tree
Hide file tree
Showing 4 changed files with 298 additions and 1 deletion.
101 changes: 101 additions & 0 deletions benchmarks/data/haloprofile_nfw_wl_numcosmo.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,101 @@
# r3d convergence shear reduced_shear magnification
1.000000000000000021e-02 1.984125132863828966e+00 3.149527996098407390e-01 -3.200332854962457674e-01 1.150341579178808304e+00
1.097498765493056146e-02 1.925797696855528818e+00 3.147006357739102200e-01 -3.399237617924414923e-01 1.319148293675306283e+00
1.204503540258782326e-02 1.867561401806574262e+00 3.144079335256003693e-01 -3.624042435162395659e-01 1.529495313993854300e+00
1.321941148466028795e-02 1.809430422725087473e+00 3.140686143452386392e-01 -3.880118729511952691e-01 1.796820837296375872e+00
1.450828778495939428e-02 1.751420947118417892e+00 3.136757921096686164e-01 -4.174435026233520585e-01 2.144813808024266688e+00
1.592282793341092198e-02 1.693551415260160953e+00 3.132216907121352922e-01 -4.516200007963957486e-01 2.611608812713636762e+00
1.747528400007683849e-02 1.635842778190667701e+00 3.126975595084897641e-01 -4.917844005373314675e-01 3.262469506632767136e+00
1.917910261672488639e-02 1.578318771921248587e+00 3.120935883473702765e-01 -5.396566798455381120e-01 4.218515846095773014e+00
2.104904144512020903e-02 1.521006205359987984e+00 3.113988245699972346e-01 -5.976873621972249584e-01 5.731373786232710543e+00
2.310129700083160542e-02 1.463935258323220356e+00 3.106010949400827559e-01 -6.694923254220292996e-01 8.420139091693549460e+00
2.535364493970111363e-02 1.407139784630655699e+00 3.096869362835614203e-01 -7.606403205338911899e-01 1.431499750998121812e+01
2.782559402207124277e-02 1.350657613702142879e+00 3.086415393711400834e-01 -8.801792041889252571e-01 3.609956838596566087e+01
3.053855508833415444e-02 1.294530842292835437e+00 3.074487114609367189e-01 -1.043859139054976071e+00 -1.285959756043685900e+02
3.351602650938842465e-02 1.238806106052092959e+00 3.060908637714297997e-01 -1.281754762605020082e+00 -2.727526094865049089e+01
3.678379771828634015e-02 1.183534818525016297e+00 3.045490309780096161e-01 -1.659352886964593710e+00 -1.693047660289394685e+01
4.037017258596553582e-02 1.128773363118800166e+00 3.028029305394867032e-01 -2.351440726605349418e+00 -1.331433188322322714e+01
4.430621457583882455e-02 1.074583221547640521e+00 3.008310701581077273e-01 -4.033495254236904692e+00 -1.177347700692085297e+01
4.862601580065353118e-02 1.021031020506283360e+00 2.986109118795262085e-01 -1.419859353902052668e+01 -1.127063151326649404e+01
5.336699231206309957e-02 9.681884769973057026e-01 2.961191010910918520e-01 9.308548385627849697e+00 -1.153741220612939244e+01
5.857020818056667133e-02 9.161322220798319149e-01 2.933317678643725945e-01 3.497550252775142976e+00 -1.265667029820050971e+01
6.428073117284321958e-02 8.649434830717620537e-01 2.902249065693771057e-01 2.148914492764445683e+00 -1.515375753374188861e+01
7.054802310718645553e-02 8.147076056875952599e-01 2.867748373365682779e-01 1.547688119637890658e+00 -2.087397847477722479e+01
7.742636826811269413e-02 7.655134004797721214e-01 2.829587496837862592e-01 1.206716077860036984e+00 -3.986972381463943549e+01
8.497534359086446332e-02 7.174523024776704450e-01 2.787553244145319908e-01 9.865779366066238509e-01 4.697767957553550673e+02
9.326033468832199691e-02 6.706173392066030470e-01 2.741454247932166854e-01 8.323007171442238228e-01 2.999649944185494732e+01
1.023531021899026366e-01 6.251019103159118906e-01 2.691128421630542555e-01 7.178293236698671809e-01 1.467849921725975193e+01
1.123324032978027659e-01 5.809983916821359617e-01 2.636450748770182706e-01 6.292221071309378466e-01 9.429192314538585862e+00
1.232846739442066547e-01 5.383965877143392920e-01 2.577341130755629317e-01 5.583453376121618739e-01 6.818917772001136690e+00
1.353047774579807516e-01 4.973820670947949885e-01 2.513771960049159815e-01 5.001357483445115770e-01 5.278875278932584791e+00
1.484968262254464932e-01 4.580344285546340544e-01 2.445775038613621777e-01 4.512786729406065844e-01 4.275175515670412985e+00
1.629750834620644351e-01 4.204255535601875704e-01 2.373447432471572638e-01 4.095155414547852479e-01 3.576870373229005917e+00
1.788649529057435017e-01 3.846179109281603337e-01 2.296955848781215004e-01 3.732568577427461398e-01 3.068095892801021130e+00
1.963040650040271395e-01 3.506629833657262663e-01 2.216539147145778210e-01 3.413541951812367747e-01 2.684505360393283357e+00
2.154434690031884481e-01 3.185998868548529095e-01 2.132508654998009623e-01 3.129598328293435872e-01 2.387602521606677186e+00
2.364489412645408295e-01 2.884542500535876575e-01 2.045246047830871772e-01 2.874370408346761940e-01 2.153005301179315989e+00
2.595024211399737379e-01 2.602374122509797938e-01 1.955198675150587040e-01 2.643008321223631607e-01 1.964556772589045375e+00
2.848035868435802032e-01 2.339459851450149952e-01 1.862872354928013230e-01 2.431776766133987300e-01 1.811148569898278549e+00
3.125715849688237014e-01 2.095618065754275283e-01 1.768821812304683661e-01 2.237773714654737178e-01 1.684905114189628073e+00
3.430469286314919430e-01 1.870522945539992987e-01 1.673639089312765060e-01 2.058729089338619811e-01 1.580095115140230488e+00
3.764935806792469308e-01 1.663711890543768668e-01 1.577940387454112581e-01 1.892857308595397303e-01 1.492452587177848455e+00
4.132012400115338546e-01 1.474596487803882250e-01 1.482351911069283157e-01 1.738746921419833391e-01 1.418738801393121918e+00
4.534878508128584174e-01 1.302476522189678632e-01 1.387495346012302955e-01 1.595276344527463419e-01 1.356451401269003210e+00
4.977023564332111460e-01 1.146556381522732121e-01 1.293973628743185222e-01 1.461548392359604542e-01 1.303626362637199554e+00
5.462277217684342601e-01 1.005963114929878804e-01 1.202357633947797322e-01 1.336838673570131530e-01 1.258700204885940144e+00
5.994842503189411476e-01 8.797653633564023190e-02 1.113174337622037241e-01 1.220554494452901900e-01 1.220412283091920846e+00
6.579332246575682053e-01 7.669923944994185006e-02 1.026896905110780828e-01 1.112201948690050568e-01 1.187734331372057195e+00
7.220809018385467848e-01 6.666525353574856427e-02 9.439370210592364085e-02 1.011359562724038069e-01 1.159818893776534354e+00
7.924828983539177196e-01 5.777503946934290507e-02 8.646396335339728845e-02 9.176573215031484654e-02 1.135961068170265875e+00
8.697490026177834288e-01 4.993023930321049270e-02 7.892801406887779703e-02 8.307601960827729060e-02 1.115569772301444829e+00
9.545484566618341882e-01 4.303502565517343392e-02 7.180639168389262028e-02 7.503554843588074175e-02 1.098145907295957180e+00
1.047615752789665233e+00 3.699722712008181313e-02 6.511279646577365632e-02 6.761433954239808031e-02 1.083265571818351036e+00
1.149756995397736903e+00 3.172922242260244052e-02 5.885443972123626660e-02 6.078303826176536390e-02 1.070567009340558595e+00
1.261856883066021062e+00 2.714860742441284336e-02 5.303254001662355971e-02 5.451247787827279145e-02 1.059740337505539198e+00
1.384886371393873050e+00 2.317864831178172505e-02 4.764293001108226011e-02 4.877343224402504174e-02 1.050519366741769378e+00
1.519911082952934755e+00 1.974854110600989324e-02 4.267673666379893233e-02 4.353651940691910527e-02 1.042674999849535045e+00
1.668100537200059241e+00 1.679350208382209345e-02 3.812109987577860792e-02 3.877222125420551335e-02 1.036009837883888007e+00
1.830738280295369780e+00 1.425471597393254165e-02 3.395989850763583839e-02 3.445098755018522613e-02 1.030353715298091943e+00
2.009233002565047776e+00 1.207916923593401089e-02 3.017445760998064067e-02 3.054339646492065449e-02 1.025559959107002594e+00
2.205130739903045534e+00 1.021939471567260595e-02 2.674421608614942561e-02 2.702034768449197363e-02 1.021502219761853070e+00
2.420128264794383366e+00 8.633151906217483365e-03 2.364733943160232318e-02 2.385326832047273207e-02 1.018071760424810623e+00
2.656087782946686904e+00 7.283064298801239511e-03 2.086126731309433102e-02 2.101431592718736377e-02 1.015175119958449557e+00
2.915053062825178731e+00 6.136232207616530016e-03 1.836319032089894082e-02 1.847656682533876590e-02 1.012732085838440987e+00
3.199267137797384475e+00 5.163836258208944451e-03 1.613045410381486189e-02 1.621418148204904763e-02 1.010673928340949601e+00
3.511191734215134641e+00 4.340743638728428860e-03 1.414089221804774518e-02 1.420254181106790983e-02 1.008941858257746071e+00
3.853528593710531247e+00 3.645126362783052566e-03 1.237309139198952516e-02 1.241835787566458738e-02 1.007485678201927470e+00
4.229242874389498752e+00 3.058098265232835674e-03 1.080659458162784928e-02 1.083974358267359231e-02 1.006262603158349256e+00
4.641588833612781961e+00 2.563375287914749449e-03 9.422048248446359547e-03 9.446262564467270209e-03 1.005236229960698058e+00
5.094138014816380178e+00 2.146961818045753714e-03 8.201300831699996977e-03 8.218946596226652662e-03 1.004375638325215370e+00
5.590810182512228721e+00 1.796864401676819207e-03 7.127459512744253357e-03 7.140289645025064795e-03 1.003654608288723438e+00
6.135907273413176100e+00 1.502833036846241629e-03 6.184912179640193221e-03 6.194221059686219263e-03 1.003050940631458499e+00
6.734150657750828550e+00 1.256129416496231744e-03 5.359321085198225434e-03 5.366061552965635181e-03 1.002545868280674712e+00
7.390722033525782386e+00 1.049320897063753645e-03 4.637594128256958097e-03 4.642465564387569799e-03 1.002123547899007150e+00
8.111308307896871739e+00 8.760985726415443852e-04 4.007839040060852225e-03 4.011353381032336379e-03 1.001770621931099114e+00
8.902150854450392004e+00 7.311175972312538044e-04 3.459305075452025029e-03 3.461836084732303025e-03 1.001475842353969581e+00
9.770099572992256398e+00 6.098577835610856975e-04 2.982316130193813676e-03 2.984136028778169598e-03 1.001229748273481990e+00
1.072267222010324339e+01 5.085024857265418646e-04 2.568198558211786092e-03 2.569505157971702699e-03 1.001024390342212467e+00
1.176811952434999142e+01 4.238338189898146926e-04 2.209206363870943809e-03 2.210143097260368946e-03 1.000853095747656019e+00
1.291549665014884063e+01 3.531423626112195082e-04 1.898445905251861660e-03 1.899116563762063030e-03 1.000710268235453171e+00
1.417474162926806258e+01 2.941496140336472844e-04 1.629801769686371445e-03 1.630281316306329354e-03 1.000591218290415529e+00
1.555676143930472222e+01 2.449416022846627201e-04 1.397865073264147225e-03 1.398207552462374078e-03 1.000492019198784943e+00
1.707352647470692020e+01 2.039122154541811568e-04 1.197865089325737776e-03 1.198109398467535813e-03 1.000409385259195982e+00
1.873817422860384951e+01 1.697149445144921530e-04 1.025604822828627338e-03 1.025778912839904336e-03 1.000340568898805627e+00
2.056512308348653661e+01 1.412218886434813325e-04 8.774009124395309089e-04 8.775248381545070310e-04 1.000283273887397018e+00
2.257019719633921540e+01 1.174890023167391259e-04 7.500280542213627740e-04 7.501161846235159100e-04 1.000235582228965248e+00
2.477076355991711409e+01 9.772668974813717899e-05 6.406679936696190640e-04 6.407306101511695516e-04 1.000195892650883600e+00
2.718588242732942817e+01 8.127496673494111722e-05 5.468630206148762055e-04 5.469074705013485076e-04 1.000162868908974767e+00
2.983647240283340452e+01 6.758251325267230093e-05 4.664738185351530327e-04 4.665053461388910211e-04 1.000135396386632403e+00
3.274549162877731590e+01 5.618903238825395105e-05 3.976414610106828720e-04 3.976638053551217610e-04 1.000112545691407950e+00
3.593813663804629499e+01 4.671041287030923514e-05 3.387533088667614998e-04 3.387691329128275490e-04 1.000093542146996972e+00
3.944206059437659917e+01 3.882626423492733852e-05 2.884125380311216517e-04 2.884237364473246179e-04 1.000077740245859159e+00
4.328761281083061618e+01 3.226945635413461537e-05 2.454110168894071487e-04 2.454189364250647289e-04 1.000064602271141467e+00
4.750810162102798273e+01 2.681734950485998396e-05 2.087052500933283725e-04 2.087108471650623025e-04 1.000053680419153546e+00
5.214008287999689628e+01 2.228444822559633360e-05 1.773951109808273162e-04 1.773990642210892079e-04 1.000044601858117055e+00
5.722367659350220492e+01 1.851625272659538867e-05 1.507050950545896890e-04 1.507078855998873746e-04 1.000037056247741862e+00
6.280291441834259558e+01 1.538411637038455094e-05 1.279678407579346977e-04 1.279698094603753228e-04 1.000030785319544702e+00
6.892612104349701951e+01 1.278094744641159475e-05 1.086096798307477975e-04 1.086110679830997887e-04 1.000025574181624632e+00
7.564633275546290747e+01 1.061761878521274297e-05 9.213799684057929497e-05 9.213897513709253819e-05 1.000021244065547821e+00
8.302175681319752698e+01 8.819970285347690018e-06 7.813019531097530469e-05 7.813088442305428220e-05 1.000017646278491767e+00
9.111627561154895716e+01 7.326307686650936971e-06 6.622408564389649691e-05 6.622457082547878523e-05 1.000014657162157317e+00
1.000000000000000000e+02 6.085316380124952305e-06 5.610952734205624707e-05 5.610986878835989866e-05 1.000012173892210088e+00
47 changes: 47 additions & 0 deletions benchmarks/test_haloprofile.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,3 +128,50 @@ def test_haloprofile(model):
tol = np.clip(np.abs(HALOPROFILE_TOLERANCE * data[:, 1]), 1e-12, np.inf)
err = np.abs(prof - data[:, 1])
assert np.all(err <= tol)


def test_weak_lensing_functions():
data = np.loadtxt("./benchmarks/data/haloprofile_nfw_wl_numcosmo.txt")
z_lens = 1.0
z_source = 2.0
a_lens = 1.0 / (1.0 + z_lens)
a_source = 1.0 / (1.0 + z_source)
halomass = 1.0e15
concentration = 5
mDelta = "vir" # 200
rmin = 0.01
rmax = 100
r = np.exp(
np.log(rmin) +
np.log(rmax/rmin) * np.arange(data.shape[0]) / (data.shape[0]-1))

r_al = r / a_lens

mdef = ccl.halos.MassDef(mDelta, 'matter')
c = ccl.halos.ConcentrationConstant(c=concentration, mdef=mdef)
p = ccl.halos.HaloProfileNFW(
c, truncated=False, projected_analytic=True, cumul2d_analytic=True
)

kappa = p.convergence(COSMO, r_al, halomass,
a_lens, a_source, mass_def=mdef)
tol = np.clip(np.abs(HALOPROFILE_TOLERANCE * data[:, 1]), 1e-12, np.inf)
err_kappa = np.abs(kappa - data[:, 1])
assert np.all(err_kappa <= tol)

gamma = p.shear(COSMO, r_al, halomass, a_lens, a_source, mass_def=mdef)
tol = np.clip(np.abs(HALOPROFILE_TOLERANCE * data[:, 2]), 1e-12, np.inf)
err_gamma = np.abs(gamma - data[:, 2])
assert np.all(err_gamma <= tol)

gt = p.reduced_shear(COSMO, r_al,
halomass, a_lens, a_source, mass_def=mdef)
tol = np.clip(np.abs(HALOPROFILE_TOLERANCE * data[:, 3]), 1e-12, np.inf)
err_gt = np.abs(gt - data[:, 3])
assert np.all(err_gt <= tol)

mu = p.magnification(COSMO, r_al,
halomass, a_lens, a_source, mass_def=mdef)
tol = np.clip(np.abs(HALOPROFILE_TOLERANCE * data[:, 4]), 1e-12, np.inf)
err_mu = np.abs(mu - data[:, 4])
assert np.all(err_mu <= tol)
35 changes: 35 additions & 0 deletions pyccl/background.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -269,3 +269,38 @@ def rho_x(cosmo, a, species, is_comoving=False):
return _vectorize_fn4(
lib.rho_x, lib.rho_x_vec, cosmo, a,
species_types[species], int(is_comoving))


def sigma_critical(cosmo, a_lens, a_source):
"""Returns the critical surface mass density.
.. math::
\\Sigma_{\\mathrm{crit}} = \\frac{c^2}{4\\pi G}
\\frac{D_{\\rm{s}}}{D_{\\rm{l}}D_{\\rm{ls}}},
where :math:`c` is the speed of light, :math:`G` is the
gravitational constant, and :math:`D_i` is the angular diameter
distance. The labels :math:`i =` s, l and ls denotes the distances
to the source, lens, and between source and lens, respectively.
Args:
cosmo (:class:`~pyccl.core.Cosmology`): A Cosmology object.
a_lens (float): lens' scale factor.
a_source (float or array_like): source's scale factor.
Returns:
float or array_like: :math:`\\Sigma_{\\mathrm{crit}}` in units
of :math:`\\M_{\\odot}/Mpc^2`
"""
physical_constants = lib.cvar.constants
Ds = angular_diameter_distance(cosmo, a_source, a2=None)
Dl = angular_diameter_distance(cosmo, a_lens, a2=None)
Dls = angular_diameter_distance(cosmo, a_lens, a_source)
A = (
physical_constants.CLIGHT**2
* physical_constants.MPC_TO_METER
/ (4.0 * np.pi * physical_constants.GNEWT
* physical_constants.SOLAR_MASS)
)

Sigma_crit = A * Ds / (Dl * Dls)
return Sigma_crit

0 comments on commit c3d04fc

Please sign in to comment.