{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "plt.rcParams.update({'font.size': 16})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set cluster name and position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_name, cluster_ra, cluster_dec = 'Abell2065', 230.62156, 27.70763\n",
    "#cluster_name, cluster_ra, cluster_dec = 'Abell2063', 230.77116,  8.60859 \n",
    "#cluster_name, cluster_ra, cluster_dec = 'Abell1795', 207.21886, 26.59160"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read CSV file containing SDSS data (assumes the file lives in same directory as your notebook, and is named cluster_name_SDSS.csv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.io import ascii\n",
    "SDSS=ascii.read(cluster_name+'_SDSS.csv')\n",
    "print('{} objects in table'.format(len(SDSS)))\n",
    "print('Column names: ',SDSS.colnames)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot a (logarithmic) histogram showing the number of objects as a function of r-band magnitude. Also plot a line at the SDSS limiting r-band magnitude of r=22.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "bins=np.arange(10,30,0.2)\n",
    "plt.hist(SDSS['r'],bins=bins)\n",
    "plt.yscale('log')\n",
    "plt.axvline(x=22.2,color='red')\n",
    "plt.xlabel('r magnitude')\n",
    "plt.ylabel('N');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot r mag uncertainty versus r mag. If you want your mags accurate to 10% or better, what is the rough magnitude limit of your analysis?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Fig,(FullPlot,ZoomPlot) = plt.subplots(1,2,figsize=(8,4))\n",
    "\n",
    "for subplot in [FullPlot,ZoomPlot]:\n",
    "    subplot.scatter(SDSS['r'],SDSS['modelMagErr_r'],s=0.1) ## Note s=0.1 -- there are many points, so make them small!\n",
    "    subplot.set_xlabel('r mag')\n",
    "    subplot.set_ylabel('r mag err')\n",
    "    subplot.axhline(y=0.1,color='red')\n",
    "    subplot.axvline(x=21.5,color='red')\n",
    "    \n",
    "FullPlot.set_ylim(0,1)\n",
    "ZoomPlot.set_ylim(0,0.2)\n",
    "Fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot dec vs ra (so a sky map) of resolved sources (type=3) brighter than r=20. Make another for unresolved sources (type=6) brighter than r=20. Think about the differences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fig,(ResPlot,UnresPlot) = plt.subplots(1,2,figsize=(8,4))\n",
    "\n",
    "wantmag=SDSS['r']<20\n",
    "\n",
    "want=np.logical_and(SDSS['type']==3,wantmag)\n",
    "ResPlot.scatter(SDSS['ra'][want],SDSS['dec'][want],s=1)\n",
    "\n",
    "want=np.logical_and(SDSS['type']==6,wantmag)\n",
    "UnresPlot.scatter(SDSS['ra'][want],SDSS['dec'][want],s=1)\n",
    "\n",
    "for subplot in [ResPlot,UnresPlot]:\n",
    "    subplot.set_aspect('equal')\n",
    "    subplot.invert_xaxis()\n",
    "    subplot.set_xlabel('RA')\n",
    "    subplot.set_ylabel('Dec')\n",
    "    \n",
    "Fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Work out the angular distance (in arcsec) of each object from the cluster center. This uses astropy's \"SkyCoord\" functionality to calculate angular distances, and makes a new column in the SDSS table called 'angdist_arcsec'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.coordinates import SkyCoord\n",
    "from astropy import units as u\n",
    "cluster_pos=SkyCoord(cluster_ra,cluster_dec,unit='deg',frame='icrs')\n",
    "gal_coo=SkyCoord(SDSS['ra'],SDSS['dec'],unit='deg')\n",
    "SDSS['angdist_arcsec']=gal_coo.separation(cluster_pos).arcsec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make a plot of positions on the sky (ra,dec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6,6))\n",
    "plt.scatter(SDSS['ra'],SDSS['dec'],s=1)\n",
    "\n",
    "# we want to do these next two things any time we make an RA/Dec plot\n",
    "# this sets the plot XY equal in coordinate degrees (but doesnt factor in the cos(dec) distortion)\n",
    "plt.gca().set_aspect('equal')\n",
    "# this inverts the ra axis so that east is to the left (ie usual astronomical orientation)\n",
    "plt.gca().invert_xaxis()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Do the same thing, but this time only for extended objects (ie galaxies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gals = SDSS['type']==3 # Extended objects are things with TYPE = 3.\n",
    "plt.figure(figsize=(6,6))\n",
    "plt.scatter(SDSS['ra'][gals],SDSS['dec'][gals],s=1)\n",
    "\n",
    "plt.gca().set_aspect('equal')\n",
    "plt.gca().invert_xaxis()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Well, that's a mess. How about doing it just for bright galaxies?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's call a galaxy 'bright' if it's brighter than a given r-magnitude\n",
    "r_maglim = 20.0\n",
    "brightgals = np.logical_and(gals,SDSS['r']<r_maglim)\n",
    "\n",
    "plt.figure(figsize=(6,6))\n",
    "plt.scatter(SDSS['ra'][brightgals],SDSS['dec'][brightgals],s=1)\n",
    "\n",
    "plt.gca().set_aspect('equal')\n",
    "plt.gca().invert_xaxis()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Redo that, overlaying with red X's the ones that have redshifts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(6,6))\n",
    "plt.scatter(SDSS['ra'][brightgals],SDSS['dec'][brightgals],s=1)\n",
    "\n",
    "# remember that if an object doesn't have a measured redshift, we had set its redshift to -999\n",
    "brightgals_withz = np.logical_and(brightgals,SDSS['redshift']>0)\n",
    "plt.scatter(SDSS['ra'][brightgals_withz],SDSS['dec'][brightgals_withz],s=20,marker='x',color='red')\n",
    "\n",
    "plt.gca().set_aspect('equal')\n",
    "plt.gca().invert_xaxis()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Let's examine the structural fits, whether a galaxy is better fit by an exponential model or a deVaucouleur model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# just select resolved objects within 30 arcmin of the cluster center\n",
    "resolved_inside30 = np.logical_and(SDSS['type']==3,SDSS['angdist_arcsec']<1800)\n",
    "\n",
    "# lnLExp_g is the log-likelihood that the object is well fit by an exponential disk model\n",
    "# lnLDeV_g is the log-likelihood that the object is well fit by a deVaucouleur elliptical galaxy model\n",
    "#  so if lnLExp_g > lnLDeV_g, the object is more likely to be a disk galaxy\n",
    "# and if lnLExp_g < lnLDeV_g, the object is more likely to be an elliptical\n",
    "# but remember, this doesn't mean it *is* a disk or an elliptical!\n",
    "\n",
    "plt.scatter(SDSS['r'][resolved_inside30],SDSS['lnLExp_g'][resolved_inside30]-SDSS['lnLDeV_g'][resolved_inside30],s=1)\n",
    "plt.xlim(12,25)\n",
    "plt.ylim(-4000,4000)\n",
    "\n",
    "plt.axhline(y=0,color='red')\n",
    "plt.axvline(x=19,color='black',ls=':')\n",
    "plt.xlabel('r mag')\n",
    "plt.ylabel('lnLExp_g - lnLDeV_g')\n",
    "plt.text(20,3500,'Disk-like')\n",
    "plt.text(20,-3500,'Elliptical-like')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Look at an interactive table of galaxies (type=3) projected inside 10 arcminutes (10*60 = 600 arcsec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "want = np.logical_and(SDSS['type']==3,SDSS['angdist_arcsec']< 600)\n",
    "# uncomment one of these two:\n",
    "#SDSS[want].show_in_notebook() # to show it in the jupyter notebook\n",
    "SDSS[want].show_in_browser(jsviewer=True) # to show it in a browser page\n",
    "\n",
    "# Note, I do not recommend doing this for the full dataset, it's too big!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot redshift versus r-band magnitude for galaxies that have a measured redshift and are within 30 arcminutes of the cluster center. Using the plot as a guide, work out a quantitative, statistical estimate of the cluster redshift. If you wanted to define a \"spectroscopically confirmed cluster member\", how might you do it?\n",
    "\n",
    "Tip: Look at the redshift plot, decide what range of redshifts define the cluster. It'll help to set your ylimits to show a redshift range of z=0-0.1, where the cluster should be. \n",
    "\n",
    "Then make a selection on galaxies with redshifts in that range, and calculate the average redshift of just those objects.\n",
    "\n",
    "Then make a new flag to indicate \"spectroscopically confirmed galaxies\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code to plot and work out redshift info goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now that we know the redshift of the cluster, let's work out how the angular size of a 1 Mpc radius, and plot things inside/outside that radius."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_redshift = 0.05 # REMEMBER, YOU NEED TO SET THIS TO THE PROPER VALUE FOR EACH CLUSTER\n",
    "\n",
    "# import the astropy cosmology package, setting the cosmology to the Planck18 result\n",
    "from astropy.cosmology import Planck18 as cosmo\n",
    "\n",
    "# calculate the angular size distance to the cluster\n",
    "DA = cosmo.angular_diameter_distance(cluster_redshift) # this is the angular diameter distance in Mpc\n",
    "\n",
    "# and now we can work out the angular size that corresponds to a cluster radius of 1 Mpc\n",
    "# Here we are just ysing the small angle approximation: \n",
    "#    angle[radians] = physical size / distance\n",
    "# (and remember physical size and distance must be in the same units!)\n",
    "# Then we convert from radians to arcseconds by multiplying by 206265\n",
    "One_Mpc_In_Arcsec = (1.0/DA.value) * 206265\n",
    "print('At z={:.4f}, 1 Mpc projects to {:.1f} arcsec'.format(cluster_redshift,One_Mpc_In_Arcsec))\n",
    "\n",
    "# now make a flag that indicates if a galaxy is inside 1 Mpc or not:\n",
    "inside1Mpc = SDSS['angdist_arcsec']<One_Mpc_In_Arcsec\n",
    "\n",
    "# while we are playing with astropy.cosmology, let's get the luminosity distance to the\n",
    "# cluster and work out the distance modulus to convert apparent magnitude to absolute\n",
    "# magnitude, in case we need it later:\n",
    "DL = cosmo.luminosity_distance(cluster_redshift) # this is the luminosity distance in Mpc\n",
    "distmod = 5.*np.log10(DL.value*1e6) - 5\n",
    "# so if we ever need an r-band absolute magnitude we can calculate it as r_absmag = SDSS['r'] - distmod\n",
    "\n",
    "# for more details on cosmological distance definitions, go to my ASTR222 website:\n",
    "#    http://burro.case.edu/Academics/Astr222/Cosmo/Observing/observing.shtml\n",
    "\n",
    "\n",
    "#OK, let's plot things\n",
    "\n",
    "plt.figure(figsize=(6,6))\n",
    "\n",
    "# first plot all bright galaxies\n",
    "plt.scatter(SDSS['ra'][brightgals],SDSS['dec'][brightgals],s=1)\n",
    "\n",
    "# now select all galaxies which are bright AND inside 1 Mpc\n",
    "want = np.logical_and(brightgals,inside1Mpc)\n",
    "plt.scatter(SDSS['ra'][want],SDSS['dec'][want],s=1,color='red')\n",
    "\n",
    "plt.gca().set_aspect('equal')\n",
    "plt.gca().invert_xaxis()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot g-r color versus r magnitude for all resolved sources projected within 1 Mpc of the cluster center -- this is a color-magnitude diagram (CMD) for galaxies.\n",
    "\n",
    "*Tip: Remember to use appropriate ranges: a reasonable r-mag range is r = 13 to 24 and a reasonable g-r color range is g-r = -1 to +2.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the CMD again, then overplot in a different color the CMD for spectroscopically confirmed cluster members. Again, restrict it to resolved sources projected within 1 Mpc of the center.\n",
    "\n",
    "*Tip: Do a dual selection: galaxies within the 1 Mpc projected radius AND within the range of redshifts that you defined for the cluster. Then plot the CMD for those galaxies on top of the one you made in the previous step for all objects projected within 1 Mpc. Use a different color/symbol for the spectroscopically confirmed galaxies*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Identify the bluest spectrosocopically confirmed galaxies, find their ra and dec, and then find them using Skyserver's \"Navigate\" function. What do they look like morphologically? Spectroscopically?\n",
    "\n",
    "*Tip: Define a new column of data for the SDSS data table to hold the g-r color: SDSS['g-r']=SDSS['g']-SDSS['r']. The do a show_in_browser call on objects within the cluster redshift range. Sort on the g-r column to find the bluest objects and look at their coordinates. Then find them using Navigate.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Identify the most luminous spectroscopically confirmed galaxies and look at them in Navigator. What do they look like morphologically? Spectroscopically?\n",
    "\n",
    "*Tip: Now sort your show_in_browser table on the r-magnitude and look at the coordinates of the brightest object. Then find it using Navigate.*\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Identify the highest redshift objects in the field (they won't be in the cluster, obviously!) and find them in Navigator. What do they look like, morphologically and spectroscopically?\n",
    "\n",
    "*Tip: Don't worry if SDSS classifies it as a galaxy or not, or if it is inside 1 Mpc or not. Just do a show_in_browser call on all objects with a redshift (whether or not they are in the cluster redshift range), and sort on the redshift to find the coordinates of the highest redshift object.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now make a CMD color-coded by disk vs elliptical likelihood, just for galaxies projected inside 1 Mpc of the cluster center.\n",
    "\n",
    "*Tip: define a flag for disky things by saying **disky = SDSS['lnLExp_g'] > SDSS['lnLDeV_g']** and vice versa for E's. Then plot the CMD for disky things in blue and elliptical things in red.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}