F14-Radial distribution.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "import sys\n",
    "sys.path.insert(0, \"/home/yymao/Dropbox/Academia/Collaborations/SAGA-local/saga\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from itertools import combinations\n",
    "from collections import Counter, defaultdict\n",
    "from astropy.table import Table, join, vstack\n",
    "from astropy.coordinates import SkyCoord, Distance\n",
    "from SAGA.utils import add_skycoord\n",
    "from easyquery import Query, QueryMaker\n",
    "from astropy.io import ascii\n",
    "\n",
    "import SAGA\n",
    "from SAGA.database import FitsTable, GoogleSheets\n",
    "from SAGA import ObjectCuts as C\n",
    "from SAGA.utils import add_skycoord\n",
    "from SAGA.utils.distance import z2v, d2m, m2d\n",
    "from SAGA.utils.display import show_images\n",
    "from SAGA.objects.object_catalog import calc_fiducial_p_sat, calc_fiducial_p_sat_corrected\n",
    "print(SAGA.__version__)\n",
    "\n",
    "from scipy.stats import norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "saga = SAGA.QuickStart()\n",
    "saga.set_default_base_version(\"paper2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base = SAGA.database.FitsTable(\"/home/yymao/Documents/Research/SAGA/PaperII/data-archive/saga_base_all.fits\").read()\n",
    "sats = Query(C.is_sat, C.r_abs_limit, saga.host_catalog.construct_host_query(\"paper2_complete\")).filter(base)\n",
    "del base"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert len(sats) == 123\n",
    "assert len(np.unique(sats[\"HOSTID\"])) == 33"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbins = np.linspace(10, 300, 1001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_cdf(x, bins):\n",
    "    return np.searchsorted(np.sort(x), bins) / len(x)\n",
    "\n",
    "def print_range(arr):\n",
    "    a = np.round(arr).astype(np.int)\n",
    "    print(\"{}^{{+{}}}_{{-{}}}\".format(a[1], a[2]-a[1], a[1]-a[0]))\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "saga_median_r = np.array([np.median(sats_this['RHOST_KPC']) for sats_this in sats.group_by(\"HOSTID\").groups])\n",
    "print_range(np.quantile(saga_median_r, norm.sf([1, 0, -1])))\n",
    "\n",
    "mask1 = saga_median_r < np.quantile(saga_median_r, 0.25)\n",
    "mask2 = saga_median_r > np.quantile(saga_median_r, 0.75)\n",
    "\n",
    "saga_radial = np.vstack([get_cdf(sats_this['RHOST_KPC'], rbins) for sats_this in sats.group_by(\"HOSTID\").groups])\n",
    "saga_radial_mean = saga_radial.mean(axis=0)\n",
    "saga_radial_std = saga_radial.std(axis=0)\n",
    "saga_radial_pct = np.quantile(saga_radial, norm.sf([1, 0, -1]), axis=0)\n",
    "\n",
    "saga_radial_mean1 = saga_radial[mask1].mean(axis=0)\n",
    "saga_radial_mean2 = saga_radial[mask2].mean(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dhalf = [np.median(Query(\"RHOST_KPC < 150\").filter(sats_this, \"RHOST_KPC\")) for sats_this in sats.group_by(\"HOSTID\").groups if Query(\"RHOST_KPC < 150\").count(sats_this)]\n",
    "print_range(np.quantile(dhalf, norm.sf([1, 0, -1])))\n",
    "\n",
    "plt.hist(dhalf, np.linspace(10, 150, 6), alpha=0.6)\n",
    "plt.xlim(0, 160)\n",
    "plt.ylim(0, 10)\n",
    "plt.scatter(dhalf, np.ones_like(dhalf)*9.5, s=5)\n",
    "plt.xlabel(\"SAGA $d_\\\\mathrm{proj, half}$ (limited to sats $< 150$ kpc)\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import quad, cumtrapz\n",
    "from math import sqrt\n",
    "\n",
    "def g(c):\n",
    "    return 1.0 / (np.log1p(c) - c / (1 + c))\n",
    "\n",
    "def Cinv(x):\n",
    "    with np.errstate(invalid=\"ignore\"):\n",
    "        return np.where(x > 1, np.arccosh(x),  np.arccos(x))\n",
    "\n",
    "def nfw_proj(x, c):\n",
    "    y = c*x\n",
    "    return g(c) * (Cinv(1/y) / np.sqrt(np.abs(y*y-1)) + np.log(y/2))\n",
    "\n",
    "def nfw_rho(x, c=3):\n",
    "    y = c*x\n",
    "    return 1/y/(1+y)**2\n",
    "\n",
    "def Sigma(R, rho_func, r_lim=1):\n",
    "    return quad(lambda r: r * rho_func(r) / sqrt(r + R), R, r_lim, weight=\"alg\", wvar=(-0.5, 0))[0]\n",
    "\n",
    "def uniform(x, lim_inner=10, lim_outer=300):\n",
    "    ao = lim_outer*lim_outer\n",
    "    ai = lim_inner*lim_inner\n",
    "    return (x*x - ai) / (ao-ai)\n",
    "\n",
    "def get_theory(rho_func):\n",
    "    Rall = np.linspace(10, 300, 10001)[:-1] / 300\n",
    "    S = np.fromiter((Sigma(R, rho_func, r_lim=1) for R in Rall), dtype=np.float, count=len(Rall))\n",
    "    S_int = cumtrapz(S*Rall, Rall, initial=0)\n",
    "    S_int /= S_int[-1]\n",
    "    return Rall*300, S_int"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_th, s_th = get_theory(lambda r: 1/r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lg = SAGA.database.FitsTable(\"../data/McConnachie2012.fit\").read()\n",
    "lg[\"Mr\"] = lg[\"vmag_lc\"] - 0.4\n",
    "lg = (Query((np.isfinite, \"vmag_lc\"), ~QueryMaker.equals(\"Name\", \"Canis Major\"), \"Mr < -12.295\") | QueryMaker.equals(\"Name\", \"Andromeda\")).filter(lg)\n",
    "lg[\"Mr\"][lg[\"Name\"] == \"Andromeda\"] = -22"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xyz = SkyCoord(lg[\"RAJ2000\"], lg[\"DEJ2000\"], lg[\"D\"], unit=(\"deg\", \"deg\", \"kpc\")).galactocentric.cartesian.xyz.value.T\n",
    "xyz = np.vstack([xyz, [0, 0, 0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx_m31 = np.argmax(lg[\"Name\"] == \"Andromeda\")\n",
    "idx_mw = len(lg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_projections(xyz, idx_host, idx_other, z_limit=500, n_needed=10):\n",
    "    realization = []\n",
    "    xyz_shifted = xyz - xyz[idx_host]\n",
    "    while len(realization) < n_needed:\n",
    "        m_rand = np.linalg.qr(np.random.randn(3, 3))[0]\n",
    "        xyz_rot = np.dot(xyz_shifted, m_rand)\n",
    "        r_proj = np.hypot(xyz_rot[:,0], xyz_rot[:,1])\n",
    "        if r_proj[idx_other] < 300:\n",
    "            continue\n",
    "        realization.append(r_proj[(r_proj > 10) & (r_proj < 300) & (np.abs(xyz_rot[:,2]) < z_limit)])\n",
    "    return realization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mw_realization = generate_projections(xyz, idx_mw, idx_m31, n_needed=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m31_realization = generate_projections(xyz, idx_m31, idx_mw, n_needed=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mw_radial = np.vstack([get_cdf(r, rbins) for r in mw_realization])\n",
    "m31_radial = np.vstack([get_cdf(r, rbins) for r in m31_realization])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mw_radial_mean = mw_radial.mean(axis=0)\n",
    "mw_radial_std = mw_radial.std(axis=0)\n",
    "mw_radial_pct = np.quantile(mw_radial, norm.sf([1, 0, -1]), axis=0)\n",
    "\n",
    "m31_radial_mean = m31_radial.mean(axis=0)\n",
    "m31_radial_std = m31_radial.std(axis=0)\n",
    "m31_radial_pct = np.quantile(m31_radial, norm.sf([1, 0, -1]), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(4.8, 3.8))\n",
    "\n",
    "plt.plot(rbins, mw_radial_pct[1], c=\"C1\", label='MW', lw=1, ls=\"--\", zorder=20);\n",
    "plt.plot(rbins, m31_radial_pct[1], c=\"C7\", label='M31', lw=1, ls=\"-.\", zorder=10);\n",
    "\n",
    "plt.fill_between(rbins, mw_radial_pct[0], mw_radial_pct[2], color=\"C1\", alpha=0.25, lw=0, zorder=5, rasterized=True);\n",
    "plt.fill_between(rbins, m31_radial_pct[0], m31_radial_pct[2], color=\"C7\", alpha=0.25, lw=0, zorder=1, rasterized=True);\n",
    "\n",
    "\n",
    "plt.plot(rbins, saga_radial_mean, c=\"C2\", label='SAGA', lw=2, zorder=50)\n",
    "plt.plot(rbins, saga_radial_mean1, c=\"C2\", label='SAGA 25%\\nmost conc.', lw=0.8, zorder=40);\n",
    "\n",
    "plt.plot(r_th, s_th, c=\"C0\", label='$\\\\rho(r) = r^{-1}$', ls=\":\", zorder=30);\n",
    "\n",
    "\n",
    "(lines, labels) = plt.gca().get_legend_handles_labels()\n",
    "l = plt.Line2D([0], [0], linestyle='none', marker='')\n",
    "for _  in range(2):\n",
    "    lines.insert(0, l)\n",
    "    labels.insert(0, \"\")\n",
    "\n",
    "plt.axhline(0.0, c='k', ls=':', lw=0.5, zorder=0);\n",
    "plt.axhline(1.0, c='k', ls=':', lw=0.5, zorder=0);\n",
    "plt.axvline(10.0, c='k', ls=':', lw=0.5, zorder=0);\n",
    "plt.axvline(300.0, c='k', ls=':', lw=0.5, zorder=0);\n",
    "plt.legend(lines, labels, loc=\"lower right\", labelspacing=0.3, fontsize=13, handlelength=1, handletextpad=0.5,ncol=2, columnspacing=1)\n",
    "plt.xlabel('$d_\\\\mathrm{proj}$ [kpc]');\n",
    "plt.ylabel('$N(<d_\\\\mathrm{proj}) \\, / \\, N(< 300\\,\\\\mathrm{kpc})$');\n",
    "plt.tight_layout()\n",
    "plt.savefig('/home/yymao/Downloads/sat_radial.pdf', bbox_inches='tight', dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}