NFW and TNFW implementation

[andreasfilipp]

The NFW profile is currently defined by the concentration c, and the TNFW by the scale radius r_s. Ideally we would like to have both parameterizations available for both lens classes.

[andreasfilipp]

Further, TNFW misbehaves for big truncation radii.

TNFW should give similar kappa map to NFW, especially for high truncation radii, but e.g.

tnfw = caustics.TNFW(cosmology=cosmology, tau=500.)

nfw = caustics.NFW(cosmology=cosmology)

z_lens = torch.tensor([0.5])
z_s = 1.0

m200 = torch.tensor([1e11])

conc = concentration(m200)

r_scale = scale_radius_fct(m200, 0.5)
r_scale = torch.tensor(r_scale)

x0 = torch.tensor([-0.5]) 
y0 = torch.tensor([0.5]) 

params_tnfw = torch.stack((z_lens, x0, y0, m200, r_scale), dim=1)
params_nfw = torch.stack((z_lens, x0, y0, m200, conc), dim=1)

alpha_x_tnfw, alpha_y_tnfw = tnfw.reduced_deflection_angle(theta_x, theta_y, z_s, params_tnfw)
alpha_x_nfw, alpha_y_nfw = nfw.reduced_deflection_angle(theta_x, theta_y, z_s, params_nfw)

kappa_tnfw = tnfw.convergence(theta_x, theta_y, z_s, params_tnfw)
kappa_nfw = nfw.convergence(theta_x, theta_y, z_s, params_nfw) 

with

 def concentration( m_200):  # m in Msun!!!
     return 6 * (m_200 / (1e12 ))**(-0.098) 

 def scale_radius_fct( m, z): # m in Msun!!!
     c = concentration(m)
     # Compute r_200 from m_200
     rho_c = cosmology.critical_density(z)
     r_200 = (3 * m / (4 * torch.pi * 200 * rho_c)) ** (1 / 3)
     # Use NFW Definition
     rs = r_200 / c
     # Convert to angular size
     D_z = cosmology.angular_diameter_distance(z) 
     rs = (rs / D_z) * 180 / torch.pi * 60**2 #*M_s**(1/3)
     return rs # arcsec

cosmology = caustics.cosmology.FlatLambdaCDM()

returns very different kappa maps. (same features, but very different scales)