{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Set 5" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Stellar Spectra, Part II: Chemical Abundances" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the last problem set, we learned how stellar spectra can be used to measure effective temperatures. In this set, we'll practice using stellar spectra to measure chemical abundances.\n", "#### Fill in any parts of the code or text that have ``FILL_ME``. Code will not work if you don't fill in the missing parts." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# An error class to catch the instances of FILL_ME -- don't change this!\n", "class SolutionMissingError(Exception):\n", " def __init__(self, text=None):\n", " Exception.__init__(self, text if text else \"You need to complete the solution for this code to work!\")\n", "def FILL_ME(words=None):\n", " raise SolutionMissingError(words)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### (a) Plot part of the solar spectrum and identify the two strongest absorption lines in this region." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First let's read in data from the solar spectrum (in the file `solspec.txt`). The spectrum is from Delbouille et al. (1972): http://bass2000.obspm.fr/solar_spect.php\n", "\n", "To read in the data, we're going to use the Python `pandas` package. If you've never used it before, it's a super useful package for dealing with data files. You may need to install the package if you don't already have it: https://pandas.pydata.org/\n", "\n", "You can learn more about `pandas` from the documentation here: https://pandas.pydata.org/pandas-docs/stable/" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Import some useful packages\n", "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "# Open the data file\n", "data = pd.read_csv('solspec.txt', delimiter='\\s+')\n", "\n", "# From the spectrum documentation, we see that the data are normalized to 10000, so let's normalize it to 1 instead\n", "flux = data['Flux']/10000. \n", "\n", "# Read in wavelength data\n", "wavelength = data['Wavelength']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's plot our observed spectrum." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(10,5)) # Make the figure bigger \n", "\n", "FILL_ME('Write plotting code here. Make sure to label your axes and include units as needed!')\n", "\n", "# Shade in the parts that are around the two strongest lines\n", "lowerlimit1 = FILL_ME('Lower wavelength limit for the first line')\n", "upperlimit1 = FILL_ME('Upper wavelength limit for the first line')\n", "lowerlimit2 = FILL_ME('Lower wavelength limit for the second line')\n", "upperlimit2 = FILL_ME('Lower wavelength limit for the second line')\n", "\n", "plt.axvspan(lowerlimit1,upperlimit1,color='r',alpha=0.3)\n", "plt.axvspan(lowerlimit2,upperlimit2,color='r',alpha=0.3)\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(1) What are the central wavelengths of the two strongest lines in this part of the spectrum?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lambda1 = FILL_ME('Central wavelength of first line')\n", "lambda2 = FILL_ME('Central wavelength of second line')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(2) Using some other resource (for example, Googling \"strong solar spectral lines\"), can you identify these lines? Which atomic species produces these lines? What electron transitions do they correspond to?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "FILL_ME('^^ Answer the question above! (You can delete this cell afterward.)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### (b) Try fitting one of the lines with two profiles: a Gaussian and Lorentzian. Comment on the fit.\n", "\n", "We're going to use the `scipy.optimize.curve_fit` package to do these non-linear fits.\n", "\n", "First start with a Gaussian profile:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Import the package\n", "from scipy.optimize import curve_fit\n", "\n", "# First define a Gaussian line profile:\n", "def gauss(x, *params):\n", " \n", " amplitude, mean, stddev = params\n", " \n", " gauss = FILL_ME('standard formula for a Gaussian curve')\n", " \n", " # Turn the standard Gaussian into an inverted Gaussian\n", " gauss = 1. - gauss\n", " \n", " return gauss\n", "\n", "# Make some initial guesses for the fitting parameters (A, mu and sigma above)\n", "p0 = [1., lambda1, 1.]\n", "\n", "# Crop the spectrum to be within the wavelength limits\n", "idx = np.where((wavelength > lowerlimit1) & (wavelength < upperlimit1))[0]\n", "\n", "# Do the actual fitting\n", "coeff, _ = curve_fit(gauss, wavelength[idx], flux[idx], p0=p0)\n", "\n", "# Get the fitted curve\n", "gauss_fit = gauss(wavelength[idx], *coeff)\n", "\n", "# Finally, lets get the fitting parameters:\n", "print('Gaussian amplitude = ', coeff[0])\n", "print('Gaussian mean = ', coeff[1])\n", "print('Gaussian standard deviation = ', coeff[2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we'll fit a Lorentzian profile:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Define a Lorentzian line profile:\n", "def lorentz(x, *params):\n", " \n", " amplitude, mean, fwhm = params # Note that the 3rd parameter is \"full-width half maximum,\" NOT standard deviation\n", " \n", " # Standard formula for a Lorentzian curve\n", " lorentz = amplitude*1./(1.+np.power((x-mean)/(fwhm/2.),2.))\n", " \n", " # Turn the standard Lorentzian into an inverted Lorentzian\n", " lorentz = 1. - lorentz\n", " \n", " return lorentz\n", "\n", "FILL_ME('Write code to fit the Lorentzian to the data')\n", "\n", "lorentz_fit = FILL_ME('Get the fitted Lorentzian curve')\n", "\n", "# Finally, lets get the fitting parameters:\n", "print('Lorentzian amplitude = ', coeff[0])\n", "print('Lorentzian mean = ', coeff[1])\n", "print('Lorentzian FWHM = ', coeff[2])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's compare both fitted profiles to the actual line." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Make a plot!\n", "plt.figure()\n", "plt.plot(wavelength[idx], flux[idx], label='Observed')\n", "plt.plot(wavelength[idx], gauss_fit, label='Gaussian')\n", "plt.plot(wavelength[idx], lorentz_fit, label='Lorentzian')\n", "plt.xlabel('Wavelength ($\\AA$)')\n", "plt.ylabel('Normalized flux')\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(1) Which line profile fits the observed data better? What does this tell us, physically, about how this line is produced?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "FILL_ME('^^ Answer the question above! (You can delete this cell afterward.)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### (c) Measure the equivalent widths of the two strongest lines." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In order to measure chemical abundances, we need to measure the equivalent width ($W_{\\lambda}$ or EW). Let's write a function to estimate the EW by simply summing the area between the spectrum and the continuum." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def measure_EW_area(flux,wavelength,wvl_lower,wvl_upper):\n", " ''' This function measures the equivalent width of a spectral line.\n", " \n", " Inputs:\n", " flux, wavelength -- spectrum to be measured\n", " wvl_lower, wvl_upper -- wavelength limits for the spectral line\n", " \n", " Outputs:\n", " EW -- equivalent width\n", " '''\n", " \n", " # Crop the spectrum to be within the wavelength limits\n", " idx = np.where((wavelength > wvl_lower) & (wavelength < wvl_upper))[0]\n", " \n", " # Define the wavelength difference between two adjacent points in the spectrum\n", " deltawavelength = wavelength[1] - wavelength[0]\n", " \n", " EW = FILL_ME('Compute EW')\n", " \n", " return EW" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now compute the EW of the two strong lines identiifed in (a)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "EW_line1 = measure_EW_area(flux,wavelength,lowerlimit1,upperlimit1)\n", "EW_line2 = measure_EW_area(flux,wavelength,lowerlimit2,upperlimit2)\n", "\n", "print(EW_line1, EW_line2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### (d) Use the curve of growth to determine the column density of the absorbing atoms." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First get the data from the file `growth.txt`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "growth = FILL_ME('Use pandas to get the data.')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot the curve of growth\n", "plt.figure()\n", "plt.plot(growth['log(Nf)'], growth['log(W/lambda)'])\n", "plt.xlabel(r'$\\mathrm{Log}(Nf)$')\n", "plt.ylabel(r'$\\mathrm{Log}(W/\\lambda)$')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have the equivalent widths ($W$) from part (c) and the wavelength ($\\lambda$) from part (a), so we can determine the y-values for our two lines. Let's determine the x-values ($\\mathrm{Log}(Nf)$) for each line.\n", "\n", "We'll use the function `scipy.interpolate.interp1d` to interpolate the data from the file. If you've never used this function before, you may want to read the documentation here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.interpolate import interp1d\n", "\n", "f = FILL_ME('Use interp1d to get a function that interpolates from the y-values to the x-values')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Now determine the values of log(Nf)\n", "logNf_1 = f(np.log10(EW_line1/lambda1))\n", "logNf_2 = f(np.log10(EW_line2/lambda2))\n", "\n", "print(logNf_1, logNf_2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(1) Are these lines in the shallow, partially saturated, or fully saturated region of the curve of growth?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "FILL_ME('^^ Answer the question above! (You can delete this cell afterward.)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the x-values ($\\mathrm{Log}(Nf)$) have a factor of $f$. This factor is known as the \"oscillator strength\" and varies from line to line. It's related to the quantum mechanical transition probability, and can be determined either theoretically or in the lab. \n", "\n", "We want to compute the column density $\\mathrm{Log}(N)$ without the oscillator strength. \n", "\n", "(2) Before we do any calculations: do you predict that the column densities $\\mathrm{Log}(N)$ measured from the two different lines will be the same or different? Why do you think that's the case? (Note that you will be graded on completion rather than accuracy for this question.)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "FILL_ME('^^ Answer the question above! (You can delete this cell afterward.)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(3) What are the values of $f$ for the two lines? (You might want to use the NIST Atomic Spectral Database: https://physics.nist.gov/PhysRefData/ASD/lines_form.html -- you'll need to \"Show advanced settings\" and add $f$ under \"Additional Criteria\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f1 = FILL_ME('Oscillator strength for first line')\n", "f2 = FILL_ME('Oscillator strength for second line')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(4) What are the values of $\\mathrm{Log}(N)$ for the two lines? (Note that the units of $N$ are atoms/cm$^2$)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "logN_1 = FILL_ME('Compute log(N) for first line')\n", "logN_2 = FILL_ME('Compute log(N) for second line')\n", "\n", "print(logN_1, logN_2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(5) Do these answers match your prediction from question (1)? If not, why not? How might our answer change if we estimated $\\mathrm{Log}(N)$ using other lines produced by this atomic species?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "FILL_ME('^^ Answer the question above! (You can delete this cell afterward.)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### (e) Some final questions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(1) The column density $\\mathrm{Log}(N)$ tells us about the number of *absorbing* atoms---that is, atoms that are in the right state to absorb light at a line frequency. Summarize how we could obtain the column density of *all* atoms. (Hint: Your answer should involve the names of two equations. Consider what we learned from Set 4.)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "FILL_ME('^^ Answer the question above! (You can delete this cell afterward.)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hydrogen is the most abundant element in stars, so abundances of other elements are usually written as a ratio relative to helium. Solar abundances are often written in \"epsilon\" notation: $\\log\\epsilon_{\\mathrm{X}}= \\log(N_{\\mathrm{X}}/N_{\\mathrm{H}}) + 12$, where $N_{\\mathrm{X}}$ is the number density of an element X.\n", "\n", "(2) The actual sodium abundance of the Sun is $\\log\\epsilon_{\\mathrm{Na}}= 6.24$. If the column density of hydrogen in the Sun is $N_{\\mathrm{H}}=6.6\\times 10^{23}~\\mathrm{cm}^{-2}$, what is the total column density of *all* sodium atoms? What is the fraction of sodium atoms that are in the right state to produce the strong spectral lines we observed?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "FILL_ME('^^ Answer the question above! (You can delete this cell afterward.)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "When considering stars other than the Sun, we normally write abundances using the following notation: $[\\mathrm{X/H}]=\\log\\frac{(N_{\\mathrm{X}}/N_{\\mathrm{H}})_{*}}{(N_{\\mathrm{X}}/N_{\\mathrm{H}})_{\\odot}}$,\n", "where $N_{\\mathrm{X}}$ is the number density of an element X, and the subscripts $*$ and $\\odot$ refer to some star and the Sun.\n", "\n", "(3) What is the value of [Na/H] for the Sun? What about [Fe/H] for the Sun? What does it mean if a star has [Fe/H] < 0?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "FILL_ME('^^ Answer the question above! (You can delete this cell afterward.)')" ] } ], "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.7.1" } }, "nbformat": 4, "nbformat_minor": 2 }