{ "cells": [ { "cell_type": "markdown", "id": "0ea9f7d1", "metadata": {}, "source": [ "======================== Import Packages ==========================" ] }, { "cell_type": "code", "execution_count": 1, "id": "ade79182", "metadata": {}, "outputs": [], "source": [ "import os, sys, pdb, glob\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from astropy.table import Table, join, MaskedColumn\n", "import matplotlib as mpl\n", "from astropy import constants as const\n", "from astropy import units as u\n", "from scipy.stats import spearmanr, linregress\n", "from matplotlib.ticker import MaxNLocator\n", "import matplotlib.patches as mpatches\n", "from matplotlib.patches import FancyBboxPatch\n", "from astroquery.vizier import Vizier\n", "import warnings\n", "from astropy.logger import AstropyWarning\n", "from IPython.display import display\n", "warnings.filterwarnings('ignore', category=AstropyWarning)" ] }, { "cell_type": "markdown", "id": "82533b69", "metadata": {}, "source": [ "========================== Define Functions ==========================" ] }, { "cell_type": "code", "execution_count": 2, "id": "940c5107", "metadata": {}, "outputs": [], "source": [ "def get_model_grid_old(f):\n", "\n", " ### LOAD FILE FROM WILLIAMS & BEST 2014 (2014ApJ...788...59W)\n", " gasgrid = Table.read(f, format='ascii.csv')\n", "\n", " ### LOAD RELEVANT MODEL OUTPUTS\n", " gasgrid['M_gas'].name = 'Mgas'\n", " gasgrid['f_3-2_13co'].name = 'F13CO32'\n", " gasgrid['f_3-2_c18o'].name = 'FC18O32'\n", " gasgrid['f_3-2_c18o_low'].name = 'FC18O32l'\n", "\n", " ### REMOVE MASKED VALUES (1 BAD END-OF-LINE IN GRID FROM JPW)\n", " gasgrid = gasgrid[np.where(~gasgrid['F13CO32'].mask)]\n", "\n", " return gasgrid" ] }, { "cell_type": "code", "execution_count": 3, "id": "c90b6798", "metadata": {}, "outputs": [], "source": [ "def get_model_grid(f):\n", "\n", " \"\"\"\n", " PURPOSE: Load gas model grid\n", " Remove unnecessary model grid points\n", "\n", " INPUT: Table 3 from Williams & Best 2014 (2014ApJ...788...59W)\n", "\n", " OUTPUT: Gas model grid\n", "\n", " \"\"\"\n", "\n", " ### LOAD FILE FROM WILLIAMS & BEST 2014 (2014ApJ...788...59W)\n", " gg = Table.read(f, format='ascii.cds')\n", "\n", " ### ONLY KEEP RELEVANT MODEL OUTPUTS TO REDUCE SIZE FOR FITTING\n", " gg = gg[np.where(gg['gamma'] == 0.)]\n", " gg = gg[np.where(gg['Mgas'] <= 0.03)]\n", "\n", " return gg" ] }, { "cell_type": "code", "execution_count": 4, "id": "fdf5d643", "metadata": {}, "outputs": [], "source": [ "def get_cal_err(f, e, e_cal=0.10):\n", "\n", " \"\"\"\n", " PURPOSE: Calculate multiplication factor for adding \n", " calibration error to flux measurement error\n", " (assumed to be 10% unless otherwise specified)\n", "\n", " INPUT: f = measured flux (float)\n", " e = measurement error (float)\n", " e_cal = calibration error fraction (float; optional)\n", "\n", " OUTPUT: Multiplication factor for adding calibration error\n", "\n", " \"\"\"\n", "\n", " mult = np.sqrt((f * e_cal)**2 + (e)**2)\n", "\n", " return mult" ] }, { "cell_type": "code", "execution_count": 5, "id": "9a61742e", "metadata": {}, "outputs": [], "source": [ "def get_model_idx(g13, g18_ism, g18_lo, f13, e13, f18, e18, d13, d18):\n", "\n", " \"\"\"\n", " PURPOSE: Index model grid points for a given flux measurement \n", "\n", " INPUT: g13 = model grid points for 13CO flux\n", " g18_hi = model grid points for C18O flux (ISM abundance)\n", " g18_lo = model grid points for C18O flux (low abundance)\n", " f13 = 13CO flux measurement (float)\n", " e13 = 13CO flux measurement error(float)\n", " f18 = C18O flux measurement (float)\n", " e18 = C18O flux measurement error(float)\n", " d13 = detection flags for 13CO flux (masked array)\n", " d18 = detection flags for C18O flux (masked array)\n", "\n", " OUTPUT: ifit = indexes of model grid\n", " inote = note indicating type of detection\n", "\n", " \"\"\"\n", "\n", " ### IF BOTH LINES DETECTED DETECTED\n", " ### INDEX GRID WITHIN ERRORS (MEASUREMENT + FLUX CAL)\n", " if (d13 == True) and (d18 == True):\n", " i13 = ( abs(g13 - f13) < get_cal_err(f13, e13) )\n", " i18_ism = ( abs(g18_ism - f18) < get_cal_err(f18, e18) )\n", " i18_lo = ( abs(g18_lo - f18) < get_cal_err(f18, e18) )\n", " inote = \"GD\"\n", "\n", " ### IF ONLY 13CO DETECTED\n", " ### INDEX GRID WITHIN UPPER LIMIT FOR C180\n", " elif (d13 == True) and (d18 == False):\n", " i13 = (abs(g13 - f13) < get_cal_err(f13, e13) )\n", " i18_ism = (g18_ism < f18)\n", " i18_lo = (g18_lo < f18)\n", " inote = \"D13\"\n", "\n", " ### IF BOTH LINES UNDETECTED\n", " ### INDEX GRID WITHIN UPPER LIMITS FOR BOTH\n", " elif (d13 == False) and (d18 == False):\n", " i13 = (g13 < f13)\n", " i18_ism = (g18_ism < f18)\n", " i18_lo = (g18_lo < f18)\n", " inote = \"ND\"\n", "\n", " ### STEP CODE IF UNKNOWN RESULT\n", " else:\n", " print(\"Unknown result\")\n", " pdb.set_trace()\n", "\n", " ### COMBINING ISM & LOW C18O ABUNDANCES\n", " ### CONSERVATIVE SINCE CONSIDERING \"ALL\" POSSIBLE CO ABUNDANCES\n", " i18 = i18_ism + i18_lo\n", "\n", " ### KEEP ONLY THOSE IN BOTH LINES\n", " ### I.E., CREATING BOX AROUND MEASUREMENT OR UPPER LIMIT\n", " ifit = i13 & i18\n", "\n", " return ifit, inote" ] }, { "cell_type": "code", "execution_count": 6, "id": "2c554ab3", "metadata": {}, "outputs": [], "source": [ "def get_model_gasmass(gm, ifit, inote):\n", "\n", " \"\"\"\n", " PURPOSE: Get gas mass from model grid \n", "\n", " INPUT: gm = all gas masses from model grid \n", " ifit = model grid indexes for a given flux measurement\n", " inote = note indicating type of detection\n", "\n", " OUTPUT: mgas_fit = gas mass based on model fit\n", " mgas_min = lower limti of gas mass based on model fit\n", " mgas_max = upper limit of gas mass based on model fit\n", " mgas_note = note indicating type of gas mass estimate\n", " \"\"\"\n", "\n", " nfit = np.count_nonzero(ifit)\n", " if (nfit > 0):\n", " \n", " mgasfit = gm[ifit]\n", "\n", " ### ONLY 13CO DETECTED\n", " if (d13==True) and (d18==True):\n", " mgas_fit = 10**np.mean(np.log10(mgasfit))\n", " mgas_min, mgas_max = np.min(mgasfit), np.max(mgasfit)\n", " mgas_note = \"GF\"\n", "\n", " ### BOTH LINES DETECTED\n", " elif (d13==True) and (d18==False):\n", " mgas_fit = 10**np.mean(np.log10(mgasfit))\n", " mgas_min, mgas_max = 0.0, np.max(mgasfit)\n", " mgas_note = \"GL\"\n", "\n", " ### BOTH LINES UNDETECTED\n", " elif (d13==False) and (d18==False):\n", " mgas_fit = -99.0\n", " mgas_min, mgas_max = 0.0, np.max(mgasfit)\n", " mgas_note = \"UL\"\n", "\n", " ### STOP CODE IF UNKNOWN RESULT\n", " else:\n", " print(\"Unknown flux measurement result\")\n", " pdb.set_trace()\n", " \n", " else:\n", "\n", " ### STOP CODE IF NO MATCHES TO MODEL GRID\n", " print(\"No matches to model grid\")\n", " pdb.set_trace()\n", "\n", " ### COMBINE NOTES\n", " mgas_note = (', ').join([inote, mgas_note])\n", " \n", " return mgas_fit, mgas_min, mgas_max, mgas_note" ] }, { "cell_type": "code", "execution_count": 7, "id": "d7efacac", "metadata": {}, "outputs": [], "source": [ "def get_gasmass(g, f13, e13, d13, f18, e18, d18):\n", "\n", " \"\"\"\n", " PURPOSE: Calculate gas mass\n", "\n", " INPUT: g = model grid \n", " f13 = 13CO flux measurement (float)\n", " e13 = 13CO flux measurement error(float)\n", " d13 = detection flags for 13CO flux (masked array)\n", " f18 = C18O flux measurement (float)\n", " e18 = C18O flux measurement error(float)\n", " d18 = detection flags for C18O flux (masked array)\n", "\n", " OUTPUT: mg_fit = gas mass based on model fit\n", " mgas_min = lower limti of gas mass based on model fit\n", " mgas_max = upper limit of gas mass based on model fit\n", " mgas_note = note indicating type of gas mass estimate\n", " \"\"\"\n", "\n", " ### INDEX MODEL GRID FOR THIS FLUX MEASUREMENT\n", " i_fit, i_note = get_model_idx(g['F13CO32'], g['FC18O32'], g['FC18O32l'], \n", " f13, e13, f18, e18, d13, d18)\n", "\n", " ### CALCULATE GAS MASS\n", " mgas_fit, mgas_lo, mgas_hi, mgas_note = get_model_gasmass(g['Mgas'], i_fit, i_note)\n", " \n", " return mgas_fit, mgas_lo, mgas_hi, mgas_note" ] }, { "cell_type": "code", "execution_count": 8, "id": "43ff2a7f", "metadata": {}, "outputs": [], "source": [ "def get_data(catalog, join_key='Name'):\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", "\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='inner', keys=join_key)\n", "\n", " return t" ] }, { "cell_type": "markdown", "id": "0285bb8c", "metadata": {}, "source": [ "============================= Code ==================================" ] }, { "cell_type": "code", "execution_count": 9, "id": "9f38c46c", "metadata": {}, "outputs": [], "source": [ "### GET LUPUS DATA\n", "T = get_data(\"J/ApJ/828/46\")" ] }, { "cell_type": "code", "execution_count": 10, "id": "2d855ca1", "metadata": {}, "outputs": [], "source": [ "### LOAD GAS MODEL GRID FROM WILLIAMS & BEST 2014 (2014ApJ...788...59W)\n", "# G = get_model_grid('../input/apj495435t3_mrt.txt')\n", "G = get_model_grid_old('../data/gasgrid.csv')" ] }, { "cell_type": "code", "execution_count": 11, "id": "e189600b", "metadata": {}, "outputs": [], "source": [ "### GET GAS MASSES\n", "mg_f, mg_m, mg_l, mg_h, mg_n = [], [], [], [], []\n", "for i, val in enumerate(T['Name']):\n", "\n", " ### GET GAS FLUXES IN JY, SCALED TO 140 PC TO MATCH MODEL GRID\n", " f2l = (T['Dist'][i] / 140.)**2 / 1000.0\n", " f13, e13 = f2l * T['F13CO'][i], f2l * T['e_F13CO'][i]\n", " f18, e18 = f2l * T['F18CO'][i], f2l * T['e_F18CO'][i]\n", "\n", " ### FLAG (NON-)DETECTIONS \n", " d13 = T['l_F13CO'][i] != '<'\n", " d18 = T['l_F18CO'][i] != '<'\n", "\n", " ### CALCULATE GAS MASSES\n", " mgf, mgl, mgh, mgn = get_gasmass(G, f13, e13, d13, f18, e18, d18)\n", "\n", " ### SAVE GAS MASSES (M_JUP), LIMITS, AND FLAGS\n", " ### WEIRD ROUNDING / FORMATTING IS TO MATCH OUTPUT IN PAPER\n", " if mgn == 'ND, UL':\n", " mg_n.append('<')\n", " mg_f.append(round(float(\"{0:.2e}\".format(mgh))*(const.M_sun.cgs/const.M_jup.cgs).value, 1))\n", " mg_l.append(float(\"{0:.2e}\".format(-99.)))\n", " mg_h.append(float(\"{0:.2e}\".format(-99.)))\n", "\n", " elif mgn == 'D13, GL':\n", " mg_n.append('')\n", " mg_f.append(round(float(\"{0:.2e}\".format(mgf))*(const.M_sun.cgs/const.M_jup.cgs).value, 1))\n", " mg_h.append(round(float(\"{0:.2e}\".format(mgh))*(const.M_sun.cgs/const.M_jup.cgs).value, 1))\n", " mg_l.append(float(\"{0:.2e}\".format(-99.)))\n", "\n", " elif mgn == 'GD, GF':\n", " mg_n.append('')\n", " mg_f.append(round(float(\"{0:.2e}\".format(mgf))*(const.M_sun.cgs/const.M_jup.cgs).value, 1))\n", " mg_l.append(round(float(\"{0:.2e}\".format(mgl))*(const.M_sun.cgs/const.M_jup.cgs).value, 1))\n", " mg_h.append(round(float(\"{0:.2e}\".format(mgh))*(const.M_sun.cgs/const.M_jup.cgs).value, 1))\n", "\n", " else:\n", " print(\"Unknown gas mass result\")\n", " pdb.set_trace()" ] }, { "cell_type": "code", "execution_count": 12, "id": "57f95408", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Table length=89\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
NameMgasl_Mgasb_MgasB_Mgas
str17float64str1float64float64
J15430131-34091531.0<----
J15430227-34440591.0<----
J15445789-34233920.3<----
J15450634-34173780.1--3.1
J15450887-34173333.21.010.5
J15592523-42350660.3<----
J16000060-42215670.3<----
J16000236-42221450.7--1.0
J16002612-41535530.3<----
J16011549-41523516.71.031.4
...............
Sz 88A1.0<----
Sz 88B1.0<----
Sz 900.3--10.5
Sz 950.2--3.1
Sz 961.0<----
Sz 971.0<----
Sz 980.4--10.5
Sz 991.0<----
V1192 Sco1.0<----
V856 Sco1.0<----
" ], "text/plain": [ "\n", " Name Mgas l_Mgas b_Mgas B_Mgas\n", " str17 float64 str1 float64 float64\n", "----------------- ------- ------ ------- -------\n", "J15430131-3409153 1.0 < -- --\n", "J15430227-3444059 1.0 < -- --\n", "J15445789-3423392 0.3 < -- --\n", "J15450634-3417378 0.1 -- 3.1\n", "J15450887-3417333 3.2 1.0 10.5\n", "J15592523-4235066 0.3 < -- --\n", "J16000060-4221567 0.3 < -- --\n", "J16000236-4222145 0.7 -- 1.0\n", "J16002612-4153553 0.3 < -- --\n", "J16011549-4152351 6.7 1.0 31.4\n", " ... ... ... ... ...\n", " Sz 88A 1.0 < -- --\n", " Sz 88B 1.0 < -- --\n", " Sz 90 0.3 -- 10.5\n", " Sz 95 0.2 -- 3.1\n", " Sz 96 1.0 < -- --\n", " Sz 97 1.0 < -- --\n", " Sz 98 0.4 -- 10.5\n", " Sz 99 1.0 < -- --\n", " V1192 Sco 1.0 < -- --\n", " V856 Sco 1.0 < -- --" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "TG = Table()\n", "TG['Name'] = np.copy(T['Name'])\n", "TG.add_column(MaskedColumn(name='Mgas', data=mg_f))\n", "TG.add_column(MaskedColumn(name='l_Mgas', data=mg_n))\n", "TG.add_column(MaskedColumn(name='b_Mgas', data=mg_l, mask=[x==-99.0 for x in mg_l]))\n", "TG.add_column(MaskedColumn(name='B_Mgas', data=mg_h, mask=[x==-99.0 for x in mg_h]))\n", "display (TG)" ] }, { "cell_type": "code", "execution_count": 13, "id": "801cd538", "metadata": {}, "outputs": [], "source": [ "### UNCOMMENT TO WRITE FILE\n", "# TG.write('../output/gasmasses.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 }