{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <p style=\"text-align: center;\">PHYS 134L Spring 2024 Lab 3</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-danger\"><b>Due date:</b> Sunday, April 28th, 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": [
    "Read through this entire lab before you start. In this lab, we will look at the brightness of stars as measured in astronomical units called ''magnitudes.'' To complete this lab you should have already read textbook Chapter 3 and Sections 2.3 and 2.4. In this lab, you will be asked to type out some equations. In a jupyter notebook you can do this using standard LaTeX math notation. If you're new to LaTeX you can use [this online equation editor](https://latex.codecogs.com/) to help you along to start. Later in this course you'll be using LaTeX for your final report. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## <p style=\"text-align: center;\">Part 1: Astronomical Magnitudes</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**We'll begin by making a two-dimensional plot of ```cluster1.fits``` and overplotting the locations of the stars given in ```cluster1.cat.```**  Refer to the first lab if you have forgotten how to do this. As in that lab, please use a scale that shows the full dynamic range of the image and that looks nice: no garish color schemes, and few (if any) pixels beyond the limits of your color scale.  The $x$ and $y$ pixel locations might be in different locations than they were in ```object.cat``` from Lab 1, so please check and make sure that you are reading the correct column.  You may open ```cluster1.cat``` in a text editor (like notepad) or on JupyterHub to double-check the contents of each column.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The file ```cluster1.cat``` has many columns besides the $x$ and $y$ positions.  We will now look at the information in the columns ```MAG_ISOCOR``` and ```MAGERR_ISOCOR```.  **First, make a scatter plot showing ```MAG_ISOCOR``` as a function of ```NUMBER```, i.e., simply plot the magnitudes as a function of the order in which they appear in the catalog.  Overplot the error bars given by ```MAGERR_ISOCOR```.**  Check the documentation for ```pyplot.errorbar``` if you are unsure of how to get started."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's hard to interpret that previous plot; it looks like a bit of a mess.  To make it easier to see what is going on, let's sort the stars by magnitude. You'll want to sort the magnitudes and their uncertainties in the same way.  In other words, you want the same uncertainty to remain associated with each measured magnitude. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume that you have read in the catalog file and set an array called ```mags``` equal to the column ```MAG_ISOCOR``` and an array called ```mag_errs``` equal to the column ```MAGERR_ISOCOR```.  Explain why the following code will produce nonsense if you use the array ```mag_errs_sorted``` for the error bars of ```mags_sorted```:\n",
    "```\n",
    "mags_sorted = np.sort(mags)\n",
    "mag_errs_sorted = np.sort(mag_errs)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To sort the arrays in the same way, try ```np.argsort```.  Experiment yourself a bit to get it working, and read the documentation page.  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Now remake your plot of magnitudes with error bars, but sorted by magnitude on the horizontal axis.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do the estimated errors varay with magnitude in a systematic way? Explain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far we have used catalogs output by Source Extractor to obtain only the most minimal information about the objects in the detector field of view--the positions and fluxes from each detected object. We will now look at Source Extractor output files for which we have asked the program to do a more thorough job. Typically we want more kinds of information about each object for one of three reasons: to represent the data in a more convenient way, to display the result of a different way to estimate a number (such as the flux) for which we already have a value, or to give information about the reliability of the results displayed in other columns.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will continue to work with the catalog given by ```cluster1.cat```, but now will load the columns corresponding to ```FLUX_ISOCOR``` and ```FLUXERR_ISOCOR```.  Open the catalog file in a text editor like notepad if you cannot find them.  As usual, remember that Source Extractor counts from 1, while python indexes from 0.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will check the relationship between flux and magnitude, if you need a reminder, check out Section 3.2 of the Burns.  **Make a scatter plot with ```MAG_ISOCOR``` on the $x$-axis and ```FLUX_ISOCOR``` on the $y$-axis.  Once you have done this, adjust the vertical axis to have a logarithmic scaling (check ```pyplot.yscale```).** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you have this scatter plot, try fitting a line by eye.  The mechanics of how to actually do this properly are beyond the scope of this class.  For now, please define a variable $x$ by\n",
    "```\n",
    "x = np.linspace(np.amin(mags,np.amax(mags)),10)\n",
    "```\n",
    "and a variable $y$ by\n",
    "```\n",
    "y = a*x + b\n",
    "```\n",
    "where you determine the best-fit ceofficients ```a``` and ```b``` by eye. Depending on how you have implemented your logarithmic vertical scale you may need to plot either \n",
    "```\n",
    "y = np.exp(a*x+b)\n",
    "```\n",
    "or \n",
    "```\n",
    "y = 10**(a*x+b)\n",
    "```\n",
    "for it to appear as a line on your plot. **Please overplot your best-fit line on the scatter plot of magnitudes against logarithmic flux.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#You code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Write down the definition of astronomical magnitudes in the space below.**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**How closely does your fitted constant ```a``` match expectations?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Approximately what magnitude would yield a flux of 1 given your values of ```a``` and ```b```?  Please show your work.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## <p style=\"text-align: center;\">Part 2: The Sizes of (Images of) Stars</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The catalog file has a column entitled ```FWHM_IMAGE``` (meaning Full Width at Half Maximum of the star image). This parameter is a measure of the width of the star images in units of pixels -- smaller numbers mean sharper images.  You can also estimate these values yourself by looking at the image in ds9.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Open the image in ds9, set the “scale” option to “zscale,” and estimate the size of the images\n",
    "of stars (in pixels). How do these sizes compare with the FWHM IMAGE values in the catalog\n",
    "file?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Now try setting the ds9 ``scale'' to ``min/max'' and ``linear'' and estimate the sizes of a few of the star images (necessarily the brightest ones, since those are all you can see this way).  Use the cursor to read the pixel values.  How do the sizes you estimate in this way compare to the FWHM_IMAGE values in the catalog file?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Set the ``Horizontal Graph'' option in the ``View'' tab, and estimate the FWHM of a bright star from this graph. How does this estimate compare with the other two estimates you gave above?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Last lab, you estimated the angular sizes of a few stars; the largest size you should have computed was around 2~mas.  You also computed the pixel scale in these images; you should have gotten a scale of around $0.^{\\prime \\prime}58$ arcsec/pixel. **Use these two numbers to estimate the angular size of a star in units of pixels.  Please show your work for the unit conversion.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The angular size of a star (in units of pixels) divided by the image FWHM gives a dimensionless number.  If this number is larger than 1, the star is said to be resolved: a larger star makes a larger image.  If this number is much less than 1, the star is said to be unresolved, because the size of the star on the image is set almost entirely by something other than its physical size.  **What is this dimensionless number, angular size/FWHM, for the current image?  Roughly how small would the FWHM have to be for any of the stars to be resolved at all?**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your anwer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**To see if all stars have the same size on the image, make a plot of FWHM ($y$-axis) against stellar magnitude ($x$-axis). As usual, please label your axes and make the plot look nice.  Do you see evidence for a dependence of FWHM on magnitude?**"
   ]
  },
  {
   "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": [
    "For telescopes on the ground, the FWHM in an image is set either  by:  a) the atmospheric turbulence, where the atmosphere blurs out light passing through it, or b) the telescope primary mirror diameter, which introduces diffraction and limits the FWHM of the point spread function to be ~$\\lambda/D$ (in radians). When limited by the atmosphere, a better site gives a shaper image; the size of the atmospheric FWHM we call the ''seeing'' (typically expressed in arcseconds). Try to find on the internet the typical seeing for a few different astronomical telescope sites: Mauna Kea, Paranal and one site of your choice. Compare this to the diffraction limit of a couple telescopes at each those sites, when observing at 600nm. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## <p style=\"text-align: center;\">Part 3: Backgrounds and Photon Noise</p>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ```BACKGROUND``` column in ```cluster1.cat``` contains estimates of the sky background flux (per pixel) in the neighborhood of each identified object. **Make a scatter plot of the sky background as a function of ```X_IMAGE```, using any plot symbols you like, but without lines connecting the points.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What happens if you do use lines to connect the data points? **Explain what the program is doing to generate this unusable plot.**"
   ]
  },
  {
   "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": [
    "**What is a typical value for the sky background? What is the approximate scatter in the reported values?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Within Python you can calculate simple statistics (min, max, mean, standard deviation) for data sets using ```numpy``` functions. Ordinarily, your visual impression of the typical value of a bunch of plotted data points will closely approximate the mean value, and your visual estimate of the peak-to-peak scatter (ignoring rare outliers) will be roughly 4 times the standard deviation. How do mean and standard deviation computed with Python compare with the visual estimates you just made?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "**Now, make a scatter plot ```MAGERR_ISOCOR``` as a function of ```MAG_ISOCOR```.** ```MAGERR_ISOCOR``` is an estimate of the error in ```MAG_ISOCOR```, based on square-root-rule counting error for the number of detected photons in the object plus the underlying sky background (derived from the Poisson distribution). The total number of sky+object photons is\n",
    "$$\n",
    "N_{\\rm phot} = {\\tt GAIN} \\cdot ({\\tt FLUX\\_ISOCOR} + {\\tt BACKGROUND} \\cdot N_{\\rm pix})\n",
    "$$\n",
    "Where ````GAIN```` is the number of detected photo-electrons per signal count (2.36 $e^-$/count for this image, determined from the ```EGAIN```\n",
    "keyword in the FITS header). $N_{\\rm pix}$ is the number of pixels occupied by the object image. With ```ISOCOR``` photometry this number is hard to know for sure, but a reasonable guess is $\\pi ({\\rm FWHM}/2)^2$. **Use this guess, the ```GAIN``` value just given, and your estimate of ```BACKGROUND``` to write an expression for $N_{\\rm phot}$ as a function of ```FLUX_ISOCOR```.**"
   ]
  },
  {
   "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": [
    "**Now use the flux-vs-magnitude expression from Part 1 to write an expression for $N_{\\rm phot}$ as a function of ```MAG_ISOCOR```.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The expected counting error (measured in photo-electrons) is the square root of $N_{\\rm phot}$. To get the noise in units of counts, we must\n",
    "divide by the ```GAIN```. **Write an expression for this counting error as a function of MAG_ISOCOR.**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the tricky part. **What error (in magnitudes) is implied by this error in the measured value of $N_{\\rm phot}$?** Assume that the error is small\n",
    "compared to $N_{\\rm phot}$ itself (i.e., that $N_{\\rm phot} \\ll 1$), and show that if\n",
    "$$\n",
    "{\\rm mag} = m_1 - 2.5 \\log_{10}({\\rm flux})\n",
    "$$\n",
    "(where $m_1$ is the magnitude that yields ${\\rm flux} = 1$), then\n",
    "$$\n",
    "\\delta {\\rm mag} \\approx 1.086 \\frac{\\delta {\\rm flux}}{\\rm flux}.\n",
    "$$\n",
    "Derive this expression below.  You may use standard propagation of uncertainties and take the square root."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your equation here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Now combine this expression and the one for counting error to write an expression for the expected error in magnitudes as a function of the object magnitude.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your equation here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Finally, define integer magnitudes covering the range of magnitudes in ```cluster1.cat```:**\n",
    "\n",
    "```x = np.arange(int(np.amin(mags)), int(np.amax(mags)) + 1)```\n",
    "\n",
    "**and evaluate your expression for magnitude uncertainties at these integer ```x``` values.  Overplot this on your plot of ```MAGERR_ISOCOR``` vs.~```MAG_ISOCOR``` using a thick, red line. Comment on whether they match or not**"
   ]
  },
  {
   "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 4: Consistency and Systematic Errors</p>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we know how to use the data from Source Extractor to estimate the fluxes and magnitudes of stars, and also the precision that simple physics says we should be achieving. But how well do we know our errors, really? Might effects other than photon counting statistics be dominant? And what about systematic errors, which give us consistent and repeatable wrong answers? How can we test for these?\n",
    "\n",
    "\n",
    "Source Extractor can estimate source fluxes in several ways – two of these are isophotal photometry (which deals well with objects having funny shapes) and aperture photometry (which works well for perfectly round objects, as star images are supposed to be). Check section 7.4 of *Source Extractor for Dummies* for an explanation of what these things mean. The picture on p.~41 is particularly helpful.\n",
    "\n",
    "\n",
    "Since there are different ways to estimate what ought to be the same quantities, let’s see if we get the same answers with aperture photometry as we did with isophotal photometry. If not, the differences tell us something about our errors. \n",
    "\n",
    "Make a figure with two graphs, one on top of the other. For this task, you'll want to use subplots sharing the $x$-axis (read the documentation page for ```pyplot.subplots```).  On the top axes, plot ```MAG_ISOCOR``` on the $x$-axis and ```MAG_APER``` on the $y$-axis. You should expect a tight but not perfect correlation between the two measurements, with a few dramatic outliers. The wide range in plotted magnitudes makes it hard to see the errors. **On the lower graph plot ```MAG_ISOCOR``` on the $x$-axis and, on the $y$-axis, plot (```MAG_ISOCOR – MAG_APER```). Choose your y-axis plot range carefully, to show the most information.\n",
    "If necessary, sacrifice a few outliers to show the typical scatter better.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the peak-to-peak scatter in the difference (```MAG_ISOCOR - MAG_APER```), if you ignore extreme outliers? How does this compare with the tabulated uncertainties ```MAGERR_ISOCOR``` and ```MAGERR_APER```? What do you conclude from this?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Identify the stars corresponding to 3 of the most extreme outliers in your upper plot, and look at them using ds9. Can you describe the potential source(s) of these systematic errors?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your answer here*"
   ]
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
