plotcmd.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Make pretty \"images\" showing the color magnitude diagrams for SDSS galaxies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The SDSS DR 7 photometry (> 100,000 randomly seelcted galaxies within the full >1 million galaxy sample) will be used to generate galaxy CMD. We don't do careful filtering here (i.e., there may still be AGN in the sample) but we use the kcorrected (i.e., galaxy SED) magnitudes where available, aiming to select mostly galaxies well-fit by galaxy templates. The goal here is to show where typical galaxies (z>0.01) reside in color-magnitude space. We will overlay our sample of Swift pair galaxies, drawn from the Barton paper of close pairs of galaxies, on the CMD as a comparison.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### initialize with import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt \n",
    "import sys, os\n",
    "import astropy.io.fits as fits\n",
    "import numpy as np\n",
    "from astropy.cosmology import FlatLambdaCDM\n",
    "cosmo = FlatLambdaCDM(H0=70, Om0=0.3, Tcmb0=2.725)\n",
    "\n",
    "# ignore runtime warnings\n",
    "\n",
    "if not sys.warnoptions:\n",
    "    import warnings\n",
    "    warnings.simplefilter(\"ignore\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdss = fits.open('data/sdssphot_random.fits')\n",
    "data=sdss[1].data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialize sample: \n",
    "Get colors and absolute magnitudes from this SDSS-DR7 sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get required data\n",
    "z=data['z']\n",
    "d_l_pc=cosmo.luminosity_distance(z).value*1.e6\n",
    "\n",
    "# get apparent magnitudes\n",
    "mag=data['kcor_mag']\n",
    "gmag=mag[:,0]\n",
    "rmag=mag[:,1]\n",
    "imag=mag[:,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 91871 galaxies included in our final analysis... \n"
     ]
    }
   ],
   "source": [
    "# determine where the values are appropriate. \n",
    "#good = np.where((gmag > 0) & (rmag > 0) & (imag > 0))\n",
    "good=np.where((z>0.01) & ( gmag > 0 ) & ( rmag > 0 ) & ( imag > 0 ) & ( d_l_pc > 0))\n",
    "print(\"There are %i galaxies included in our final analysis... \" % np.size(good))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for good galaxies, calculate remaining useful data:\n",
    "dm=5*(np.log10(d_l_pc[good])-1)  # distance modulus\n",
    "\n",
    "# absolute magnitudes\n",
    "g_absmag = gmag[good] - dm\n",
    "r_absmag = rmag[good] - dm\n",
    "i_absmag = imag[good] - dm\n",
    "\n",
    "col_g_r = g_absmag - r_absmag\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make Plot:\n",
    "Now that we have all the data for the plot, make 2D histogram showing color-magnitude diagram for the ramdom SDSS galaxies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAHwCAYAAAC7apkrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xt4XHd17vF3WVfrYjuxLBHbccIjHE6SKiWtmkJamrQkbUJPE0qhTWhLOIWG9pDynNJySeGES08p0Au0JZSmlMP1kFLKJSWhLuYSWpzQKJDGOCaplYsjbDS+yrJkXSyt88eMyw9Fkr2Xpb1Hmu/nefxEM7PXrDV7z4zebO2Zbe4uAAAAAGUrih4AAAAAqCYEZAAAACBBQAYAAAASBGQAAAAgQUAGAAAAEgRkAAAAIEFABrBsmdnlZjZQ9BwLycy+YGY3LOD9fdXMXrFQ9zfL/e8ws8sX6/4BYDEQkAFUPTN7iZn1mdlRM9tbCYk/WcAcXzUzN7MfnnH9ZyvXX77YM7j71e7+4Urfl5nZvy12z9mY2bmVx3y08m/QzD5vZlfOmPdCd/9qETMCQBQBGUBVM7PXSHqPpLdL6pK0SdL7JF27yH3r5rjpEUkvTZZbK+nZkvYt5jxVbI27t0n6YUlflPQZM3vZYjc1s/rF7gGgdhGQAVQtM1st6W2SXuXun3b3EXefdPd/cvfXVpZpMrP3mNmeyr/3mFnTHPd3fmUv8OHKn/6vSW77kJn9tZndZWYjkn56jrE+LulXkgB9vaTPSJpI7usSM7un0mevmb3XzBqT23/WzB42syEze5+Z3X3iMIcTe4XN7E/N7JCZPWZmVye1XzWzV5jZ+ZLeL+k5lT24h9Pbk+V/YC+zmV1pZt+p9H6vJJuxjn7DzHZWem8xs3Pm3EAJd/+eu/+FpLdIeqeZrajc3+NmdsUCrZevm9m7zeygpLeYWbeZfdnMDpjZfjP7uJmtSe7vcTN7rZk9aGYjZvZ3ZtZV+QvEsJltNbMzTuXxAagtBGQA1ew5kppVDqBzeaPKe3CfpfJezEskvWnmQmbWIOmfJP2LpE5JvyPp42b2zGSxl0j6I0ntkuY6dGGPpIck/Wzl8kslfWTGMlOSfldSR+UxPE/S/6zM0SHpU5JulrRW0sOSLp1R/+OV6zskvUvS35nZDwRZd98p6bck3ePube6+RidR6f2PKq+fDkn9kn4iuf0Fkv5A0gslrZP0r5I+cbL7neHTKq/fZ85y20Ksl0cr9/9HKof7P5a0XtL5ks5WOaCnfknSlZLOk/QLkr5QeYwdKv8OfHXGxwegBhCQAVSztZL2u/vxeZb5VUlvc/eSu++T9FZJvz7Lcs+W1CbpHe4+4e5flvR5lfcAn/A5d/+6u0+7+9g8PT8i6aWVcL3G3e9Jb3T3+939Xnc/7u6PS/obSZdVbn6+pB2VPeLHJf2lpO/NuP8n3P1v3X1K0oclnaXy4SWn6/mSHnL3T7n7pMqHrqS9Xynpj919Z2W2t0t61qnuRa7YU/nvmTNvWID1ssfd/6pSf8zdd7n7F919vLLt/zy5vxP+yt0H3f27Kgf+b7j7t9x9XOX/8bo4w2MDUCM4hgtANTsgqcPM6ucJyeslPZFcfqJy3WzLPenu0zOW3ZBcfvIU5/q0pD+rzPfRmTea2Xkqh7VeSS0qv9fen85xYll3d3vqN218L7l9tLLzuO0UZ5vPbL3Tx3yOpL8wsz9LrjOV11G6judzYn0enHnDAqyXH9g+ZtapcpB+rsp7/VdIOjSjZjD5+dgslxdivQJYZtiDDKCa3SNpTNIL5llmj8rB7oRN+v5ezJnLnX3i2Nhk2e8ml/1UhnL3UZX/VP/bmiUgS/prSd+RtNndV6n8J/0Th0jslbTxxIKVQyc2PuUeTs1s846oHD5PeFry816VD0NIe5+d3P6kpFe6+5rk30p335Zhpl+UVFL5EImZTne9zHy8f1y57qLK/f2aZhxTDQARBGQAVcvdhyTdIulWM3uBmbWYWYOZXW1m76os9glJbzKzdZXjWG+R9LFZ7u4bKofH11Xu43KVj0m9PTjeH0i6rHKowEztko5IOmpm/03lIH3CnZJ6Ko+nXtKr9IMhNotBSRvTD7pJekDSCyvr6hmSXj6j94Vm9sJK71fP6P1+STeb2YVS+UOSZvbiUxmk8uG3myS9WdLNM/bUn7DQ66Vd0lFJh81sg6TXnsqsAHAyBGQAVc3d/1zSa1T+YNk+lfdy3iTps5VF/o+kPkkPStou6ZuV62bez4SkayRdLWm/yl8V91J3/05wrj3uPtcH+X5f5Q/8DUv6W0l/n9Ttl/RilT98d0DSBZX5xwNjfFnSDknfM7P9leverfI3agyqfPzyx2fp/Y5K782Svp7c/hlJ75R0u5kdkfRtldfXfA5b+Vs/tqt8HPGL3f2Dcyy70OvlrZJ+RNKQygH70yeZFQBOibmf0l8UAQCLoHLIx4CkX3X3rxQ9T7VgvQAoEnuQASBnZvZzZrbGyt/XfOI43HsLHqtwrBcA1YKADAD5e47K30G8X+XjoF/g7seKHakqsF4AVAUOsQAAAAAS7EEGAAAAEsvuRCEdHR1+7rnnFj0GAAAACnD//ffvd/d1p3Mfyy4gn3vuuerr6yt6DAAAABTAzE71zJ9z4hALAAAAIEFABgAAABIEZAAAACBBQAYAAAASBGQAAAAgQUAGAAAAEgRkAAAAIEFABgAAABIEZAAAACBBQAYAAAASBGQAAAAgQUAGAAAAEgRkAAAAIEFABgAAABIEZAAAACBBQAYAAAASBGQAAAAgQUAGAAAAEgRkAAAAIFFoQDazD5pZycy+fZLlfszMpszsRXnNBgAAgNpU9B7kD0m6ar4FzKxO0jslbcljIAAAANS2QgOyu39N0sGTLPY7kv5RUmnxJwIAAECtK3oP8rzMbIOkX5T0/qJnAQAAQG2o6oAs6T2SXu/uU/MtZGY3mlmfmfXt27cvp9EAAACwHNUXPcBJ9Eq63cwkqUPS883suLt/Nl3I3W+TdJsk9fb2eu5TAgAAYNmo6oDs7k8/8bOZfUjS52eGYwAAAGAhFRqQzewTki6X1GFmA5LeLKlBktyd444BAACQu0IDsrtfn2HZly3iKAAAAICk6v+QHgAAAJArAjIAAACQICADAAAACQIyAAAAkCAgAwAAAAkCMgAAAJAgIAMAAAAJAjIAAACQICADAAAACQIyAAAAkCAgAwAAAAkCMgAAAJAgIAMAAAAJAjIAAACQICADAAAACQIyAAAAkCAgAwAAAAkCMgAAAJAgIAMAAAAJAjIAAACQICADAAAACQIyAAAAkCAgAwAAAAkCMgAAAJAgIAMAAAAJAjIAAACQICADAAAACQIyAAAAkCAgAwAAAAkCMgAAAJAgIAMAAAAJAjIAAACQqC96AABYLKVSKVTX2dm5wJNguVnOz63+/v5QXXd39wJPAhSHPcgAAGQQDccAlg4CMgAAAJAgIAMAAAAJAjIAAACQICADAAAACQIyAAAAkCAgAwAAAIlCA7KZfdDMSmb27Tlu/1Uze7Dyb5uZ/XDeMwIAAKC2FL0H+UOSrprn9sckXebuF0n6Q0m35TEUAAAAalehZ9Jz96+Z2bnz3L4tuXivpI2LPRMAAABqW9F7kLN4uaQvzHaDmd1oZn1m1rdv376cxwIAAMBysiQCspn9tMoB+fWz3e7ut7l7r7v3rlu3Lt/hAAAAsKwUeojFqTCziyR9QNLV7n6g6HkAoGilUilU19nZucCTzG/79u2hup6entx6tbS0ZK5pb28P9YqKPLbIOjwd3d3dufbr7+/Ppc/pPK7ojHmvS8yuqgOymW2S9GlJv+7ujxQ9D4ClJe9AWO1KpVJu6yQaWJeC4eHh3ENyVtu3b1dXV1du/aJhsNrXY39/f9XPiMVRaEA2s09IulxSh5kNSHqzpAZJcvf3S7pF0lpJ7zMzSTru7r3FTAsAAIBaUPS3WFx/kttfIekVOY0DAAAALI0P6QEAAAB5ISADAAAACQIyAAAAkCAgAwAAAAkCMgAAAJAgIAMAAAAJAjIAAACQICADAAAAiao+1TQA5K1UKoXqlspprZfCKYG3bNmSuWb9+vWLMAmAWkVABoAFUCqVcgvJ0T7RcBzR1dWl3bt359Zvz549Wrt2bea64eHhUL/R0dFQXWTbDQ4Ohnq1tLSE6g4dOhSqO3z4cOaaK664ItTrwQcfDNVF1snAwECo15lnnhmqQ3XgEAsAAAAgQUAGAAAAEgRkAAAAIEFABgAAABIEZAAAACBBQAYAAAASBGQAAAAgQUAGAAAAEgRkAAAAIEFABgAAABKcahoACrZ9+/bMNQcPHsxcs3Hjxsw1APK1bdu2zDWXXnrpIkxS2wjIAJDo7OwseoRFMzAwoMsuuyxU29/fn7lmZGQk1Gt0dDRUd+jQoVBdRENDQ6hux44dmWsmJiZCvYaGhkJ1q1evDtW1tLRkrtm6dWuoV/SxRera2tpCvYaHh0N1qA4cYgEAAAAkCMgAAABAgoAMAAAAJAjIAAAAQIKADAAAACQIyAAAAECCgAwAAAAkCMgAAABAgoAMAAAAJAjIAAAAQIKADAAAACTqix4AAKpNX19f5pre3t5FmAQAUAQCMgAsgL6+Pp1xxhm59RsaGgrV3XnnnaG6iYmJzDV79+4N9TKzUN3k5GTmmrGxsVCv5ubmUN3w8HDmmpaWllCvJ554IlTX2toaqotsN3cP9Vq9enWo7tChQ5lrxsfHQ706OjpCdRs3bgzVYWFxiAUAAACQICADAAAACQIyAAAAkCAgAwAAAIlCA7KZfdDMSmb27TluNzP7SzPbZWYPmtmP5D0jAAAAakvRe5A/JOmqeW6/WtLmyr8bJf11DjMBAACghhUakN39a5IOzrPItZI+4mX3SlpjZmflMx0AAABqUdF7kE9mg6Qnk8sDlet+gJndaGZ9Zta3b9++3IYDAADA8lPtAXm2bx1/yreKu/tt7t7r7r3r1q3LYSwAAAAsV9UekAcknZ1c3ihpT0GzAAAAoAZU+6mm75B0k5ndLunHJQ25e+zcpQCqQn9/f6iuu7t7gScBAGB2hQZkM/uEpMsldZjZgKQ3S2qQJHd/v6S7JD1f0i5Jo5L+RzGTAlgI0XCct5GRkVDdgQMHQnV1dXW51EjS448/Hqqbnp7OXDM2NhbqFV3/kXUS/dzKihWxP8BOTk7m1sv9KUcknpLItpZic65ZsybUa3BwMFR37NixzDXR9XHo0KFQ3fDwcOaaTZs2hXr19PSE6mpBoQHZ3a8/ye0u6VU5jQMAAABU/THIAAAAQK4IyAAAAECCgAwAAAAkCMgAAABAgoAMAAAAJAjIAAAAQIKADAAAACQIyAAAAECCgAwAAAAkCMgAAABAgoAMAAAAJOqLHgBA7eju7g7XlkqlBZxkfmNjY6G60dHRUF1jY2Pmmu9973uhXhMTE6E6d89cs3///lCvw4cPh+qmp6cz10S32dTUVKhufHw8t17Rx2ZmobqIyHP/dOoijy26/qOvtSNHjmSu2bRpU6hXT09PqK4WsAcZAAAASBCQAQAAgAQBGQAAAEgQkAEAAIAEARkAAABIEJABAACABAEZAAAASBCQAQAAgAQBGQAAAEgQkAEAAIAEARkAAABI1Bc9AICFUSqVwrWDg4Ohuq6urtx69fT0hOq2bt2auaa9vT3U6+jRo6G6/v7+zDVDQ0OhXitWxPaLHD58OHPN+Ph4qFe0bnR0NHPNoUOHQr2iM46MjGSuOXbsWKjX2NhYqG5ycjJUF1FXVxeqq6+PxZeWlpbMNe4e6hWti2y3tra2UK/zzz8/VBd9P15K2IMMAAAAJAjIAAAAQIKADAAAACQIyAAAAECCgAwAAAAkCMgAAABAgoAMAAAAJAjIAAAAQIKADAAAACQIyAAAAECCgAwAAAAkYiczB6pAqVQK13Z2di7gJPOLzpnnjAAA4PsIyMAycTqBenh4OLeeg4ODoV5bt24N1U1NTWWuefzxx0O9Dh06FKqLrP+xsbFQryNHjoTqJicnM9fs378/1Cv62CLrP7o+RkdHQ3WRbR3tNTExEaqLbOu8NTU1heoaGxsz1zQ0NIR6rVgR+yN95D1r3759oV6PPPJIqK6npydUt5RwiAUAAACQICADAAAACQIyAAAAkCg0IJvZVWb2sJntMrM3zHL7JjP7ipl9y8weNLPnFzEnAAAAakdhAdnM6iTdKulqSRdIut7MLpix2JskfdLdL5Z0naT35TslAAAAak2Re5AvkbTL3R919wlJt0u6dsYyLmlV5efVkvbkOB8AAABqUJEBeYOkJ5PLA5XrUm+R9GtmNiDpLkm/M9sdmdmNZtZnZn3RrzoBAAAApGIDss1ync+4fL2kD7n7RknPl/RRM3vKzO5+m7v3unvvunXrFmFUAAAA1IoiA/KApLOTyxv11EMoXi7pk5Lk7vdIapbUkct0AAAAqElFBuT7JG02s6ebWaPKH8K7Y8YyuyU9T5LM7HyVAzLHUAAAAGDRFBaQ3f24pJskbZG0U+Vvq9hhZm8zs2sqi/2epN80s/+Q9AlJL3P3mYdhAAAAAAumvsjm7n6Xyh++S6+7Jfn5IUk/kfdcWBo6OzuLHuGU5Dlnf39/qK69vT1Ut23btsw1+/fvD/WKGh0dzVxz6NChUK+DBw+G6oaGhjLXRB6XJI2MjITqjhw5krnm6NGjoV7RGSPrZGJiItRrbGwsVBfR0NAQqjOb7aM+1dUv2quxsTFU19LSkrmmtbU11Cv62CKfperq6gr1am5uDtXVAs6kBwAAACQIyAAAAECCgAwAAAAkCMgAAABAgoAMAAAAJAjIAAAAQIKADAAAACQIyAAAAECCgAwAAAAkCMgAAABAgoAMAAAAJOqLHgDAU919992huqNHj4bqpqamQnWHDx/OXDM5ORnqNTIyEqobGxvLXFMqlUK9oo9teno6c83o6GioV2R9RJlZqK6+PvarKVLX2NgY6rV69epQXWRbu3uoV/R1He3X1NSUuaauri7Uq6WlJVTX2tqauWbt2rWhXm1tbaG6rq6uzDWXXnppqNfP//zPh+pqAXuQAQAAgAQBGQAAAEgQkAEAAIAEARkAAABIEJABAACABAEZAAAASBCQAQAAgAQBGQAAAEgQkAEAAIAEARkAAABIEJABAACAROyE90CVKJVKmWs6OzsXYZK53X333bn2AwAAp4eAjKoQCbqn0ysSkrds2RLq19zcHKrbu3dvqC5ieno6t7ojR46Eeh0/fjxUd+zYscw1K1bE/rg2OTmZW11dXV2o18TERKguMmN0PUbrGhoacusVFXnNuHuoV7Suvj4WDSLPyZUrV4Z6RetWrVqVuaatrS3Uq6urK1T3jGc8I3PN+vXrQ736+/tDdd3d3aG6pYRDLAAAAIAEARkAAABIEJABAACABAEZAAAASBCQAQAAgAQBGQAAAEiEArKZtS70IAAAAEA1mDcgm9kGM+s1s8bK5U4ze7uk/8xlOgAAACBncwZkM/tfkh6Q9FeS7jWzGyTtlLRS0o/mMx4AAACQr/lOl3OjpGe6+0Ez2yRpl6Sfcvd78xkNAAAAyN98h1iMuftBSXL33ZIeIRwDAABguZtvD/JGM/vL5HJnetndX714Y6HWdHZ2Zq7Ztm1buN+uXbsy1zQ3N4d6PfHEE6G61tbsn4UdHx8P9RoZGQnVHT16NHNNdMahoaFQ3eTkZOaaurq6UK/p6elQ3fHjx0N1EfX1873tz83dM9dMTU2FekVfaxF5rntJWrEi+2fjo8/HqMiMUuw9K/p8bGhoCNW1tLRkrjnrrLNCvc4555xQXWTGiYmJUK/u7u5QXS2Y75n52hmX71/MQQAAAIBqMGdAdvcP5zkIAAAAUA0KPVGImV1lZg+b2S4ze8Mcy/yymT1kZjvM7P/lPSMAAABqS+zgnwVgZnWSbpV0paQBSfeZ2R3u/lCyzGZJN0v6CXc/ZGbZD1QFAAAAMihyD/Ilkna5+6PuPiHpdknXzljmNyXd6u6HJMndSznPCAAAgBqTKSCb2TcXsPcGSU8mlwcq16XOk3SemX3dzO41s6vmmOtGM+szs759+/Yt4IgAAACoNVn3INsC9p7tvmZ+j1C9pM2SLpd0vaQPmNmapxS53+buve7eu27dugUcEQAAALUma0C+cwF7D0g6O7m8UdKeWZb5nLtPuvtjkh5WOTADAAAAiyJTQHb3Ny1g7/skbTazp5tZo6TrJN0xY5nPSvppSTKzDpUPuXh0AWcAAAAAfkBhH9Jz9+OSbpK0RdJOSZ909x1m9jYzu6ay2BZJB8zsIUlfkfRadz9QzMQAAACoBYV9zZskuftdku6acd0tyc8u6TWVfwAAAMCiK/REIQAAAEC1mXMPspldK2mju99aufwNSSe+IuJ17v6pHOZDgUql2NdOd3bmcz6XSy+9VNu2bQvVTk5OZq5paGgI9VqxIvb/oSMjI5lrDhyIHYG0du3aUN13v/vdzDVHjhwJ9Yo6fvx45pr6+tgf16ampkJ1kefW6OhoqNf09HSorq6uLnNNY2NjqNfY2FiorqmpKXNN5HFJ8dd1ZJ2U/5ianVnsi6dWrlyZW11bW1uoV3S7rVnzlC/COqnoe3+0LvLY2tvbQ72q/fd8keZ7hb9OP/ihuSZJP6byV6799iLOBAAAABRmvt0kje6ensjj3yofkDtgZq2LPBcAAABQiPn2IJ+RXnD3m5KLnI0DAAAAy9J8AfkbZvabM680s1dK+vfFGwkAAAAoznyHWPyupM+a2UskfbNy3Y+qfCzyCxZ7MAAAAKAIcwZkdy9JutTMfkbShZWr73T3L+cyGQAAAFCAk36XUSUQE4oBAABQEzhRCAAAAJAgIAMAAAAJAjIAAACQICADAAAAiZN+SA+1K3qu9f7+/sw1o6OjoV5dXV2husceeyxzzeDgYKjXxMREqO7IkSOZa44dOxbqtXfv3lDd5ORk5pqmpqZQr+hji6z/qampUC93z60u+rwys1Bdnr0aGxsXeJK5tbS0hOoaGhpCddPT05lrVqyI7ctatWpVqC4qsi6j67+1NXZC3/b29sw169evD/WKvkaf+9znZq4ZHh4O9cLc2IMMAAAAJAjIAAAAQIKADAAAACQIyAAAAECCgAwAAAAkCMgAAABAgoAMAAAAJAjIAAAAQIKADAAAACQIyAAAAECCgAwAAAAkCMgAAOC0tLS0FD0CsKDqix4A1Wn79u3h2p6ensw1/f39oV6PPfZYqG50dDRzzfT0dKjX8ePHQ3VmlrnG3UO9pqamQnWRGaO9ous/MuPY2FioV0NDQ6iurq4uc000kETX/8TEROaa+vrYr5gVK2L7biKPLdqrubk5VBd5ja5cuTLUK1rX1NSUW11ra2uoV7Qu8pyMvq7PP//8UN2uXbsy13R1dYV6dXZ2hupqAXuQAQAAgAQBGQAAAEgQkAEAAIAEARkAAABIEJABAACABAEZAAAASBCQAQAAgAQBGQAAAEgQkAEAAIAEARkAAABIEJABAACARPaTkqMm9PT0qFQqhWq3bt2auSZ6HvmDBw+G6lpbWzPXnHXWWaFeo6OjobrGxsbMNW1tbaFeR44cCdW1t7dnrjlw4ECoV0tLS6ju6NGjmWs6OjpCvYaGhkJ1U1NTmWui62NiYiJU19zcnLlmeno61MvMQnX19dl/pR0/fjzUK7I+pNjruqGhIdQrsj4kacWK2L6zyDpx91CvVatWheoi6z/6WtuzZ0+o7uKLL85cMzg4GOqFuRW6B9nMrjKzh81sl5m9YZ7lXmRmbma9ec4HAACA2lNYQDazOkm3Srpa0gWSrjezC2ZZrl3SqyV9I98JAQAAUIuK3IN8iaRd7v6ou09Iul3StbMs94eS3iVpLM/hAAAAUJuKDMgbJD2ZXB6oXPdfzOxiSWe7++fnuyMzu9HM+sysb9++fQs/KQAAAGpGkQF5tk9h/NfR+ma2QtK7Jf3eye7I3W9z91537123bt0CjggAAIBaU2RAHpB0dnJ5o6T0I5/tkn5I0lfN7HFJz5Z0Bx/UAwAAwGIqMiDfJ2mzmT3dzBolXSfpjhM3uvuQu3e4+7nufq6keyVd4+59xYwLAACAWlBYQHb345JukrRF0k5Jn3T3HWb2NjO7pqi5AAAAUNsKPVGIu98l6a4Z190yx7KX5zETAAAAahunmgYAAAASBGQAAAAgQUAGAACnpbW1tegRTqqxsbHoEbCEFHoMMhZfqVTKvecVV1yRuWb79u2hXl1dXaG6/fv3Z66JnoRmamoqVNfU1JS5Znh4ONSrra0tVDc2lv0El6tXrw71mpiYCNVFfikeO3Ys1KujoyNUF3ls7n7yhWbR3NwcqjOb7avr51dXVxfqNT4+HqqLrv+I6HMksv6j67G+PvYrPlo3PT2duWbt2rWhXk972tNCdS0tLZlrRkdHQ70ivwujv7OjvwsxN/YgAwAAAAkCMgAAAJAgIAMAAAAJAjIAAACQICADAAAACQIyAAAAkCAgAwAAAAkCMgAAAJAgIAMAAAAJAjIAAACQ4FTTAAAACyx62mhUBwLyMtfZ2Rmu7e/vD9Xt3r07c80ZZ5yRWy9Jamtry1zT3Nwc6hV9bI888kjmms2bN4d6DQwMhOpWrlyZuWZ6ejrUa9WqVaG6+vrsb3PHjh0L9YrWNTY2Zq4xs1Cv8fHxUF1dXV3mmug2GxkZCdU1NDRkrmltbQ31ilq9enVuvaLrf3JyMlQXeT9YsSL2h+w1a9aE6pqamjLXdHV1hXpFnc7vbSwcDrEAAAAAEgRkAAAAIEFABgAAABIEZAAAACBBQAYAAAASBGQAAAAgQUAGAAAAEgRkAAAAIEFABgAAABIEZAAAACBBQAYAAAAS9UUPgOrV3d0dqiuVSgs8ydw2bdoUqhsfH8+t186dO0N1Z555Zuaa0dHRUK/169eH6oaGhjLXbNiwIdTrwIEDobqJiYlcaiSpo6MjVBfp5+6hXvX1sbf9tra2zDVjY2OhXmvXrg3VTU1NZa6JznjOOeeE6lasyL5f6tixY6Fe69atC9V1dXWF6h599NHMNatXrw716u3tDdVFfj/t3r071Cv6OxQIV6obAAAPEUlEQVTVgT3IAAAAQIKADAAAACQIyAAAAECCgAwAAAAkCMgAAABAgoAMAAAAJAjIAAAAQIKADAAAACQIyAAAAECCgAwAAAAkCMgAAABAgoAMAAAAJOqLbG5mV0n6C0l1kj7g7u+YcftrJL1C0nFJ+yT9hrs/kfugS1x/f3+orr29PVTX2dmZuaZUKoV6tbS0hOoiBgcHQ3VnnnnmAk8yt9WrV4fqhoaGQnXd3d2huoi2trZQXeQ5En3uj4+Ph+rq6/N7K44+Hx977LHMNc3NzaFe69evD9WNjo5mrjl69Gio14oVsf1L09PTmWs6OjpCvaIOHDgQqou8/0S3ddTu3bsz1zQ1NS3CJKh2he1BNrM6SbdKulrSBZKuN7MLZiz2LUm97n6RpE9Jele+UwIAAKDWFHmIxSWSdrn7o+4+Iel2SdemC7j7V9z9xC6BeyVtzHlGAAAA1JgiA/IGSU8mlwcq183l5ZK+MNsNZnajmfWZWd++ffsWcEQAAADUmiIDss1ync+6oNmvSeqV9Cez3e7ut7l7r7v3rlu3bgFHBAAAQK0p8kN6A5LOTi5vlLRn5kJmdoWkN0q6zN1jn4ABAAAATlGRe5Dvk7TZzJ5uZo2SrpN0R7qAmV0s6W8kXePusa85AAAAADIoLCC7+3FJN0naImmnpE+6+w4ze5uZXVNZ7E8ktUn6BzN7wMzumOPuAAAAgAVR6Pcgu/tdku6acd0tyc9X5D4UAAAAahpn0gMAAAASBGQAAAAgQUAGAAAAEoUeg4x8dHd359qvVMrvC0dGR0dPvtAC6erqyq2XJLW0tGSuaW9vD/XauXNnqG5ycjJzTeRxSdKFF14YqtuxY0fmmtWrV4d6RZ8jPT09mWu2bt0a6hVd/+edd17mmqNHj4Z6RUUeW3Rb5/nYmpqaQnXR5+P4eOwbVQ8fPpy5JvLcPx29vb2Za/L8nYbqwR5kAAAAIEFABgAAABIEZAAAACBBQAYAAAASBGQAAAAgQUAGAAAAEgRkAAAAIEFABgAAABIEZAAAACBBQAYAAAASBGQAAAAgQUAGAAAAEvVFDwBIUmdnZ251pVIp1Csqr8e2fft2jY6OhnpddtllobrIuhwcHAz1am9vD9VdeOGFmWu6u7tDvbZt2xaq6+vry1zT1dUV6tXS0hKqy1Oez5Hx8fFQr4aGhlBd5DWa9zbr7e0N1fX39y/wJNUh+h6OpY09yAAAAECCgAwAAAAkCMgAAABAgoAMAAAAJAjIAAAAQIKADAAAACQIyAAAAECCgAwAAAAkCMgAAABAgoAMAAAAJDjVNKpC9PTPnAIUAJauyHs/7/vIg7l70TMsqN7eXu/r6yt6DGQQDccSb5RL0fbt20N1PT09CzxJ7Tmd19rw8HDmmu7u7sw1/f39mWtOaG9vD9dmled7T/Q109XVtcCTzG9wcDBzTXRG3vsxHzO73917T+c+OMQCAAAASBCQAQAAgAQBGQAAAEgQkAEAAIAEARkAAABIEJABAACABAEZAAAASBCQAQAAgAQBGQAAAEgQkAEAAIBEoQHZzK4ys4fNbJeZvWGW25vM7O8rt3/DzM7Nf0oAAADUkvqiGptZnaRbJV0paUDSfWZ2h7s/lCz2ckmH3P0ZZnadpHdK+pX8pwVQtFKpFKobHh7OXNPd3R3qlafo+shTf39/0SMgJ4ODg0WPACyowgKypEsk7XL3RyXJzG6XdK2kNCBfK+ktlZ8/Jem9Zmbu7nkOisXV2dlZ9AjLwukEpjy3QU9PT6huKQTCpSC6rSN10YAc/R+USL+l8D9D0ddMVOS11tXVtQiTAMUp8hCLDZKeTC4PVK6bdRl3Py5pSNLamXdkZjeaWZ+Z9e3bt2+RxgUAAEAtKDIg2yzXzdwzfCrLyN1vc/ded+9dt27dggwHAACA2lRkQB6QdHZyeaOkPXMtY2b1klZLOpjLdAAAAKhJRQbk+yRtNrOnm1mjpOsk3TFjmTsk3VD5+UWSvszxxwAAAFhMhX1Iz92Pm9lNkrZIqpP0QXffYWZvk9Tn7ndI+jtJHzWzXSrvOb6uqHkBAABQG4r8Fgu5+12S7ppx3S3Jz2OSXpz3XAAAAKhdnEkPAAAASBCQAQAAgAQBGQAAAEgQkAEAAIBEoR/SAwDkK3IaYU4H/1TRU5+zLoGlgYAMLBOdnZ3hX9pLQTRYLNdAEnlc0edHqVTS8PBw5rru7u5Qv6i8+uX9Otu+fXuorqenJ1S3XF8zQBYcYgEAAAAkCMgAAABAgoAMAAAAJAjIAAAAQIKADAAAACQIyAAAAECCgAwAAAAkCMgAAABAgoAMAAAAJAjIAAAAQIKADAAAACTqix4AAKpNqVQK1XV2di7wJHOLzggAODkCMrCM5BnQlqvlHjyjz5GsdaVSKfd1mdfzv7OzU/39/eHarHp6ekK9AMRxiAUAAACQICADAAAACQIyAAAAkCAgAwAAAAkCMgAAAJAgIAMAAAAJAjIAAACQICADAAAACQIyAAAAkCAgAwAAAAkCMgAAAJCoL3oAAKgmnZ2dRY9wSqp9zmqf73R1d3dnrimVSiqVSpnrlvu6BKoRe5ABAACABAEZAAAASBCQAQAAgAQBGQAAAEgQkAEAAIAEARkAAABIEJABAACABAEZAAAASBQSkM3sTDP7opn9Z+W/Z8yyzLPM7B4z22FmD5rZrxQxKwAAAGpLUXuQ3yDpS+6+WdKXKpdnGpX0Une/UNJVkt5jZmtynBEAAAA1qKhTTV8r6fLKzx+W9FVJr08XcPdHkp/3mFlJ0jpJh/MZEQBOXeQUwics11MJn846yUt03S+FxwYgrqiA3OXueyXJ3fea2bzvUGZ2iaRGSf15DAcAwGJYrv8zBCw3ixaQzWyrpKfNctMbM97PWZI+KukGd5+eY5kbJd0oSZs2bco4KQAAAPB9ixaQ3f2KuW4zs0EzO6uy9/gsSbP+rcrMVkm6U9Kb3P3eeXrdJuk2Sert7fXTmxwAAAC1rKgP6d0h6YbKzzdI+tzMBcysUdJnJH3E3f8hx9kAAABQw4oKyO+QdKWZ/aekKyuXZWa9ZvaByjK/LOmnJL3MzB6o/HtWMeMCAACgVhTyIT13PyDpebNc3yfpFZWfPybpYzmPBgAAgBrHmfQAAACABAEZAAAASBCQAQAAgAQBGQAAAEgQkAEAAIAEARkAAABIFPI1bwAAoPqUSrOe2PakOjs7F3gSoFgEZABYAASEp1rO62Q5P7aIUqnEOsGywiEWAAAAQIKADAAAACQIyAAAAECCgAwAAAAkCMgAAABAgoAMAAAAJAjIAAAAQIKADAAAACQIyAAAAECCgAwAAAAkCMgAAABAwty96BkWlJkNS3q46Dlwyjok7S96CJwSttXSwbZaOthWSwfbaul4pru3n84d1C/UJFXkYXfvLXoInBoz62N7LQ1sq6WDbbV0sK2WDrbV0mFmfad7HxxiAQAAACQIyAAAAEBiOQbk24oeAJmwvZYOttXSwbZaOthWSwfbauk47W217D6kBwAAAJyO5bgHGQAAAAgjIAMAAACJZROQzezFZrbDzKbNrDe5/kozu9/Mtlf++zNFzom5t1Vy+yYzO2pmv1/EfPi++baVmV1kZvdUbt9uZs1FzYl53wMbzOzDlW2008xuLnJOSGb2J2b2HTN70Mw+Y2ZrkttuNrNdZvawmf1ckXOibK7tRb6oPvO9tiq3n3K+WDYBWdK3Jb1Q0tdmXL9f0i+4e4+kGyR9NO/B8BRzbasT3i3pC/mNg3nMuq3MrF7SxyT9lrtfKOlySZO5T4fUXK+rF0tqqrwH/qikV5rZufmOhhm+KOmH3P0iSY9IulmSzOwCSddJulDSVZLeZ2Z1hU2JE2bdXiJfVKO5ttUJp5wvls2JQtx9pySZ2czrv5Vc3CGp2cya3H08x/GQmGtbVa57gaRHJY3kPBZmMc+2+llJD7r7f1SWO5DzaJhhnm3lklor/1OzUtKEpCP5ToeUu/9LcvFeSS+q/HytpNsrv58eM7Ndki6RdE/OIyIx1/YiX1SfeV5bmfPFctqDfCp+SdK3ePJWJzNrlfR6SW8tehac1HmS3My2mNk3zex1RQ+EOX1K5V8IeyXtlvSn7n6w2JGQ+A19f4/WBklPJrcNVK5D9Ui3V4p8UX3+a1tF8sWS2oNsZlslPW2Wm97o7p87Se2Fkt6p8p4vLLLgtnqrpHe7+9HZ9i5jcQS3Vb2kn5T0Y5JGJX3JzO539y8t0phQeFtdImlK0npJZ0j6VzPb6u6PLtKY0KltKzN7o6Tjkj5+omyW5fku1hwEt9eJWvJFjoLbKnO+WFIB2d2viNSZ2UZJn5H0UnfvX9ipMJvgtvpxSS8ys3dJWiNp2szG3P29CzsdUsFtNSDpbnffL0lmdpekH5FEQF5EwW31Ekn/7O6Tkkpm9nVJvSr/qRGL5GTbysxukPTfJT3Pv39CggFJZyeLbZS0Z3EmRCq4vcgXBQhuq8z5YtkfYlH5BOOdkm52968XPQ/m5u7Pdfdz3f1cSe+R9HbCcdXaIukiM2upHNt6maSHCp4Js9st6WesrFXSsyV9p+CZapqZXaXyn3uvcffR5KY7JF1nZk1m9nRJmyX9exEz4vvm2l7ki+oz17aK5ItlE5DN7BfNbEDScyTdaWZbKjfdJOkZkv63mT1Q+ddZ2KCYb1uhysy1rdz9kKQ/l3SfpAckfdPd7yxuUszzurpVUpvK33Jxn6T/6+4PFjQmyt4rqV3SFyu/k94vSe6+Q9InVf6fzX+W9Cp3nypuTFTMur1EvqhGc22rzDjVNAAAAJBYNnuQAQAAgIVAQAYAAAASBGQAAAAgQUAGAAAAEgRkAAAAILGkThQCAJidmU1J2q7y+/pjkn7d3Q8XOxUALE3sQQaA5eGYuz/L3X9I0kFJryp6IABYqgjIALD83CNpQ9FDAMBSRUAGgGXEzOokPU/l0xYDAAIIyACwPKw0swckHZB0pqQvFjwPACxZBGQAWB6OufuzJJ0jqVEcgwwAYebuRc8AADhNZnbU3dsqP18s6XOSut19stjJAGDpYQ8yACwz7v4tSf8h6bqiZwGApYg9yAAAAECCPcgAAABAgoAMAAAAJAjIAAAAQIKADAAAACQIyAAAAECCgAwAAAAkCMgAAABA4v8DlWBaYmta7qgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x504 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import random\n",
    "n = 100\n",
    "x = r_absmag \n",
    "y = col_g_r  \n",
    "  \n",
    "\n",
    "# Creating bins \n",
    "x_min = -24  #np.min(x) \n",
    "x_max = -12  #np.max(x) \n",
    "  \n",
    "y_min = -0.3   #np.min(y) \n",
    "y_max = 1.5   #np.max(y) \n",
    "  \n",
    "nbin_x = round(abs(x_max-x_min)/0.25)\n",
    "nbin_y = round(abs(y_max-y_min)/0.05)\n",
    "\n",
    "x_bins = np.linspace(x_min, x_max, nbin_x) \n",
    "y_bins = np.linspace(y_min, y_max, nbin_y) \n",
    "  \n",
    "fig, ax = plt.subplots(figsize =(10, 7)) \n",
    "# Creating plot \n",
    "plt.hist2d(x, y, bins =[x_bins, y_bins], norm=mpl.colors.LogNorm(), cmap=mpl.cm.gray_r) \n",
    "plt.title(\"Color Magnitude Diagram\") \n",
    "  \n",
    "ax.set_xlabel('R')  \n",
    "ax.set_ylabel('G - R')  \n",
    "  \n",
    "ax.invert_xaxis()\n",
    "\n",
    "# show plot \n",
    "plt.tight_layout()  \n",
    "plt.show() \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Next steps: <font color=red> TO DO! </font> \n",
    "I would like Emma to add the points for the full Barton sample with small (blue?)points and then our specific sample of ~ 30 pairs with larger points (red?) to show how these compare with parent sample. Note, we might have to change the magnitude choices used since GRI may not be available for the Barton sample. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": 4
}