PHY134L/Lab_5/Lab_5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <p style=\"text-align: center;\">PHYS 134L Spring 2024 Lab 5</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\"><b>Due date:</b> Sunday, May 12th, 2024 by 11:59pm, submitted through Gradescope.</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Names: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Enter your name and your partner's name here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this lab we will start exploring how photometry can teach us about the fundamental properties of stars. To do this we'll start by considering their spectra, and what we can learn about the spectra from two images in different filters. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## <p style=\"text-align: center;\">Part 1: Stellar Spectra</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first part of this lab will be an introduction to spectra.  Astronomical data typically comes in two basic forms: images and spectra.  You have already seen images.  Spectra are a bit different: we lose at least some of our ability to resolve an object and see its shape, but we spread out its light into its wavelengths.  The intro slides to this lab included a spectrum as it appears on a detector.  It is a colored spectrum of the Sun, showing many black parts where there is less light.  These are absorption lines and give a whole lot of information about the physics of the Sun. \n",
    "\n",
    "Just like Source Extractor takes the images and extracts the positions and brightnesses of stars, we can extract the flux of a star as a function of wavelength.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this part of the lab, you will use three spectra contained in the files ```spectrum1.dat```, ```spectrum2.dat```, and ```spectrum3.dat```.  The spectra are of three stars.  The data files have two columns: the first is wavelength in Angstroms ($1~{\\rm \\AA}=10^{-10}$~m), and the second is flux density in units of erg/s/cm<sup>2</sup>/$\\AA$.  \n",
    "You may read the files using, e.g.,\n",
    "```\n",
    "spectrum = np.loadtxt(specfile)\n",
    "lam = spectrum[:, 0]\n",
    "flam = spectrum[:, 1]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please be aware that ```lambda``` has a special meaning within python, so don't try to define it as a variable!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a plot of each star's spectrum as a function of wavelength.  Please label the figures nicely using appropriate labels. For example: \n",
    "```\n",
    "plt.xlabel(r'Wavelength ( $ \\AA $ )'\n",
    "```\n",
    "\n",
    "The ```r``` before the string tells matplotlib that it should use Latex to render this.  Please get the units right for the vertical axis.  Hint: if you are unsure of how to write a symbol in Latex, try [http://detexify.kirelabs.org](http://detexify.kirelabs.org). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's try to learn something about these stars.  We can try to approximate the stars as blackbodies using the Planck function (look it up in Wikipedia).  **Write down the Planck function in a markdown cell in your notebook so that it has as close as possible to the same units as the flux density of your spectra (i.e.~make sure it's per Angstrom, not per Hertz or per nm).** The units won't be completely identical: the Planck function will still be per steradian.  This will be important soon."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**By eye, find the temperature and normalization of the Planck function that best fits each stellar spectrum.**  You should choose the temperature that gives a good fit to the peak wavelength (a blackbody isn't a great fit to ```spectrum2.dat``` for reasons that we will discuss later), and multiply the Planck function by a number to put it on the same scale as your spectrum.  This number will be very small.  Make a plot with the three stellar spectra again, with your best-fit Planck functions overplotted.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have fit two parameters by eye to each star's spectrum: the temperature and a normalization factor to the Planck function.  What units does your normalization constant need to have in order for the units to work out? What is the physical meaning of this value?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## <p style=\"text-align: center;\">Part 2: Stellar Colors</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many factors influence the appearance of a star’s emitted energy spectrum, but the most important is the stellar temperature. You already saw this in Part 1: the three stars have almost the same chemical composition, but their spectra look very different.  Their surface temperatures range from around 4000K to around 10000K: Spectrum 1 is from a $\\sim$ 5000K star, Spectrum 2 is from a $\\sim$ 9500K star, and Spectrum 3 is from a $\\sim$ 6200K star."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cluster images that you have been looking at, however, are not spectra.  The telescope takes the light from the stars, focuses it, and uses a filter to let only some colors through.  Only photons within a specified wavelength range actually make it onto the detector.  Here we will investigate how the visible-light spectrum changes with temperature, and how these changes are manifested in the stars’ photometric colors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Lab 5 directory contains three filter files and three additional stellar spectra.  The stars are at three different temperatures:\n",
    "\n",
    "```f05.dat``` = 9550K\n",
    "\n",
    "```f10.dat``` = 6110K\n",
    "\n",
    "```f15.dat``` = 4350K\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These new files contain the energy flux $F_\\lambda$ (proportional to ${\\rm erg}\\,{\\rm s}^{-1}\\,{\\rm cm}^{-2}\\,{\\rm \\AA}^{-1}$), where wavelength $\\lambda$ is measured in Angstroms (${\\rm \\AA}$), and $1~{\\rm Angstrom} = 0.1~{\\rm nm}$. The flux values have been arbitrarily normalized so they have a value of 100 units near 5500\\,${\\rm \\AA}$.  Remember from Part 1 that the normalization of the stellar spectrum told you the angular diameter of the star on the sky.  The shape of the spectrum tells you about the physical conditions of the star, mainly its temperature and chemical composition. \n",
    "\n",
    "Take the spectra of the three stars ```f05.dat```, ```f10.dat```, and ```f15.dat```, and plot them together on the same graph, using different colors.  Also plot the Planck functions at these three temperatures on the same axes.  You will notice that the hottest star is a particularly poor match to the blackbody; this is due to the hydrogen in the star's atmosphere."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using a CCD detector with colored filters, we measure the energy flux from stars integrated over the bandpass of each filter:\n",
    "\n",
    "\\begin{equation*}\n",
    "    F = \\int_0^\\infty T[\\lambda] F_\\lambda d\\lambda,\n",
    "\\end{equation*}\n",
    "\n",
    "where $T[\\lambda]$ is called the *transfer function* and describes how efficiently we record energy from photons of a given wavelength.  $T[\\lambda] = 1$ means that our instrument records every photon at a given wavelength, while $T[\\lambda] = 0$ means that it never records a single photon.  **We will do the integrals for the filters $B$ and $V$.**  It's worth spending a bit of time on how to do an integral numerically, as this is something many of you will have to do again and again."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The integral is defined as the limit of the Riemann sum,\n",
    "\n",
    "\\begin{equation*}\n",
    "    \\int_a^b f(x) dx = \\lim_{N \\rightarrow \\infty} \\sum_{i = 0}^{N - 1} f[x_i] \\Delta x_i,\n",
    "\\end{equation*}\n",
    "where $\\Delta x_i$ is the spacing between values of $x_i$, $x_0 = a$, and $x_N = b$.  Let's see how we can write this in python.  First, we'll define $x$ to run from $a$ to $b$ with $N + 1$ points in total.  For now, define $a = 0$, $b = 1$, and $N = 1000$ (we'll change these later when we actually integrate over the band passes).\n",
    "```\n",
    "x = np.linspace(a, b, N + 1)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now check: use the print function to verify that ```x[0] = a``` and ```x[N] = b```.  In your case, the spacing is uniform, so $\\Delta x_i$ is the same for all points $i$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use your value for $\\Delta x$, the ```numpy``` function ```sum```, and your values of ```x``` to numerically evaluate the integral\n",
    "\n",
    "\\begin{equation*}\n",
    "    \\int_0^1 \\sin x dx\n",
    "\\end{equation*}\n",
    "\n",
    "and compare to the actual value of the integral.  Remember that your (left) Riemann sum should stop at $N - 1$, which you can write using index notation as ```[:N]```.  Show that your answer is closer to the truth if $N = 2000$, but is farther from the truth if $N = 500$.\n",
    "\n",
    "Now we can do an integral numerically, but how do we integrate \n",
    "\n",
    "\\begin{equation*}\n",
    "    F = \\int_0^\\infty T[\\lambda] F_\\lambda d\\lambda? \n",
    "\\end{equation*}\n",
    "\n",
    "We need to define a wavelength array, and then be able to call $T[\\lambda]$ and $F[\\lambda]$ just as we did above with sine.  We can do this by first reading in the data files and then *interpolating* between the data points.  We will use the package ```scipy.interpolate```, and we will use cubic polynomials (called *splines*) to interpolate:\n",
    "\n",
    "```from scipy import interpolate```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this point, ```python``` actually makes our lives pretty easy.  First, let's read in one of the spectral files (first column is wavelength, second column is flux density $F_\\lambda$):\n",
    "```\n",
    "spec_data = np.loadtxt('f05.dat')\n",
    "```\n",
    "and then\n",
    "```\n",
    "Flam = interpolate.interp1d(spec_data[:, 0], spec_data[:, 1], kind='cubic')\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now ```Flam``` is a function that you can call just like ```np.sin```.  Calling this function will return an error if you ask it for the value outside the range of $x$-values you supplied to ```interpolate.interp1d```. **For the transfer function, we may assume that it is zero outside of the wavelength range given.  Read the documentation page (check the keywords ```bounds_error``` and ```fill_value```) to make the function behave in this way.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to compute magnitudes.  **Use the tools you have developed above to integrate the three spectra over the transfer functions given by ```Bfilt.dat``` and ```Vfilt.dat```.**  The definition of magnitude includes a zero point:\n",
    "\n",
    "\\begin{equation*}\n",
    "    M = -2.5 \\log_{10} F - Z.\n",
    "\\end{equation*}\n",
    "\n",
    "For now, ignore the zero points $Z$ (let $Z=0$).  What are the magnitudes of the three stars in the two bands?  Please fill in the following table.  Also, please compute the difference between the $B$ and $V$ band magnitudes (the $B-V$ color).  This is the ratio between the flux in the $B$ band and the flux in the $V$ band."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---| ```f05.dat``` | ```f10.dat``` | ```f15.dat1```\n",
    "---|---|---|---\n",
    "B |---|---|---\n",
    "V |---|---|---\n",
    "B-V |---|---|---\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Does a high value of $B-V$ indicate a star that is very blue (puts out lots of blue light relative to yellow/green light), or a star that is very red?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the qualitative relation between temperature and B − V color?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you can imagine that only using images of a star in two filters, we can learn a lot about its fundamental properties! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So colors provide a way to estimate stellar temperatures. Suppose we know the temperature -- then what else do we know?  Stellar evolution theory provides relations between the characteristics of stars at the time they are born (mass and composition), and their properties (including temperature and luminosity) at any time in their later lives. We'll explore this more in Lab 6. For a quick introduction to (or refresher on)\n",
    "many of the basic notions, see [here](http://www.tim-thompson.com/hr.html). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}