Skip to content

Commit

Permalink
incorporating original demo_mock_cluster.ipynb from main to remove ch…
Browse files Browse the repository at this point in the history 
…anges to the file in this branch that should not have been committed
  • Loading branch information
@cavestruz
cavestruz committed Jul 19, 2024
1 parent 7fe44dc commit 18325ff
Showing 1 changed file with 2 additions and 308 deletions.
310 changes: 2 additions & 308 deletions examples/demo_mock_cluster.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -858,318 +858,12 @@
" ngals=ngals\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ideal_data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from astropy import wcs\n",
"dx = -1.0 / 60.0\n",
"dy = 1.0 / 60.0\n",
"\n",
"w = wcs.WCS(naxis = 2)\n",
"\n",
"w.wcs.crval = [cluster_ra, cluster_dec]\n",
"w.wcs.cd = np.array([[dx, 0.0], [0.0, dy]])\n",
"w.wcs.ctype = ['RA---TAN', 'DEC--TAN']\n",
"w.wcs.crpix = [0, 0]\n",
"\n",
"coord = np.stack([ideal_data['ra'], ideal_data['dec']], axis = 1)\n",
" \n",
"xy = w.wcs_world2pix(coord, 1)\n",
"\n",
"x = xy[:,0]\n",
"y = xy[:,1]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"ltmin = 0.0\n",
"ltmax = 1.28\n",
"dlt = 0.16\n",
"\n",
"nbin = int((ltmax - ltmin + 1.0e-30) / dlt)\n",
"\n",
"lt1 = np.linspace(ltmin, ltmax - dlt, num = nbin)\n",
"lt2 = np.linspace(ltmin + dlt, ltmax, num = nbin)\n",
"\n",
"tbin = 10 ** (0.5 * (lt1 + lt2))\n",
"\n",
"t = np.sqrt(x * x + y * y)\n",
"lt = np.log10(t)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cos2 = (x * x - y * y) / (t * t + 1.0e-30)\n",
"sin2 = 2.0 * x * y / (t * t + 1.0e-30)\n",
"\n",
"gt = (-1.0) * cos2 * ideal_data['e1'] - sin2 * ideal_data['e2']\n",
"gx = sin2 * ideal_data['e1'] - cos2 * ideal_data['e2']"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"gt_ave = np.zeros(nbin)\n",
"gx_ave = np.zeros(nbin)\n",
"gt_sig = np.zeros(nbin)\n",
"gx_sig = np.zeros(nbin)\n",
"weight = np.ones(np.shape(gt)[0])\n",
"for j in range(nbin):\n",
" mask = np.logical_and((lt > lt1[j]), (lt < lt2[j]))\n",
" ct = t[mask]\n",
" cgt = gt[mask]\n",
" cgx = gx[mask]\n",
" cwei = weight[mask]\n",
"\n",
" wsum1 = np.sum(cwei)\n",
" wsum2 = np.sum(cwei * cwei)\n",
"\n",
" gt_ave[j] = np.sum(cgt * cwei) / wsum1\n",
" gx_ave[j] = np.sum(cgx * cwei) / wsum1\n",
"\n",
" gt_sig[j] = np.sqrt((np.sum(cgt * cgt * cwei) / wsum1) * (wsum2 / (wsum1 * wsum1)))\n",
" gx_sig[j] = np.sqrt((np.sum(cgx * cgx * cwei) / wsum1) * (wsum2 / (wsum1 * wsum1)))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib as mpl"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"xmin = 1.0\n",
"xmax = 20.0\n",
"ymin = -0.08\n",
"ymax = 0.15\n",
"\n",
"plt.figure(figsize = (8, 6))\n",
"\n",
"\n",
"plt.plot(tbin * 1.005, gt_ave, '.', markersize = 10, color = 'r')\n",
"plt.errorbar(tbin * 1.005, gt_ave, yerr = gt_sig, fmt = 'none', ecolor = 'r', label = '_nolegend_')\n",
"\n",
"plt.plot(tbin * 0.995, gx_ave, 'x', markersize = 10, color = 'b')\n",
"plt.errorbar(tbin * 0.995, gx_ave, yerr = gt_sig, fmt = 'none', ecolor = 'b', label = '_nolegend_')\n",
"\n",
"plt.legend(['$\\gamma_+$', '$\\gamma_\\\\times$'], markerscale = 1, frameon = False, labelspacing = 0.3, numpoints = 1, loc = 'upper right')\n",
"\n",
"plt.hlines(0, xmin, xmax, linestyles = 'dotted', color = 'black')\n",
"ax = plt.gca()\n",
"ax.set_xscale('log')\n",
"\n",
"ax.yaxis.set_major_locator(mpl.ticker.MultipleLocator(0.1))\n",
"ax.yaxis.set_minor_locator(mpl.ticker.MultipleLocator(0.02))\n",
"\n",
"plt.xlim([xmin,xmax])\n",
"plt.ylim([ymin,ymax])\n",
"\n",
"plt.xlabel('$\\\\theta$ [arcmin]')\n",
"plt.ylabel('$\\gamma_+$ or $\\gamma_\\\\times$')\n",
"#plt.title('ideal data')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from scipy.optimize import curve_fit"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def fit_sis(x, t_ein):\n",
" return 0.5 * (t_ein / x)\n",
"\n",
"# fit parameter t_ein (Einstein radius) with observed shear profile\n",
"tein, cov = curve_fit(fit_sis, tbin, gt_ave, sigma = gt_sig)\n",
"\n",
"print('Einstein radius: %f arcmin' % tein[0])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# conpute best-fit profile\n",
"tt1 = 10 ** np.linspace(-1.0, 2.0, num = 300)\n",
"yy1 = fit_sis(tt1, tein[0])\n",
"\n",
"xmin = 1.0\n",
"xmax = 20.0\n",
"ymin = 0.005\n",
"ymax = 0.3\n",
"\n",
"plt.figure(figsize = (8, 6))\n",
"\n",
"plt.hlines(0, xmin, xmax, linestyles = 'dotted')\n",
"\n",
"plt.plot(tt1, yy1, '-', color = 'b', lw = 2.)\n",
"\n",
"plt.plot(tbin, gt_ave, '.', markersize = 10, color = 'r')\n",
"plt.errorbar(tbin, gt_ave, yerr = gt_sig, fmt = 'none', ecolor = 'r', label = '_nolegend_')\n",
"\n",
"plt.legend(['SIS fit'], markerscale = 1, frameon = False, labelspacing = 0.3, numpoints = 1, loc = 'upper right')\n",
"\n",
"ax = plt.gca()\n",
"ax.set_xscale('log')\n",
"ax.set_yscale('log')\n",
"\n",
"plt.xlim([xmin,xmax])\n",
"plt.ylim([ymin,ymax])\n",
"\n",
"plt.xlabel('$\\\\theta$ [arcmin]')\n",
"plt.ylabel('$\\gamma_+$')\n",
"plt.title('ideal data')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from astropy import constants as const\n",
"from astropy import units as u\n",
"from astropy.cosmology import FlatLambdaCDM\n",
"cosmo = FlatLambdaCDM(H0 = 70, Om0 = 0.27 - 0.045)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"tein_rad = tein[0] * u.arcmin\n",
"dos = cosmo.angular_diameter_distance(src_z)\n",
"dol = cosmo.angular_diameter_distance(cluster_z)\n",
"dls = cosmo.angular_diameter_distance_z1z2(cluster_z, src_z)\n",
"sig = const.c * np.sqrt((tein_rad * dos).to(u.Mpc, u.dimensionless_angles()) / (4.0 * np.pi * dls))\n",
"\n",
"print('velocity dispersion: ', sig.to(u.km/u.s))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from colossus.halo import concentration\n",
"from colossus.halo import profile_nfw\n",
"from colossus.cosmology import cosmology\n",
"params = {'flat': True, 'H0': 70.0, 'Om0': 0.27 - 0.045, 'Ob0': 0.045, 'sigma8': 0.81, 'ns': 0.95}\n",
"cosmology.addCosmology('myCosmo', **params)\n",
"cosmo_c = cosmology.setCosmology('myCosmo')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def calc_R200m(halo_mass,redshift=0.25):\n",
" concen = concentration.concentration(halo_mass,'200m',redshift,model='diemer15')\n",
" \n",
" profile = profile_nfw.NFWProfile(M=halo_mass,mdef='200m',z=redshift,c=concen)\n",
" R200m = profile.RDelta(redshift,'200m')\n",
" R200m /= 1000. #ts: convert to Mpc/h \n",
" return concen,R200m"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def convert_M200m_to_Mvir(halo_mass,redshift=0.25):\n",
" concen = concentration.concentration(halo_mass,'200m',redshift,model='diemer15')\n",
" profile = profile_nfw.NFWProfile(M=halo_mass,mdef='200m',z=redshift,c=concen)\n",
" Rvir = profile.RDelta(redshift,'vir')\n",
" Mvir = profile.enclosedMass(Rvir)\n",
" return concen,Rvir/1000.,Mvir"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"concen,R200m = calc_R200m(10**15.,redshift=cluster_z)\n",
"cvir,Rvir,Mvir = convert_M200m_to_Mvir(10**15.,redshift=cluster_z)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"G = 4.299*10**(-9) #unit:[Mpc(km/s)^2 Msun^-1] \n",
"Mhalo = 10**15.\n",
"hubble = cosmo.H(cluster_z)\n",
"a_z = 1./(1.+cluster_z)\n",
"sigma = G*(Mhalo)/R200m/a_z**2. #TS:[(Mpc/h)^2]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"np.sqrt(sigma)"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
Expand All @@ -1183,7 +877,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.9"
"version": "3.12.4"
}
},
"nbformat": 4,
Expand Down

0 comments on commit 18325ff

Please sign in to comment.