SAGA-halo-ellipticities.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import os\n",
    "from astropy import units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "import SAGA\n",
    "from SAGA import ObjectCuts as C\n",
    "from dotenv import load_dotenv\n",
    "\n",
    "load_dotenv()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "saga_dir = os.getenv('SAGA_DIR')\n",
    "saga_db = os.getenv('SAGA_DB')\n",
    "\n",
    "print(saga_db)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load SAGA database\n",
    "saga = SAGA.QuickStart(local_dir=saga_dir)\n",
    " \n",
    "saga.database[\"combined_base\"].remote.path = saga_db\n",
    "saga.database[\"combined_base\"].download(overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load catalogs\n",
    "base = saga.object_catalog.load_combined_base_catalog()\n",
    "base = saga.host_catalog.construct_host_query(\"paper3\").filter(base)\n",
    "\n",
    "sats = C.is_sat.filter(base)\n",
    "hosts = C.is_host.filter(base)\n",
    "\n",
    "df = pd.DataFrame(dict(sats))\n",
    "\n",
    "saga_sats = df\n",
    "# print(df.columns.tolist())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sns.set(style=\"darkgrid\", palette=\"viridis\")\n",
    "custom_style = {\n",
    "    \"figure.facecolor\": \"212946\",\n",
    "    \"axes.facecolor\": \"212946\",\n",
    "    \"savefig.facecolor\": \"212946\", \n",
    "    \"grid.color\": \"2A3459\",\n",
    "    \"text.color\": \"0.9\",\n",
    "    \"axes.labelcolor\": \"0.9\",\n",
    "    \"xtick.color\": \"0.9\",\n",
    "    \"ytick.color\": \"0.9\",\n",
    "    \"grid.linestyle\": \"-\",\n",
    "    \"lines.solid_capstyle\": \"round\"\n",
    "}\n",
    "\n",
    "sns.set_style(\"darkgrid\", custom_style)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Ellipticity calculations</h2>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>Get cartesian data points</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['X'] = SkyCoord(ra=df.RA*u.degree, dec=df.DEC*u.degree, frame='icrs').transform_to(SkyCoord(ra=df.HOST_RA*u.degree, dec=df.HOST_DEC*u.degree, frame='icrs').skyoffset_frame()).cartesian.y\n",
    "df['Y'] = SkyCoord(ra=df.RA*u.degree, dec=df.DEC*u.degree, frame='icrs').transform_to(SkyCoord(ra=df.HOST_RA*u.degree, dec=df.HOST_DEC*u.degree, frame='icrs').skyoffset_frame()).cartesian.z\n",
    "df['HOST_X'] = SkyCoord(ra=df.HOST_RA*u.degree, dec=df.HOST_DEC*u.degree, frame='icrs').transform_to(SkyCoord(ra=df.HOST_RA*u.degree, dec=df.HOST_DEC*u.degree, frame='icrs').skyoffset_frame()).cartesian.y\n",
    "df['HOST_Y'] = SkyCoord(ra=df.HOST_RA*u.degree, dec=df.HOST_DEC*u.degree, frame='icrs').transform_to(SkyCoord(ra=df.HOST_RA*u.degree, dec=df.HOST_DEC*u.degree, frame='icrs').skyoffset_frame()).cartesian.z\n",
    "df['X_REFL'] = -1*df['X']\n",
    "df['HOST_X_REFL'] = -1*df['HOST_X']\n",
    "\n",
    "\n",
    "df[['OBJID', 'RA', 'DEC', 'X', 'Y', 'HOSTID', 'HOST_RA', 'HOST_DEC', 'HOST_X', 'HOST_Y']]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h3>Ellipticity calculations</h3>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>Filter data</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = df\n",
    "data = df.groupby(['HOSTID']).filter(lambda x: len(x) > 2).copy()\n",
    "\n",
    "saga_sats = df\n",
    "%store saga_sats\n",
    "\n",
    "len(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parent_ids = data['HOSTID'].unique()\n",
    "avg_num_satellites = len(data)/len(parent_ids)\n",
    "print(\"Average # satellites per SAGA system: {}\".format(avg_num_satellites))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>Quadrupole moments</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ellipticity_df = data.copy()\n",
    "ellipticity_df['x_diff_sq'] = (ellipticity_df['X'] - ellipticity_df['HOST_X'])**2\n",
    "ellipticity_df['y_diff_sq'] = (ellipticity_df['Y'] - ellipticity_df['HOST_Y'])**2\n",
    "ellipticity_df['x_diff_y_diff'] = (ellipticity_df['Y'] - ellipticity_df['HOST_Y'])*(ellipticity_df['X'] - ellipticity_df['HOST_X'])\n",
    "ellipticity_df = ellipticity_df.groupby('HOSTID', as_index = False)['x_diff_sq', 'y_diff_sq', 'x_diff_y_diff'].mean().rename(columns={'x_diff_sq':'Q_xx', 'y_diff_sq':'Q_yy', 'x_diff_y_diff':'Q_xy'})\n",
    "ellipticity_df"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h4>Ellipticity components</h4>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def e_1 (Q_xx, Q_yy, Q_xy):\n",
    "    return (( Q_xx - Q_yy )/(Q_xx + Q_yy + 2*(Q_xx*Q_yy-Q_xy**2)**(1/2)))\n",
    "def e_2 (Q_xx, Q_yy, Q_xy):\n",
    "    return (2*Q_xy/(Q_xx + Q_yy + 2*(Q_xx*Q_yy-Q_xy**2)**(1/2)))\n",
    "def e_3 (e_1, e_2):\n",
    "    return ( e_1**2 + e_2**2 )**(1/2)\n",
    "\n",
    "\n",
    "ellipticity_df['e_1'] = e_1(ellipticity_df['Q_xx'], ellipticity_df['Q_yy'], ellipticity_df['Q_xy'])\n",
    "ellipticity_df['e_2'] = e_2(ellipticity_df['Q_xx'], ellipticity_df['Q_yy'], ellipticity_df['Q_xy'])\n",
    "ellipticity_df['e_3'] = e_3(ellipticity_df['e_1'], ellipticity_df['e_2'])\n",
    "\n",
    "ellipticity_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = sns.histplot(ellipticity_df[\"e_3\"], bins=7, element='poly', stat=\"density\", fill=False)\n",
    "ax.set_xlabel('$\\longleftarrow$ more elliptical         Ellipticity (2D)         more spherical $\\longrightarrow$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ellipticity_df_SAGA = ellipticity_df\n",
    "%store ellipticity_df_SAGA"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h2>Highest luminosity analysis</h2>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "highest_luminosity_sats = df.groupby('HOSTID', as_index = False)['Mr'].min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "highest_luminosity_ellipticity = pd.merge(highest_luminosity_sats, ellipticity_df, how=\"left\", left_on=\"HOSTID\", right_on=\"HOSTID\")[['HOSTID', 'Mr', 'e_3']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "highest_luminosity_ellipticity_SAGA = highest_luminosity_ellipticity[~highest_luminosity_ellipticity['e_3'].isna()]\n",
    "%store highest_luminosity_ellipticity_SAGA\n",
    "highest_luminosity_ellipticity_SAGA"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "satellite-ellipticities",
   "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.9.16"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "593019684f7021e8fe6a4c3db3cb5af2bec7f7dac2e37296bac94ddb3e8f915d"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}