diff --git a/46.png b/46.png
new file mode 100644
index 0000000..c4ebf48
Binary files /dev/null and b/46.png differ
diff --git a/818.png b/818.png
new file mode 100644
index 0000000..5423dd5
Binary files /dev/null and b/818.png differ
diff --git a/black_box.py b/black_box.py
new file mode 100644
index 0000000..856b961
--- /dev/null
+++ b/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()
+
diff --git a/fun_with_bagpipes.ipynb b/fun_with_bagpipes.ipynb
new file mode 100644
index 0000000..94ce93c
--- /dev/null
+++ b/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",
+ "![](\"files/46.png\")
\n",
+ "\n",
+ "### Galaxy 818:\n",
+ "RA = 37.9299, Dec = 0.9044\n",
+ "\n",
+ "![](\"files/818.png\")