{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <p style=\"text-align: center;\">PHYS 134L Spring 2024 Lab 4</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\"><b>Due date:</b> Sunday, May 5th, 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": [
    "The previous lab explored the systematic differences that result from doing a measurement in two slightly different ways. Next, we will try to replicate the measurements made by Source Extractor using ```Astropy``` packages. The main new package we will use is called ```photutils```, an astropy package used for photometry--the measurement of the brightness of astronomical objects. You can find all the documentation for ```photutils``` [here.](https://photutils.readthedocs.io/en/stable/)\n",
    "\n",
    "We'll again be using ```cluster.fits``` from our previous labs. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## <p style=\"text-align: center;\">Part 1: Background Estimation</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step will be determining the background level, so that we don't overestimate the amount of light coming from any given star. We can do this in several ways. First, open ```cluster.fits``` in DS9. **Hover your mouse around the background and try to estimate by eye what you think the background level is, enter it here:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to estimate the background here is by looking at the median of the image. **Load the fits file into a new jupyter notebook and get the median value of the whole array.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**What kind of probablems can you think that might make a simple median a poor estimate?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the photutils background estimation [documentation](https://photutils.readthedocs.io/en/stable/background.html) and test out several methods of background subtraction. Determine the best method to use on this data. Play with image scaling and colorbars to help you decide the best method. Show your work below. Describe below your final decision on the best method here and why. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## <p style=\"text-align: center;\">Part 2: Source Identification</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we've settled on a method to extract the background we can move on to automatically detecting the sources. ```photutils``` has two options for this: ```DAOStarFinder``` and ```IRAFStarFinder```. See the [documentation here](https://photutils.readthedocs.io/en/stable/detection.html). We'll stick with ```DAOStarFinder```, but feel free to read up both methods and experiment with them. \n",
    "\n",
    "**Follow the example in the documentation and generate a list of detected sources.** There are several knobs you can easily turn to change your results: the FWHM, threshhold (i.e. how many ```stds``` above the background a source needs to be before being considered real) and how the ```std``` is calculated. Play with these until you are satisfied. Try to maximize the number of detected sources. **Display the cluster image and overplot the locations of the stars that you found to confirm that things are behaving as you expect.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You'll notice some obviously bright stars don't get identified. Go look at the image in DS9. Do the unidentified stars have anything in common? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do these same stars get picked up by the Source Extractor? If so, how reliable do you think their aperture photometry will be?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here (if need be)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Source Extractor may have methods to compensate for these kinds of issues (such as the \"FLAG\" column), but this kind of thing should raise your suspicion and motivate you to read the documentation to find out. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## <p style=\"text-align: center;\">Part 3: Aperture Photometry</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we're moving on to extracting the photometry. Here we'll be carrying out photometry using a technique called \"aperture photometry\", where we simply add up all the counts within a circular aperture. There are other types of photometry, such as \"weighted PSF\" photometry, that may be useful in your final project. Another type of photometry is \"Isophotal photometry\", which we encountered earlier in the course. You can read about that in the Sextractor documentation or elsewhere online. \n",
    "\n",
    "See the ```photutils``` documentation [here](https://photutils.readthedocs.io/en/stable/aperture.html) on aperture photometry. **Read up on the documentation and try it out yourself.**\n",
    "\n",
    "To get started quickly you can ignore some of these sections in the documentation: sky apertures, pixel masking, aperture masks, aperture photometry using Sky Coordinates. \n",
    "\n",
    "**Extract photometry using circular apertures for all the sources you found using the DAOStarFinder. Don't forget to subtract off the background (using whatever method you settled on earlier) before calculating the photometry. What do you think the radius of your circular apertures should be and why?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Match your star list with the Source Extractor star list (```cluster1.cat```) by matching the closest x and y positions from the two catalogs. Plot your newly extracted photometry against the ```FLUX_APER``` value from Source Extractor. If they don't match, try to explain a few reasons you think this might be the case.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your Code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Play with the various input parameters to the photometric extraction funtion until you get the two catalogs to match as well as you can. One key parameter you can change is the aperture radius, which can make a difference especially if the background isn't sufficiently subtracted. What radius did you have to choose to match the two catalogs? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If changing the radius has such an important impact on the final answer, how does one choose an optimal radius? \n",
    "\n",
    "One way to do this is by optimizing the radius to maximize the signal to noise ratio (SNR) in your extracted photometry. To do this we need to estimate the noise in our apertures. The simplest way to do this is to pass to the ```aperture_photometry``` function a map of the expected noise in our image. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To calculate the noise, we'll first assume that the entire image is limited by the photon noise from our stars and background. This might not be strictly true, as the fainter parts of the image might be limited by the read-out noise of the detector. \n",
    "\n",
    "To create the error map, you'll want to manipulate the original data array, before background subtraction. Discuss with your partner why you want to do this before background subtraction. To calculate the error at each pixel, consider the value at that pixel to be the mean of a Poisson distribution, and the error at that pixel will be the square root of the variance of that distribution. Keep in mind that the Poisson distribution applies to the number of *photons* the detector sees, not the number of counts in your image (recall the previous lab's discussion of the detector gain). **Write down the equation for calculating the error in each pixel here:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using that equation to create an error map for the whole data array and pass that into the ```aperture_photometry``` function so that you can get an estimate of the flux error on each star. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pick 3 random stars from your list and make plots of the SNR as a function of aperture radius (i.e. the aperture sum divided by its error). At what radius do these stars maximize their SNRs?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If they are different, why do you think that is the case?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**For each star calculate its optimal aperture sum and error by maximizing the signal to noise. Convert this to a magnitude as: $$m = -2.5log_{10}(N_{aperture\\_sum}),$$ where $N_{aperture\\_sum}$ is the optimal aperture sum. For each magnitude calculate an error.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## <p style=\"text-align: center;\">Part 4: Photometric Calibration</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Leaving aside the issues of error and uncertainty, we still have to put our measured magnitudes on a standard scale, such that the numbers we derive can be directly compared with other work. This comes down to measuring the constant $m_1$:\n",
    "${\\rm mag} = m_1 - 2.5 \\log_{10}({\\rm flux}),$\n",
    "so that it is reasonably consistent with other work. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the [Aladin Lite](https://aladin.cds.unistra.fr/AladinLite/) online tool to determine the magnitudes of the same 4 stars that you used for the image-scale calculation in the previous lab. Find your stars in Aladin, and then use the \"SIMBAD\" catalog (the little box on the right-hand side) to find the magnitude of your targets. You'll have to pick the appropriate filter (either the same one as the data, or not that covers similar wavelengths). Make a table below that shows the star number from the Figure in Lab 3, its magnitude from *Aladin*, and (by matching $\\{x,y\\}$ coordinates) the magnitude that you calculated with ```photutils```. Calculate the average difference between your magnitudes and those from *Aladin*, and the standard deviation of this difference:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Star| Aladin Mag | ```photutils``` Mag | difference \n",
    "---|---|---|---\n",
    "1 |---|---|---\n",
    "2 |---|---|---\n",
    "3 |---|---|---\n",
    "4 |---|---|---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Write an expression for magnitudes in the *Aladin/SIMBAD* catalog (assumed here to be the true on-sky magnitudes) in terms of your derived magnitudes.**  This expression will convert your measured magnitudes (based on image counts) into on-sky magnitudes. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume the brightest of the stars in your list of 4 has an **absolute** magnitude of 2.75.  What is its distance (in parsecs) based on its apparent magnitude? What uncertainty do you assign to this distance estimate, given only the uncertainty that you just calculated in estimating *Aladin* magnitudes from the magnitudes you calculated? Justify your answers."
   ]
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
