F15-Plane of satellites.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",
    "\n",
    "from astropy.table import Table\n",
    "from astropy.coordinates import SkyCoord, Distance\n",
    "from SAGA.utils import add_skycoord\n",
    "from easyquery import Query, QueryMaker\n",
    "\n",
    "import SAGA\n",
    "from SAGA.database import FitsTable\n",
    "from SAGA import ObjectCuts as C\n",
    "from SAGA.utils import add_skycoord\n",
    "from SAGA.utils.distance import z2v\n",
    "\n",
    "print(SAGA.__version__)"
   ]
  },
  {
   "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",
    "base = saga.host_catalog.construct_host_query(\"paper2_complete\").filter(base)\n",
    "assert len(np.unique(base[\"HOSTID\"])) == 36\n",
    "\n",
    "sats = Query(C.is_sat).filter(base)\n",
    "del base"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert len(sats) == 127\n",
    "assert C.r_abs_limit.count(sats) == 123\n",
    "assert len(np.unique(sats[\"HOSTID\"])) == 34"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sats = add_skycoord(sats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdss_x = np.linspace(0, 180, 19)[:-1] + 5\n",
    "sdss_y = np.array([0.64268367, 0.44155338, 0.4841401 , 0.4633828 , 0.4685856 ,\n",
    "       0.58841093, 0.54086203, 0.64018368, 0.48222288, 0.51515562,\n",
    "       0.52488057, 0.55817479, 0.51704362, 0.47159272, 0.44232658,\n",
    "       0.47630996, 0.58469492, 0.49163754])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "X = []\n",
    "\n",
    "for sats_this in sats.group_by('HOSTID').groups:\n",
    "    if len(sats_this) < 2:\n",
    "        continue\n",
    "\n",
    "    host_this = str(sats_this['HOSTID'][0])\n",
    "    host_coord = SkyCoord(sats_this['HOST_RA'][0], sats_this['HOST_DEC'][0], unit=\"deg\")\n",
    "    \n",
    "    sep_deg = host_coord.position_angle(sats_this['coord']).degree\n",
    "    dv = z2v(sats_this['SPEC_Z']) - sats_this['HOST_VHOST']\n",
    "    \n",
    "    for i, j in combinations(range(len(sats_this)), 2):\n",
    "        sat1 = sats_this[i]\n",
    "        sat2 = sats_this[j]\n",
    "        X.append((host_this, sep_deg[i], sep_deg[j], dv[i], dv[j], sat1[\"RHOST_KPC\"], sat2[\"RHOST_KPC\"], \n",
    "                  1000*sat1[\"HOST_DIST\"]*sat1['coord'].separation(sat2['coord']).radian))\n",
    "    \n",
    "X = Table(np.array(X, dtype=np.dtype([('HOSTID', '<U12'), ('pa1', np.float), ('pa2', np.float), ('dv1', np.float), ('dv2', np.float), ('d1', np.float), ('d2', np.float), ('sep', np.float)])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X[\"sign\"] = np.sign(X[\"dv1\"]) * np.sign(X[\"dv2\"])\n",
    "dpa = np.abs(X[\"pa1\"] - X[\"pa2\"])\n",
    "X[\"opening_angle\"] = np.where(dpa < 180, 180-dpa, dpa-180)\n",
    "X[\"corotating\"] = np.where(X[\"opening_angle\"] < 90, X[\"sign\"] == -1, X[\"sign\"] == 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dv_cut = Query('abs(dv1) > 25', 'abs(dv2) > 25')\n",
    "#dv_cut = Query()\n",
    "radius_cut = Query(\"d1<150\", \"d2<150\")\n",
    "\n",
    "print(dv_cut.count(X))\n",
    "print((dv_cut & radius_cut).count(X))\n",
    "print()\n",
    "print((~dv_cut).count(X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(4.8,4))\n",
    "\n",
    "bins = np.linspace(0, 180, 11)\n",
    "hist_co = np.histogram(Query('corotating', dv_cut).filter(X, 'opening_angle'), bins)[0]\n",
    "hist_all = np.histogram(dv_cut.filter(X, 'opening_angle'), bins)[0]\n",
    "\n",
    "p = (hist_co+1)/(hist_all+2)\n",
    "perr = np.sqrt(p*(1-p) / hist_all)\n",
    "p1 = np.minimum(1, p + perr)\n",
    "p2 = np.maximum(0, p - perr)\n",
    "plt.fill_between((bins[1:]+bins[:-1])*0.5, p1, p2, lw=0, color=\"C0\", alpha=0.2, label=\"SAGA sats$<300$ kpc\");\n",
    "\n",
    "bins = np.linspace(0, 180, 11)\n",
    "hist_co_inner = np.histogram(Query('corotating', dv_cut, radius_cut).filter(X, 'opening_angle'), bins)[0]\n",
    "hist_all_inner = np.histogram(Query(dv_cut, radius_cut).filter(X, 'opening_angle'), bins)[0]\n",
    "\n",
    "p = (hist_co_inner+1)/(hist_all_inner+2)\n",
    "perr = np.sqrt(p*(1-p) / hist_all_inner)\n",
    "p1 = np.minimum(1, p + perr)\n",
    "p2 = np.maximum(0, p - perr)\n",
    "\n",
    "plt.errorbar((bins[1:]+bins[:-1])*0.5, p, yerr=[p-p2, p1-p], lw=1, ls=\"\", color='C2', label=\"SAGA sats$<150$ kpc\", marker='s');\n",
    "plt.errorbar(sdss_x, sdss_y, c=\"C7\", alpha=0.8, label=\"SDSS (Phillips+ 2015, Fig. 3b)\")\n",
    "\n",
    "plt.text(3, 1.04, r\"Opposite side\")\n",
    "plt.text(177, 1.04, r\"Same side\", ha=\"right\")\n",
    "\n",
    "plt.legend(fontsize='small', loc=\"lower left\")\n",
    "plt.xlim(0, 180.0);\n",
    "plt.xticks([0,30,60,90,120,150,180])\n",
    "plt.axhline(0.5, c='k', ls='--', lw=0.5);\n",
    "plt.xlabel('Opening angle [deg]');\n",
    "plt.ylabel('Co-rotating fraction')\n",
    "plt.ylim(-0.15, 1.15)\n",
    "plt.yticks([0, 0.2, 0.4, 0.6, 0.8, 1.0])\n",
    "plt.tight_layout()\n",
    "plt.savefig('/home/yymao/Downloads/sat_coroating_pairs.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## additional plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(ncols=3, gridspec_kw={'wspace':0.1, 'width_ratios':(30,1,1)})\n",
    "t = Query('corotating', dv_cut, 'abs(dv1-dv2) > 30').filter(X) \n",
    "CS = ax[0].scatter(t[\"opening_angle\"], t[\"sep\"], c=np.abs(t[\"dv1\"]-t[\"dv2\"]), vmin=0, vmax=400, cmap=\"Blues_r\", label=\"Co-rotating pair\")\n",
    "plt.colorbar(CS, ax[1])\n",
    "\n",
    "t = Query('corotating == 0', dv_cut, 'abs(dv1-dv2) > 30').filter(X)\n",
    "CS = ax[0].scatter(t[\"opening_angle\"], t[\"sep\"], c=np.abs(t[\"dv1\"]-t[\"dv2\"]), vmin=0, vmax=400, cmap=\"Oranges_r\", label=\"Counter-rotating pair\")\n",
    "cbar = plt.colorbar(CS, ax[2])\n",
    "ax[1].set_yticklabels([])\n",
    "ax[2].set_xlabel(r\"$|\\Delta v|$ [km s$^{-1}$]\", fontsize='small')\n",
    "\n",
    "lgnd = ax[0].legend(markerfirst=False, scatterpoints=3, frameon=True)\n",
    "lgnd.legendHandles[0].set_color(plt.cm.Blues_r([0.2, 0.5, 0.7]))\n",
    "lgnd.legendHandles[1].set_color(plt.cm.Oranges_r([0.2, 0.5, 0.7]))\n",
    "ax[0].axvline(90, color=\"grey\", ls='--')\n",
    "ax[0].set_xlim(-5, 185.0);\n",
    "ax[0].set_xlabel('Opening angle [deg]');\n",
    "ax[0].set_ylabel('Projected separation [kpc]')\n",
    "ax[0].set_ylim(0, 600.0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}