First version

black_box.py

@@ -0,0 +1,152 @@
+"""
+
+ariel@iag
+04/02/2020
+
+"Hidden functions" to use in LAPIS classes
+
+"""
+
+import bagpipes as pipes
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+def make_bagpipes_model(age, tau, mass, metallicity, Av, z):
+    """
+
+    Plots a BAGPIPES model for S-PLUS filters
+
+    """
+
+    exp = {}
+    exp["age"] = age
+    exp["tau"] = tau
+    exp["massformed"] = mass
+    exp["metallicity"] = metallicity
+
+    dust = {}
+    dust["type"] = "Calzetti"
+    dust["Av"] = Av
+    dust["eta"] = 2.
+
+    nebular = {}
+    nebular["logU"] = -2.
+
+    model_components = {}
+    model_components["redshift"] = z
+    model_components["exponential"] = exp
+    model_components["dust"] = dust
+    model_components["t_bc"] = 0.01
+    model_components["nebular"] = nebular
+
+    filter_list_splus = ['./filters/F378.dat',
+                         './filters/F395.dat',
+                         './filters/F410.dat',
+                         './filters/F430.dat',
+                         './filters/F515.dat',
+                         './filters/F660.dat',
+                         './filters/F861.dat',
+                         './filters/u_SPLUS.dat',
+                         './filters/g_SPLUS.dat',
+                         './filters/r_SPLUS.dat',
+                         './filters/i_SPLUS.dat',
+                         './filters/z_SPLUS.dat']
+
+    model = pipes.model_galaxy(model_components, filt_list=filter_list_splus)
+
+    return model
+
+
+def read_splus_photometry(galaxy_table, galaxy_index):
+    """
+
+    Reads S-PLUS photometry
+
+    """
+
+    magnitudes = np.array([galaxy_table['F378_auto'][galaxy_index],
+                           galaxy_table['F395_auto'][galaxy_index],
+                           galaxy_table['F410_auto'][galaxy_index],
+                           galaxy_table['F430_auto'][galaxy_index],
+                           galaxy_table['F515_auto'][galaxy_index],
+                           galaxy_table['F660_auto'][galaxy_index],
+                           galaxy_table['F861_auto'][galaxy_index],
+                           galaxy_table['uJAVA_auto'][galaxy_index],
+                           galaxy_table['g_auto'][galaxy_index],
+                           galaxy_table['r_auto'][galaxy_index],
+                           galaxy_table['i_auto'][galaxy_index],
+                           galaxy_table['z_auto'][galaxy_index]])
+
+    magnitude_errors = np.array([galaxy_table['eF378_auto'][galaxy_index],
+                                 galaxy_table['eF395_auto'][galaxy_index],
+                                 galaxy_table['eF410_auto'][galaxy_index],
+                                 galaxy_table['eF430_auto'][galaxy_index],
+                                 galaxy_table['eF515_auto'][galaxy_index],
+                                 galaxy_table['eF660_auto'][galaxy_index],
+                                 galaxy_table['eF861_auto'][galaxy_index],
+                                 galaxy_table['euJAVA_auto'][galaxy_index],
+                                 galaxy_table['eg_auto'][galaxy_index],
+                                 galaxy_table['er_auto'][galaxy_index],
+                                 galaxy_table['ei_auto'][galaxy_index],
+                                 galaxy_table['ez_auto'][galaxy_index]])
+
+    wls = np.array([3772.55925537, 3940.44301183, 4094.05290923, 4291.28342875,
+                    5132.43006063, 6613.57355705, 8605.51162257, 3527.31414,
+                    4716.26251124, 6222.43591412, 7644.41843457, 8912.95141873])
+
+    f_nu = 3631 * 10 ** (-0.4 * magnitudes)
+    f_lamb = (1 / (3.34e4 * (wls ** 2))) * f_nu
+
+    f_error = (magnitude_errors / 1.086) * f_lamb
+
+    return wls, f_lamb, f_error
+
+
+def plot_model_and_observation(galaxy_table, galaxy_index, age, tau, mass, metallicity, Av):
+    """
+
+    Overplots BAGPIPES model and S-PLUS observation
+
+    """
+
+    model = make_bagpipes_model(age=age, tau=tau, mass=mass, metallicity=metallicity, Av=Av,
+                                z=galaxy_table['z'][galaxy_index])
+
+    wl, f_lamb, f_error = read_splus_photometry(galaxy_table, galaxy_index)
+
+    z = galaxy_table['z'][galaxy_index]
+
+    wave_limit = (model.wavelengths*(1 + z) > 3000) & (model.wavelengths*(1 + z) < 10000)
+
+    model.sfh.plot()
+
+    plt.figure(figsize=(13, 7))
+    plt.plot(model.wavelengths[wave_limit]*(1 + z), model.spectrum_full[wave_limit]/(1 + z), color='darkcyan',
+             label='Model Spectrum')
+    plt.plot(model.filter_set.eff_wavs, model.photometry/(1 + z), 'o', ms=10, color='navy', label='Model Photometry')
+    plt.errorbar(wl, f_lamb, yerr=f_error, fmt='o', ms=10, color='forestgreen', label='Observed Photometry')
+
+    plt.legend(frameon=False, fontsize=14)
+
+    plt.xlabel(r'$\lambda \mathrm{[\AA]}$', fontsize=16)
+    plt.ylabel(r'$F_\lambda \mathrm{[erg \, cm^2 \, s^{-1} \, \AA^{-1}]}$', fontsize=16)
+
+    return model
+
+
+if __name__ == '__main__':
+    from astropy.table import Table
+
+    catalog = Table.read('splus_laplata.txt', format='ascii')
+
+    # Select only galaxies with no missing bands and r magnitude < 17
+    catalog = catalog[(catalog['nDet_auto'] == 12) & (catalog['r_auto'] < 17) & (catalog['class_2'] == 'GALAXY')]
+
+    # ETG : 10
+    # Star-forming :
+
+    model_galaxy = plot_model_and_observation(galaxy_table=catalog, galaxy_index=10, age=6, tau=2, mass=11,
+                                              metallicity=1, Av=1)
+    plt.show()
+

fun_with_bagpipes.ipynb

@@ -0,0 +1,237 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Making models with Bagpipes\n",
+    "\n",
+    "Let's import bagipes:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import bagpipes as pipes\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Creating a model galaxy"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# The star-formation history:\n",
+    "\n",
+    "exp = {}                          \n",
+    "exp[\"age\"] = 6.                   # Age in Gyr\n",
+    "exp[\"tau\"] = 1.                   # Rate at which SFR decreases\n",
+    "exp[\"massformed\"] = 9.            # Mass, log_10(M*/M_solar)\n",
+    "exp[\"metallicity\"] = 0.5          # Metallicity, Z/Z_oldsolar\n",
+    "\n",
+    "#Dust:\n",
+    "\n",
+    "dust = {}                         \n",
+    "dust[\"type\"] = \"Calzetti\"         # Define the shape of the attenuation curve\n",
+    "dust[\"Av\"] = 0.2                  # V-band attenuation\n",
+    "\n",
+    "\n",
+    "# Compiling the parameters in the model_components dictionary:\n",
+    "\n",
+    "model_components = {}                   \n",
+    "model_components[\"redshift\"] = 0.0      # Observed redshift  \n",
+    "model_components[\"exponential\"] = exp   # The star-formation history dictionary\n",
+    "model_components[\"dust\"] = dust         # The dust component\n",
+    "\n",
+    "# Wavelength array to plot spectra.\n",
+    "wl = np.arange(1000, 10000, 1)\n",
+    "\n",
+    "# Creating a model galaxy\n",
+    "model = pipes.model_galaxy(model_components, spec_wavs=wl)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Let's see what we have done"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig = model.sfh.plot()\n",
+    "fig = model.plot()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Adding nebular component\n",
+    "\n",
+    "Will add emission lines to the spectra, equivalent to adding gas to the galaxy."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Add a nebular component\n",
+    "\n",
+    "nebular = {}                      \n",
+    "nebular[\"logU\"] = -2            # log_10(ionization parameter)\n",
+    "\n",
+    "# The young stars that are ionizing the gas are also more attenuated by dust\n",
+    "\n",
+    "dust[\"eta\"] = 2.                        # Extra dust for young stars: multiplies Av\n",
+    "model_components[\"t_bc\"] = 0.01         # Lifetime of birth clouds (Gyr)\n",
+    "\n",
+    "model_components[\"nebular\"] = nebular\n",
+    "\n",
+    "model = pipes.model_galaxy(model_components, spec_wavs=wl)\n",
+    "\n",
+    "fig = model.plot()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Let's plot the entire model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "model.plot_full_spectrum()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Playing with S-PLUS data\n",
+    "\n",
+    "First, lets load our data and select only bright galaxies."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from astropy.table import Table\n",
+    "\n",
+    "catalog = Table.read('splus_laplata.txt', format='ascii')\n",
+    "\n",
+    "# Select only galaxies with no missing bands and r magnitude < 17\n",
+    "catalog = catalog[(catalog['nDet_auto'] == 12) & (catalog['r_auto'] < 17) & (catalog['class_2'] == 'GALAXY')]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The file black_box.py contains \"hidden functions\" that we will use in this example.\n",
+    "\n",
+    "The function plot_model_and_observation will plot the S-PLUS fluxes of a galaxy in our sample along with a model generated by bagpipes. \n",
+    "\n",
+    "Let's import the function."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from black_box import plot_model_and_observation"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For example, here we are plotting photometry from galaxy 869 in our sample with a bagpipes model:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot_model_and_observation(galaxy_table=catalog, galaxy_index=869, age=4, tau=2.5, mass=10.4\n",
+    "                           , metallicity=1, Av=0.6)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Let's see if you can find a good model for galaxies 46 and 818!\n",
+    "\n",
+    "### Galaxy 46:\n",
+    "RA = 53.5671, Dec = -1.1914\n",
+    "\n",
+    "<img src=\"files/46.png\">\n",
+    "\n",
+    "### Galaxy 818:\n",
+    "RA = 37.9299, Dec = 0.9044\n",
+    "\n",
+    "<img src=\"files/818.png\">"
+   ]
+  },
+  {
+   "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.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}