color-excess.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "color-excess.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true
    },
    "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.6.7"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5fGL2VYp13Vi"
      },
      "source": [
        "# Analizando o avermelhante interestelar e calculando fotometria sintética"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6UzIN7ds1-G1"
      },
      "source": [
        "## Authors\n",
        "\n",
        "Kristen Larson, Lia Corrales, Stephanie T. Douglas, Kelle Cruz\n",
        "\n",
        "Input from Emir Karamehmetoglu, Pey Lian Lim, Karl Gordon, Kevin Covey"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Btb9Da5P4nLr"
      },
      "source": [
        "## Learning Goals\n",
        "- Investigate extinction curve shapes\n",
        "- Deredden spectral energy distributions and spectra\n",
        "- Calculate photometric extinction and reddening\n",
        "- Calculate synthetic photometry for a dust-reddened star by combining `dust_extinction` and `synphot`\n",
        "- Convert from frequency to wavelength with `astropy.unit` equivalencies\n",
        "- Unit support for plotting with `astropy.visualization`\n",
        "\n",
        "\n",
        "## Keywords\n",
        "dust extinction, synphot, astroquery, units, photometry, extinction, physics, observational astronomy\n",
        "\n",
        "## Companion Content\n",
        "\n",
        "* [Bessell & Murphy (2012)](https://ui.adsabs.harvard.edu/#abs/2012PASP..124..140B/abstract)\n",
        "\n",
        "\n",
        "\n",
        "## Summary\n",
        "\n",
        "In this tutorial, we will look at some extinction curves from the literature, use one of those curves to deredden an observed spectrum, and practice invoking a background source flux in order to calculate magnitudes from an extinction model.\n",
        "\n",
        "The primary libraries we'll be using are [dust_extinction](https://dust-extinction.readthedocs.io/en/latest/) and [synphot](https://synphot.readthedocs.io/en/latest/), which are [Astropy affiliated packages](https://www.astropy.org/affiliated/). \n",
        "\n",
        "We recommend installing the two packages in this fashion:\n",
        "```\n",
        "pip install synphot\n",
        "pip install dust_extinction\n",
        "```\n",
        "This tutorial requires v0.7 or later of `dust_extinction`. To ensure that all commands work properly, make sure you have the correct version installed. If you have v0.6 or earlier installed, run the following command to upgrade\n",
        "```\n",
        "pip install dust_extinction --upgrade\n",
        "```"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vFDq1xGXz_t4"
      },
      "source": [
        "import matplotlib.pyplot as plt\n",
        "%matplotlib inline\n",
        "\n",
        "import numpy as np\n",
        "import astropy.units as u\n",
        "from astropy.table import Table\n",
        "from dust_extinction.parameter_averages import CCM89, F99\n",
        "from synphot import units, config\n",
        "from synphot import SourceSpectrum,SpectralElement,Observation,ExtinctionModel1D\n",
        "from synphot.models import BlackBodyNorm1D\n",
        "from synphot.spectrum import BaseUnitlessSpectrum\n",
        "from synphot.reddening import ExtinctionCurve\n",
        "from astroquery.simbad import Simbad\n",
        "from astroquery.mast import Observations\n",
        "import astropy.visualization"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ssir1UjtD_Wn"
      },
      "source": [
        "# Introduction"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1rvwpEs946ju"
      },
      "source": [
        "Dust in the interstellar medium (ISM) extinguishes background starlight.  The wavelength dependence of the extinction is such that short-wavelength light is extinguished more than long-wavelength light, and we call this effect *reddening*.\n",
        "\n",
        "If you're new to extinction, here is a brief introduction to the types of quantities involved.\n",
        "The fractional change to the flux of starlight is \n",
        "$$\n",
        "\\frac{dF_\\lambda}{F_\\lambda} = -\\tau_\\lambda\n",
        "$$\n",
        "\n",
        "where $\\tau$ is the optical depth and depends on wavelength.  Integrating along the line of sight, the resultant flux is an exponential function of optical depth,\n",
        "$$\n",
        "\\tau_\\lambda = -\\ln\\left(\\frac{F_\\lambda}{F_{\\lambda,0}}\\right).\n",
        "$$\n",
        "\n",
        "With an eye to how we define magnitudes, we usually change the base from $e$ to 10,  \n",
        "$$\n",
        "\\tau_\\lambda = -2.303\\log\\left(\\frac{F_\\lambda}{F_{\\lambda,0}}\\right),\n",
        "$$\n",
        "\n",
        "and define an extinction $A_\\lambda = 1.086 \\,\\tau_\\lambda$ so that\n",
        "$$\n",
        "A_\\lambda = -2.5\\log\\left(\\frac{F_\\lambda}{F_{\\lambda,0}}\\right).\n",
        "$$\n",
        "\n",
        "\n",
        "There are two basic take-home messages from this derivation:\n",
        "\n",
        "* Extinction introduces a multiplying factor $10^{-0.4 A_\\lambda}$ to the flux.\n",
        "* Extinction is defined relative to the flux without dust, $F_{\\lambda,0}$.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "K9gwBSy_R_2E"
      },
      "source": [
        "Once astropy and the affiliated packages are installed, we can import from them as needed:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NYgj1w2X7gXc"
      },
      "source": [
        "# Example 1: Investigate Extinction Models"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ifscJWwyEZtS"
      },
      "source": [
        "The `dust_extinction` package provides various models for extinction $A_\\lambda$ normalized to $A_V$.  The shapes of normalized curves are relatively (and perhaps surprisingly) uniform in the Milky Way.  The little variation that exists is often parameterized by the ratio of extinction ($A_V$) to reddening  in the blue-visual ($E_{B-V}$),\n",
        "$$\n",
        "R_V \\equiv \\frac{A_V}{E_{B-V}}\n",
        "$$\n",
        "\n",
        "where $E_{B-V}$ is differential extinction $A_B-A_V$.  In this example, we show the $R_V$-parameterization for the Clayton, Cardelli, & Mathis (1989, CCM) and the Fitzpatrick (1999) models. [More model options are available in the `dust_extinction` documentation.](https://dust-extinction.readthedocs.io/en/latest/dust_extinction/model_flavors.html)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "b1uhXbRl79FR"
      },
      "source": [
        "# Create wavelengths array.\n",
        "wav = np.arange(0.1, 3.0, 0.001)*u.micron\n",
        "\n",
        "for model in [CCM89, F99]:\n",
        "    for R in (2.0,3.0,4.0):\n",
        "        # Initialize the extinction model\n",
        "        ext = model(Rv=R)\n",
        "        plt.plot(1/wav, ext(wav), label=model.name+' R='+str(R))\n",
        "        \n",
        "plt.xlabel('$\\lambda^{-1}$ ($\\mu$m$^{-1}$)')\n",
        "plt.ylabel('A($\\lambda$) / A(V)')\n",
        "plt.legend(loc='best')\n",
        "plt.title('Some Extinction Laws')\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "d-1XRyovX028"
      },
      "source": [
        "Astronomers studying the ISM often display extinction curves against inverse wavelength (wavenumber) to show the ultraviolet variation, as we do here.  Infrared extinction varies much less and approaches zero at long wavelength in the absence of wavelength-independent, or grey, extinction."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TMEWEKJf-iSL"
      },
      "source": [
        "# Example 2: Deredden a Spectrum"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zTeDJ-dI1cQj"
      },
      "source": [
        "Here we deredden (unextinguish) the IUE ultraviolet spectrum and optical photometry of the star $\\rho$ Oph (HD 147933).\n",
        "\n",
        "First, we will use astroquery to fetch the archival [IUE spectrum from MAST](https://archive.stsci.edu/iue/):"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Zdwr_mdZcDeh"
      },
      "source": [
        "obsTable = Observations.query_object(\"HD 147933\",radius=\"1 arcsec\")\n",
        "obsTable_spec=obsTable[obsTable['dataproduct_type']=='spectrum']\n",
        "obsTable_spec.pprint()\n",
        "\n",
        "obsids = obsTable_spec[39]['obsid']\n",
        "dataProductsByID = Observations.get_product_list(obsids)\n",
        "manifest = Observations.download_products(dataProductsByID)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MevN670kcI4B"
      },
      "source": [
        "We read the downloaded files into an astropy table:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "THilzBfVcDb-"
      },
      "source": [
        "t_lwr = Table.read('./mastDownload/IUE/lwr05639/lwr05639mxlo_vo.fits')\n",
        "print(t_lwr)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IiQ9u8QJcKtC"
      },
      "source": [
        "The `.quantity` extension in the next lines will read the Table columns into Quantity vectors.  Quantities keep the units of the Table column attached to the numpy array values."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "P54jimyWcDZa"
      },
      "source": [
        "wav_UV = t_lwr['WAVE'][0,].quantity\n",
        "UVflux = t_lwr['FLUX'][0,].quantity"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WCIgXzhdcUn7"
      },
      "source": [
        "Now, we use astroquery again to fetch photometry from Simbad to go with the IUE spectrum:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lYPVhZrXcDWw"
      },
      "source": [
        "custom_query = Simbad()\n",
        "custom_query.add_votable_fields('fluxdata(U)','fluxdata(B)','fluxdata(V)')\n",
        "phot_table=custom_query.query_object('HD 147933')\n",
        "Umag=phot_table['FLUX_U']\n",
        "Bmag=phot_table['FLUX_B']\n",
        "Vmag=phot_table['FLUX_V']"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OAxOt2AwcVfF"
      },
      "source": [
        "To convert the photometry to flux, we look up some [properties of the photometric passbands](http://ned.ipac.caltech.edu/help/photoband.lst), including the flux of a magnitude zero star through the each passband, also known as the zero-point of the passband."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cKMXywyCcDUU"
      },
      "source": [
        "wav_U = 0.3660 * u.micron \n",
        "zeroflux_U_nu = 1.81E-23 * u.Watt/(u.m*u.m*u.Hz)\n",
        "wav_B = 0.4400 * u.micron\n",
        "zeroflux_B_nu = 4.26E-23 * u.Watt/(u.m*u.m*u.Hz)\n",
        "wav_V = 0.5530 * u.micron\n",
        "zeroflux_V_nu = 3.64E-23 * u.Watt/(u.m*u.m*u.Hz)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3cFcAlSOcjht"
      },
      "source": [
        "The zero-points that we found for the optical passbands are not in the same units as the IUE fluxes.  To make matters worse, the zero-point fluxes are $F_\\nu$ and the IUE fluxes are $F_\\lambda$.  To convert between them, the wavelength is needed.  Fortunately, astropy provides an easy way to make the conversion with *equivalencies*:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pOFMtrNAcDRa"
      },
      "source": [
        "zeroflux_U = zeroflux_U_nu.to(u.erg/u.AA/u.cm/u.cm/u.s, \n",
        "                              equivalencies=u.spectral_density(wav_U))\n",
        "zeroflux_B = zeroflux_B_nu.to(u.erg/u.AA/u.cm/u.cm/u.s, \n",
        "                              equivalencies=u.spectral_density(wav_B))\n",
        "zeroflux_V = zeroflux_V_nu.to(u.erg/u.AA/u.cm/u.cm/u.s, \n",
        "                              equivalencies=u.spectral_density(wav_V))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7XlrfQtjctNG"
      },
      "source": [
        "Now we can convert from photometry to flux using the definition of magnitude:\n",
        "$$\n",
        "F=F_0\\ 10^{-0.4\\, m}\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "p_DaxRzjcDOZ"
      },
      "source": [
        "Uflux = zeroflux_U * 10.**(-0.4*Umag)\n",
        "Bflux = zeroflux_B * 10.**(-0.4*Bmag)\n",
        "Vflux = zeroflux_V * 10.**(-0.4*Vmag)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ODhSujFjc0hM"
      },
      "source": [
        "Using astropy quantities allow us to take advantage of astropy's unit support in plotting. [Calling `astropy.visualization.quantity_support` explicitly turns the feature on.](http://docs.astropy.org/en/stable/units/quantity.html#plotting-quantities)  Then, when quantity objects are passed to matplotlib plotting functions, the axis labels are automatically labeled with the unit of the quantity.  In addition, quantities are converted automatically into the same units when combining multiple plots on the same axes.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "x31EhibTcDFg"
      },
      "source": [
        "astropy.visualization.quantity_support()\n",
        "\n",
        "plt.plot(wav_UV,UVflux,'m',label='UV')\n",
        "plt.plot(wav_V,Vflux,'ko',label='U, B, V')\n",
        "plt.plot(wav_B,Bflux,'ko')\n",
        "plt.plot(wav_U,Uflux,'ko')\n",
        "plt.legend(loc='best')\n",
        "plt.ylim(0,3E-10)\n",
        "plt.title('rho Oph')\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vCfk0D87dDE8"
      },
      "source": [
        "Finally, we initialize the extinction model, choosing values $R_V = 5$ and $E_{B-V} = 0.5$.  This star is famous in the ISM community for having large-$R_V$ dust in the line of sight."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wdp5ERuqcC8_"
      },
      "source": [
        "Rv = 5.0  # Usually around 3, but about 5 for this star.\n",
        "Ebv = 0.5\n",
        "ext = F99(Rv=Rv)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cpN5n_MIdv26"
      },
      "source": [
        "To extinguish (redden) a spectrum, multiply by the `ext.extinguish` function.  To unextinguish (deredden), divide by the same `ext.extinguish`, as we do here:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DyxkYpdxdvfc"
      },
      "source": [
        "plt.semilogy(wav_UV,UVflux,'m',label='UV')\n",
        "plt.semilogy(wav_V,Vflux,'ko',label='U, B, V')\n",
        "plt.semilogy(wav_B,Bflux,'ko')\n",
        "plt.semilogy(wav_U,Uflux,'ko')\n",
        "\n",
        "plt.semilogy(wav_UV,UVflux/ext.extinguish(wav_UV,Ebv=Ebv),'b',\n",
        "             label='dereddened: EBV=0.5, RV=5')\n",
        "plt.semilogy(wav_V,Vflux/ext.extinguish(wav_V,Ebv=Ebv),'ro',\n",
        "             label='dereddened: EBV=0.5, RV=5')\n",
        "plt.semilogy(wav_B,Bflux/ext.extinguish(wav_B,Ebv=Ebv),'ro')\n",
        "plt.semilogy(wav_U,Uflux/ext.extinguish(wav_U,Ebv=Ebv),'ro')\n",
        "\n",
        "plt.legend(loc='best')\n",
        "plt.title('rho Oph')\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "L-yuz7YfFMxm"
      },
      "source": [
        "Notice that, by dereddening the spectrum, the absorption feature at 2175 Angstrom is removed.  This feature can also be seen as the prominent bump in the extinction curves in Example 1.  That we have smoothly removed the 2175 Angstrom feature suggests that the values we chose, $R_V = 5$ and $E_{B-V} = 0.5$, are a reasonable model for the foreground dust.\n",
        "\n",
        "Those experienced with dereddening should notice that that `dust_extinction` returns $A_\\lambda/A_V$, while other routines like the IDL fm_unred procedure often return $A_\\lambda/E_{B-V}$ by default and need to be divided by $R_V$ in order to compare directly with `dust_extinction`."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WqtsQTZbp9nz"
      },
      "source": [
        "# Example 3: Calculate Color Excess with `synphot`"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WwsM9_56FVYr"
      },
      "source": [
        "Calculating broadband *photometric* extinction is harder than it might look at first.  All we have to do is look up $A_\\lambda$ for a particular passband, right?  Under the right conditions, yes.  In general, no.\n",
        "\n",
        "Remember that we have to integrate over a passband to get synthetic photometry,\n",
        "$$\n",
        "A = -2.5\\log\\left(\\frac{\\int W_\\lambda F_{\\lambda,0} 10^{-0.4A_\\lambda} d\\lambda}{\\int W_\\lambda F_{\\lambda,0} d\\lambda} \\right),\n",
        "$$\n",
        "\n",
        "where $W_\\lambda$ is the fraction of incident energy transmitted through a filter.  See the detailed appendix in [Bessell & Murphy (2012)](https://ui.adsabs.harvard.edu/#abs/2012PASP..124..140B/abstract)\n",
        " for an excellent review of the issues and common misunderstandings in synthetic photometry.\n",
        "\n",
        "There is an important point to be made here. The expression above does not simplify any further. Strictly speaking, it is impossible to convert spectral extinction $A_\\lambda$ into a magnitude system without knowing the wavelength dependence of the source's original flux across the filter in question.  As a special case, if we assume that the source flux is constant in the band (i.e. $F_\\lambda = F$), then we can cancel these factors out from the integrals, and extinction in magnitudes becomes the weighted average of the extinction factor across the filter in question. In that special case, $A_\\lambda$ at $\\lambda_{\\rm eff}$ is a good approximation for magnitude extinction.\n",
        "\n",
        "In this example, we will demonstrate the more general calculation of photometric extinction.  We use a blackbody curve for the flux before the dust, apply an extinction curve, and perform synthetic photometry to calculate extinction and reddening in a magnitude system.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "avGGI5fVX2wV"
      },
      "source": [
        "First, let's get the filter transmission curves:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wO-wwoOnp2B3"
      },
      "source": [
        "# Optional, for when the STScI ftp server is not answering:\n",
        "config.conf.vega_file='http://ssb.stsci.edu/cdbs/calspec/alpha_lyr_stis_008.fits'\n",
        "config.conf.johnson_u_file='http://ssb.stsci.edu/cdbs/comp/nonhst/johnson_u_004_syn.fits'\n",
        "config.conf.johnson_b_file='http://ssb.stsci.edu/cdbs/comp/nonhst/johnson_b_004_syn.fits'\n",
        "config.conf.johnson_v_file='http://ssb.stsci.edu/cdbs/comp/nonhst/johnson_v_004_syn.fits'\n",
        "config.conf.johnson_r_file='http://ssb.stsci.edu/cdbs/comp/nonhst/johnson_r_003_syn.fits'\n",
        "config.conf.johnson_i_file='http://ssb.stsci.edu/cdbs/comp/nonhst/johnson_i_003_syn.fits'\n",
        "config.conf.bessel_j_file='http://ssb.stsci.edu/cdbs/comp/nonhst/bessell_j_003_syn.fits'\n",
        "config.conf.bessel_h_file='http://ssb.stsci.edu/cdbs/comp/nonhst/bessell_h_004_syn.fits'\n",
        "config.conf.bessel_k_file='http://ssb.stsci.edu/cdbs/comp/nonhst/bessell_k_003_syn.fits'\n",
        "\n",
        "u_band = SpectralElement.from_filter('johnson_u')\n",
        "b_band = SpectralElement.from_filter('johnson_b')\n",
        "v_band = SpectralElement.from_filter('johnson_v')\n",
        "r_band = SpectralElement.from_filter('johnson_r')\n",
        "i_band = SpectralElement.from_filter('johnson_i')\n",
        "j_band = SpectralElement.from_filter('bessel_j')\n",
        "h_band = SpectralElement.from_filter('bessel_h')\n",
        "k_band = SpectralElement.from_filter('bessel_k')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hekvrcEGpvYd"
      },
      "source": [
        "If you are running this with your own python, see the [synphot documentation](https://synphot.readthedocs.io/en/latest/#installation-and-setup) on how to install your own copy of the necessary files."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vQTSrYpvJeaY"
      },
      "source": [
        "Next, let's make a background flux to which we will apply extinction.  Here we make a 10,000 K blackbody using the model mechanism from within `synphot` and normalize it to $V$ = 10 in the Vega-based magnitude system."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NuPS0Ij2ncC0"
      },
      "source": [
        "# First, create a blackbody at some temperature.\n",
        "sp = SourceSpectrum(BlackBodyNorm1D, temperature=10000)\n",
        "# sp.plot(left=1, right=15000, flux_unit='flam', title='Blackbody')\n",
        "\n",
        "# Get the Vega spectrum as the zero point flux.\n",
        "vega = SourceSpectrum.from_vega()\n",
        "# vega.plot(left=1, right=15000)\n",
        "\n",
        "# Normalize the blackbody to some chosen magnitude, say V = 10.\n",
        "vmag = 10.\n",
        "v_band = SpectralElement.from_filter('johnson_v')\n",
        "sp_norm = sp.normalize(vmag * units.VEGAMAG, v_band, vegaspec=vega)\n",
        "sp_norm.plot(left=1, right=15000, flux_unit='flam', title='Normed Blackbody')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "H4jYc_F8CiHi"
      },
      "source": [
        "Now we initialize the extinction model and choose an extinction of $A_V$ = 2.  To get the `dust_extinction` model working with `synphot`, we create a wavelength array and make a spectral element with the extinction model as a lookup table."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9H4VHd0VCCBG"
      },
      "source": [
        "# Initialize the extinction model and choose the extinction, here Av = 2.\n",
        "ext = CCM89(Rv=3.1)\n",
        "Av = 2.\n",
        "\n",
        "# Create a wavelength array. \n",
        "wav = np.arange(0.1, 3, 0.001)*u.micron\n",
        "\n",
        "# Make the extinction model in synphot using a lookup table.\n",
        "ex = ExtinctionCurve(ExtinctionModel1D, \n",
        "                     points=wav, lookup_table=ext.extinguish(wav, Av=Av))\n",
        "sp_ext = sp_norm*ex\n",
        "sp_ext.plot(left=1, right=15000, flux_unit='flam',\n",
        "            title='Normed Blackbody with Extinction')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "h0dxfniFF-Bf"
      },
      "source": [
        "Synthetic photometry refers to modeling an observation of a star by multiplying the theoretical model for the astronomical flux through a certain filter response function, then integrating."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_1m5IFT4E2v_"
      },
      "source": [
        "# \"Observe\" the star through the filter and integrate to get photometric mag.\n",
        "sp_obs = Observation(sp_ext, v_band)\n",
        "sp_obs_before = Observation(sp_norm, v_band)\n",
        "# sp_obs.plot(left=1, right=15000, flux_unit='flam',\n",
        "#             title='Normed Blackbody with Extinction through V Filter')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vFtJs0PhGLzx"
      },
      "source": [
        "Next, `synphot` performs the integration and computes magnitudes in the Vega system."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ycvVHlOGCRur"
      },
      "source": [
        "sp_stim_before = sp_obs_before.effstim(flux_unit='vegamag', vegaspec=vega)\n",
        "sp_stim = sp_obs.effstim(flux_unit='vegamag', vegaspec=vega)\n",
        "print('before dust, V =', np.round(sp_stim_before,1))\n",
        "print('after dust, V =', np.round(sp_stim,1))\n",
        "\n",
        "# Calculate extinction and compare to our chosen value.\n",
        "Av_calc = sp_stim - sp_stim_before\n",
        "print('$A_V$ = ', np.round(Av_calc,1))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "asE7JUPuHddK"
      },
      "source": [
        "This is a good check for us to do.  We normalized our spectrum to $V$ = 10 mag and added 2 mag of visual extinction, so the synthetic photometry procedure should reproduce these chosen values, and it does.  Now we are ready to find the extinction in other passbands. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "woWiGJKwGI2v"
      },
      "source": [
        "We calculate the new photometry for the rest of the Johnson optical and the Bessell infrared filters. We calculate extinction $A = \\Delta m$ and plot color excess, $E(\\lambda - V) = A_\\lambda - A_V$.  \n",
        "\n",
        "Notice that `synphot` calculates the effective wavelength of the observations for us, which is very useful for plotting the results.  We show reddening with the model extinction curve for comparison in the plot."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0PMDsoC68w4b"
      },
      "source": [
        "bands = [u_band,b_band,v_band,r_band,i_band,j_band,h_band,k_band]\n",
        "\n",
        "for band in bands:\n",
        "    # Calculate photometry with dust:\n",
        "    sp_obs = Observation(sp_ext, band, force='extrap')\n",
        "    obs_effstim = sp_obs.effstim(flux_unit='vegamag', vegaspec=vega)\n",
        "    # Calculate photometry without dust:\n",
        "    sp_obs_i = Observation(sp_norm, band, force='extrap')\n",
        "    obs_i_effstim = sp_obs_i.effstim(flux_unit='vegamag', vegaspec=vega)\n",
        "  \n",
        "    # Extinction = mag with dust - mag without dust\n",
        "    # Color excess = extinction at lambda - extinction at V\n",
        "    color_excess = obs_effstim - obs_i_effstim - Av_calc\n",
        "    plt.plot(sp_obs_i.effective_wavelength(), color_excess,'or')\n",
        "    print(np.round(sp_obs_i.effective_wavelength(),1), ',', \n",
        "          np.round(color_excess,2))\n",
        "\n",
        "# Plot the model extinction curve for comparison \n",
        "plt.plot(wav,Av*ext(wav)-Av,'--k')\n",
        "plt.ylim([-2,2])\n",
        "plt.xlabel('$\\lambda$ (Angstrom)')\n",
        "plt.ylabel('E($\\lambda$-V)')\n",
        "plt.title('Reddening of T=10,000K Background Source with Av=2')\n",
        "plt.show()  "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9B4gdXzn8LtA"
      },
      "source": [
        "## Exercise\n",
        "Try changing the blackbody temperature to something very hot or very cool.  Are the color excess values the same?  Have the effective wavelengths changed?\n",
        "\n",
        "Note that the photometric extinction changes because the filter transmission is not uniform. The observed throughput of the filter depends on the shape of the background source flux."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "F8ngm7OttA6D"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}