{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8fab8521",
   "metadata": {},
   "source": [
    "======================== Import Packages =========================="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "38beb221",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys,os,pdb,glob\n",
    "import numpy as np\n",
    "from astropy import constants as const\n",
    "from astropy import units as u\n",
    "from astropy.table import Table, join, MaskedColumn\n",
    "from astroquery.vizier import Vizier\n",
    "import warnings\n",
    "from astropy.logger import AstropyWarning\n",
    "warnings.filterwarnings('ignore', category=AstropyWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da5042d3",
   "metadata": {},
   "source": [
    "===================== Define Functions ==================="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "571ead0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_data(catalog, join_key='Name', join_type='inner'):\n",
    "\n",
    "    \"\"\"\n",
    "    PURPOSE:    Get data from literature with Vizier\n",
    "\n",
    "    INPUT:      catalog = ctalog name on Vizier (str)\n",
    "                join_key = column header to join tables, if multiple (str; optional)\n",
    "                join_type = way to join tables, if multiple (str; optional)\n",
    "\n",
    "    OUTPUT:     t = data table (AstroPy Table)\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    ### GET FULL CATALOG (ALL COLUMNS, ALL ROWS)\n",
    "    viz = Vizier(catalog=catalog, columns=['**'])\n",
    "    viz.ROW_LIMIT = -1\n",
    "    tv = viz.get_catalogs(catalog)\n",
    "\n",
    "    ### IF MULTIPLE TABLES, JOIN THEN\n",
    "    for i, val in enumerate(tv.keys()):\n",
    "        if i == 0:\n",
    "            t = tv[val]\n",
    "        else:\n",
    "            tt = tv[val]\n",
    "            if join_key in tt.columns:\n",
    "                t = join(t, tt, join_type=join_type, keys=join_key)\n",
    "\n",
    "    return t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ebd8e8f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fix_units(b, l, f, d, t):\n",
    "    \n",
    "    \"\"\"\n",
    "    PURPOSE:    Change units of inputs to work with\n",
    "                dust mass calculations\n",
    "\n",
    "    INPUT:      b = dust opacity power-law index (unitless; float)\n",
    "                l = wavelength of observations (microns; float)\n",
    "                f = observed flux (mJy; float)\n",
    "                d = istance to disk (pc; float)\n",
    "                t = isk temperature (K; float)\n",
    "\n",
    "    OUTPUT:     Inputs but now with correct units\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    b = float(b) \n",
    "    t = float(t) * u.K \n",
    "    d = float(d) * u.pc.to(u.cm) * u.cm\n",
    "    l = float(l) * u.um.to(u.cm) * u.cm \n",
    "    f = float(f) * u.mJy\n",
    "\n",
    "    return b, l, f, d, t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7eccf9a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "05b3faa0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_planck(l, t):\n",
    "\n",
    "    \"\"\"\n",
    "    PURPOSE:    Calculate Planck function\n",
    "\n",
    "    INPUT:      l = wavelength of observations in cm (float)\n",
    "                t = disk temperature in K (float)\n",
    "\n",
    "    OUTPUT:     planck = Planck function (float)\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    ### CONVERT WAVELENGTH TO FREQUENCY\n",
    "    nu = (const.c.cgs / l)\n",
    "\n",
    "    ### CALCULATE PLANCK FUNCTION\n",
    "    plank = 2 * const.h.cgs * nu**3 / const.c.cgs**2 / (np.exp(const.h.cgs *nu / (const.k_B.cgs * t)) -1)\n",
    "\n",
    "    return plank"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "063bf06b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_kappa(l, b):\n",
    "\n",
    "    \"\"\"\n",
    "    PURPOSE:    Calculate kappa parameter for dust mass calculation\n",
    "\n",
    "    INPUT:      l = wavelength of observations in cm (float)\n",
    "                b = dust opacity power-law index (float)\n",
    "\n",
    "    OUTPUT:     kappa = kappa parameter for dust mass calculation (float)\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    ### CONVERT WAVELENGTH TO FREQUENCY\n",
    "    nu = (const.c.cgs / l)\n",
    "\n",
    "    ### ASSUMES NORMALIZATION OF 2.3 cm^2/g at 230 GHz\n",
    "    nf = 2.3 / ( (230. / 1000.)**b )\n",
    "\n",
    "    ### CALCULATE KAPPA PARAMETER\n",
    "    kappa = (nf * (nu / (1e9 * 1000 * u.Hz))**b ) * ((u.cm**2)/(u.g))\n",
    "        \n",
    "    return kappa "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "451718d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_mult(d, kappa, b_nu):\n",
    "\n",
    "    \"\"\"\n",
    "    PURPOSE:    Calculate multiplication factor that produces dust mass\n",
    "                when applied to source flux in mJy\n",
    "                (matches Eq. 1 in Ansdell+2016 when d=150pc)\n",
    "\n",
    "    INPUT:      d = distance to source in cm (float)\n",
    "                kappa = kappa parameter from get_kappa() (float)\n",
    "                b_nu = Planck function from get_planck() (float)\n",
    "\n",
    "    OUTPUT:     mult = mltiplication factor that produces dust mass\n",
    "                       when applied to source flux in mJy (float)\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    mult = (1e-26) * (d**2) / (kappa * b_nu) / const.M_sun.cgs\n",
    "        \n",
    "    return mult"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "de7c3fb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_dustmass(b, l, f, d, t):\n",
    "    \n",
    "    \"\"\"\n",
    "    PURPOSE:    Calculate disk dust mass using prescription from\n",
    "                Hildebrand 1983 (1983QJRAS..24..267H)\n",
    "\n",
    "    INPUT:      b = dust opacity power-law index (unitless; float)\n",
    "                l = wavelength of observations (microns; float)\n",
    "                f = observed flux (mJy; float)\n",
    "                d = distance to disk (pc; float)\n",
    "                t = disk temperature (K; float)\n",
    "\n",
    "    OUTPUT:     Disk dust mass in Earth masses (float)\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    b, l, f, d, t = fix_units(b, l, f, d, t)\n",
    "    k_nu = get_kappa(l, b)\n",
    "    b_nu = get_planck(l, t)\n",
    "    mult = get_mult(d, k_nu, b_nu)\n",
    "    mdust = mult * f * (const.M_sun.cgs / const.M_earth.cgs)\n",
    "\n",
    "    return round(mdust.value, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be496c78",
   "metadata": {},
   "source": [
    "============================= Code =================================="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a3332e0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#### LOAD IN LUPUS DATA\n",
    "T = get_data(\"J/ApJ/828/46\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b1ae113a",
   "metadata": {},
   "outputs": [],
   "source": [
    "### GET DUST MASSES\n",
    "DM, e_DM = [], []\n",
    "for i, val in enumerate(T['Name']):\n",
    "\n",
    "    ### SET INPUTS TO DUST MASS CALCULATION\n",
    "    Beta_Index = 1.                   # DUST OPACITY POWER-LAW INDEX (UNITLESS)\n",
    "    Lambda_Obs = 890.                 # WAVELENGTH OF OBSERVATIONS (MICRONS)\n",
    "    Temp_Disk = 20.                   # DISK TEMPERATURE (KELVIN)\n",
    "    Dist_Disk = T['Dist'][i]           # DISTANCE TO DISK (PC)\n",
    "    Flux_Disk = T['F890'][i]         # OBSERVED FLUX (mJy)\n",
    "    e_Flux_Disk = T['e_F890'][i]     # ERROR ON OBSERVED FLUX (mJy)\n",
    "\n",
    "    ### CALCULATE DUST MASS\n",
    "    DM.append(get_dustmass(Beta_Index, Lambda_Obs, Flux_Disk, Dist_Disk, Temp_Disk))\n",
    "    e_DM.append(get_dustmass(Beta_Index, Lambda_Obs, e_Flux_Disk, Dist_Disk, Temp_Disk))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "ced753fd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><i>Table length=89</i>\n",
       "<table id=\"table4722590576\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>Name</th><th>MDust</th><th>e_MDust</th></tr></thead>\n",
       "<thead><tr><th>str17</th><th>float64</th><th>float64</th></tr></thead>\n",
       "<tr><td>J15430131-3409153</td><td>0.00235</td><td>0.07288</td></tr>\n",
       "<tr><td>J15430227-3444059</td><td>0.05172</td><td>0.06347</td></tr>\n",
       "<tr><td>J15445789-3423392</td><td>-0.01175</td><td>0.04232</td></tr>\n",
       "<tr><td>J15450634-3417378</td><td>3.52637</td><td>0.09404</td></tr>\n",
       "<tr><td>J15450887-3417333</td><td>10.87769</td><td>0.11755</td></tr>\n",
       "<tr><td>J15592523-4235066</td><td>-0.00705</td><td>0.04467</td></tr>\n",
       "<tr><td>J16000060-4221567</td><td>0.56422</td><td>0.04467</td></tr>\n",
       "<tr><td>J16000236-4222145</td><td>28.17573</td><td>0.14811</td></tr>\n",
       "<tr><td>J16002612-4153553</td><td>0.28211</td><td>0.04467</td></tr>\n",
       "<tr><td>J16011549-4152351</td><td>19.29397</td><td>0.20923</td></tr>\n",
       "<tr><td>...</td><td>...</td><td>...</td></tr>\n",
       "<tr><td>Sz 88A</td><td>3.73639</td><td>0.12538</td></tr>\n",
       "<tr><td>Sz 88B</td><td>-0.08359</td><td>0.07941</td></tr>\n",
       "<tr><td>Sz 90</td><td>9.12365</td><td>0.19225</td></tr>\n",
       "<tr><td>Sz 95</td><td>1.7052</td><td>0.07523</td></tr>\n",
       "<tr><td>Sz 96</td><td>1.7052</td><td>0.11702</td></tr>\n",
       "<tr><td>Sz 97</td><td>1.93924</td><td>0.07523</td></tr>\n",
       "<tr><td>Sz 98</td><td>99.17315</td><td>0.59348</td></tr>\n",
       "<tr><td>Sz 99</td><td>-0.01672</td><td>0.07523</td></tr>\n",
       "<tr><td>V1192 Sco</td><td>0.38033</td><td>0.07941</td></tr>\n",
       "<tr><td>V856 Sco</td><td>23.32527</td><td>0.11702</td></tr>\n",
       "</table></div>"
      ],
      "text/plain": [
       "<Table length=89>\n",
       "       Name        MDust   e_MDust\n",
       "      str17       float64  float64\n",
       "----------------- -------- -------\n",
       "J15430131-3409153  0.00235 0.07288\n",
       "J15430227-3444059  0.05172 0.06347\n",
       "J15445789-3423392 -0.01175 0.04232\n",
       "J15450634-3417378  3.52637 0.09404\n",
       "J15450887-3417333 10.87769 0.11755\n",
       "J15592523-4235066 -0.00705 0.04467\n",
       "J16000060-4221567  0.56422 0.04467\n",
       "J16000236-4222145 28.17573 0.14811\n",
       "J16002612-4153553  0.28211 0.04467\n",
       "J16011549-4152351 19.29397 0.20923\n",
       "              ...      ...     ...\n",
       "           Sz 88A  3.73639 0.12538\n",
       "           Sz 88B -0.08359 0.07941\n",
       "            Sz 90  9.12365 0.19225\n",
       "            Sz 95   1.7052 0.07523\n",
       "            Sz 96   1.7052 0.11702\n",
       "            Sz 97  1.93924 0.07523\n",
       "            Sz 98 99.17315 0.59348\n",
       "            Sz 99 -0.01672 0.07523\n",
       "        V1192 Sco  0.38033 0.07941\n",
       "         V856 Sco 23.32527 0.11702"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "TD = Table()\n",
    "TD['Name'] = np.copy(T['Name'])\n",
    "TD.add_column(MaskedColumn(name='MDust', data=DM))\n",
    "TD.add_column(MaskedColumn(name='e_MDust', data=e_DM))\n",
    "display(TD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77e4909b",
   "metadata": {},
   "outputs": [],
   "source": [
    "### UNCOMMENT TO WRITE FILE \n",
    "# TD.write('../output/dustmasses.txt', format='ascii.ipac') "
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "main_language": "python",
   "notebook_metadata_filter": "-all"
  },
  "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}