{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "78618b33-bf3f-4e7f-bf2d-99ad7486584e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from astropy import cosmology\n",
    "from astropy import units as u\n",
    "\n",
    "plt.rcParams.update({'font.size': 16})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3b0bd20f-2d9a-4254-a1ff-36fd5a997dd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up the universe with the parameters you want:\n",
    "universe=cosmology.LambdaCDM(H0=72, Om0=0.3, Ode0=0.7)\n",
    "\n",
    "# for full documentation of the cosmology package, see\n",
    "# https://docs.astropy.org/en/stable/cosmology/index.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b7c919e7-6d34-4c83-8ce9-62e2b9fa79c6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Redshift = 0.00\n",
      "  Luminosity distance is 0.00 Mpc\n",
      "  Angular diameter distance is 0.00 Mpc\n",
      "  Universe age is 13.09 Gyr\n",
      "  Lookback time is 0.00 Gyr\n",
      "Redshift = 0.10\n",
      "  Luminosity distance is 447.51 Mpc\n",
      "  Angular diameter distance is 369.85 Mpc\n",
      "  Universe age is 11.83 Gyr\n",
      "  Lookback time is 1.27 Gyr\n",
      "Redshift = 0.50\n",
      "  Luminosity distance is 2754.25 Mpc\n",
      "  Angular diameter distance is 1224.11 Mpc\n",
      "  Universe age is 8.19 Gyr\n",
      "  Lookback time is 4.90 Gyr\n"
     ]
    }
   ],
   "source": [
    "# set up some redshifts of interest\n",
    "zs = [0.0, 0.1, 0.5]\n",
    "\n",
    "# calculate interesting quantities at those redshifts\n",
    "for z in zs:\n",
    "    \n",
    "    print ('Redshift = {:.2f}'.format(z))\n",
    "    \n",
    "    ######################\n",
    "    # LUMINOSITY DISTANCE (DL) is what you use when calculating apparent magnitudes,\n",
    "    # i.e., m - M = 5*np.log10(DL) - 5 (and remember DL should have units of parsecs!)\n",
    "    DL=universe.luminosity_distance(z)\n",
    "    print('  Luminosity distance is {:.2f}'.format(DL))\n",
    "\n",
    "    # but before you stuff it into the mag-distance equation, make sure DL is in parsecs, \n",
    "    # and that you’ve converted it to a regular python number:\n",
    "    DL_pc = DL.to(u.pc).value\n",
    "    # and now you can use DL_pc to convert between apparent and absolute magnitudes\n",
    "\n",
    "    ######################\n",
    "    # ANGULAR DIAMETER DISTANCE (DA) is what you use when converting between \n",
    "    # angular size (alpha) and physical size (d_phys):\n",
    "    # d_phys = alpha[arcsec]*DA/206265\n",
    "    # or\n",
    "    # d_phys =  DA * np.tan(alpha[radians])\n",
    "    DA=universe.angular_diameter_distance(z)\n",
    "    print('  Angular diameter distance is {:.2f}'.format(DA))\n",
    "\n",
    "    # again, you can convert it to a regular number and unit, for example kpc\n",
    "    DA_kpc = DA.to(u.kpc).value\n",
    "    # and now you can use DA_kpc to convert between angular size and physical size\n",
    "    # in kpc\n",
    "\n",
    "    ######################\n",
    "    # AGES AND LOOKBACK TIMES can be worked out this way:\n",
    "    print('  Universe age is {:.2f}'.format(universe.age(z)))\n",
    "    print('  Lookback time is {:.2f}'.format(universe.lookback_time(z)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7db8b50b-8697-45f7-a8c3-89f1ead44c88",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Scale Factor')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEQCAYAAACeDyIUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzcklEQVR4nO3dd3hUZfbA8e+h916U3qsKYkTEBliwIKyLvYC6im0tq/52XVFUxI6KZUVRLIu6dgVFijQLigoWNJCQAKFIh0gggZByfn+8NzqOk2SSTM2cz/PME3Lvfe89l5uZM/e+TVQVY4wxpqyqRDsAY4wx8ckSiDHGmHKxBGKMMaZcLIEYY4wpF0sgxhhjyqVatAOIpGbNmmmHDh2iHYYxxsSNZcuW7VDV5oHWJVQC6dChA0uXLo12GMYYEzdEZF1x6+wRljHGmHKJeAIRkTYi8pSIfCUiOSKiItIhyLK1ROQREdksIvu8fRwf5pCNMcYEEI07kC7AuUAm8HkZy04FrgTGAcOAzcAcEekbygCNMcaULhp1IJ+paksAEbkCOCWYQiLSB7gQuFxVX/KWfQokA+OB4eEJ1xhjTCARvwNR1cJyFh0O5AFv+uwrH3gDGCoiNUMQnjHGmCDFUyV6b2Ctqub4LU8GauAejRljjImQeEogTXD1Jv52+az/ExEZIyJLRWTp9u3bwxacMcYkmnhKIAIEGnteSiqkqlNUNUlVk5o3D9gXxhhjKq2v1+zk2U9Xh2Xf8dSRcBfQLsDyxj7rjTHGADv35vLArBTeWbaRdk3qMOro9tSpEdqP/HhKIMnAWSJSx68epBdwAEiPTljGGBM7CguVt5Zu4IFZKWTn5nPtoM5cP6QrtWtUDfmx4imBzADuAc4BXgEQkWrAecBcVc2NYmzGGBN1KzdnMfb9n/hu/a/079iE+/5yCF1b1g/b8aKSQETkbO+fR3g/TxOR7cB2Vf1URNoDq4HxqjoeQFV/EJE3gUkiUh1YC1wDdAQuiuwZGGNM7MjOzWfSvFW8uDiDhrWrM/GcPozs1xqREquIKyxadyBv+/3+jPfzU2AQrmK8Kn+u5L8MuA+YADQCfgROVdXvwhWoMcbEKlVlTvJW7vkwmc2793NB/7b869QeNKpTIyLHj0oCUdXSWk5lEKB1laruA272XsYYk7A27MrhrhnJLEjZRo+D6vP0hf04on3j0guGUDzVgRhjTMI7kF/I85+v4akFaVQR4Y4zenLpwA5Uqxr5XhmWQIwxJk4sWbOTOz74mfRtezm190GMO7MXrRrVjlo8lkCMMSbG7dyby/0fp/Dudxtp07g2L16axJAeLaMdliUQY4yJVQWFyv++Wc8jc1LJOZDPdYM78/fB4enTUR6WQIwxJgZ9vz6TO6f/zM+/ZDGgUxMm/OUQurQIX5+O8rAEYowxMWTn3lwenp3Km0s30LJBTZ684HDOPOzgsPfpKA9LIMYYEwMKCpXXv1nPxDmpZOfmc9Xxnbj+xK7Uqxm7H9OxG5kxxiSI79ZnMs57XDWwc1PuGd47rEOQhIolEGOMiZKde3N5aHYKby3dSMsGNXnqgsMZFqOPqwKxBGKMMRFWUKi8/vU6r3VVQVw8rgokvqI1xpg45/u46pgu7nFVrLWuCpYlEGOMiQDfx1UHNajFfy7sx+mHHhQ3j6sCsQRijDFhVFCovPb1OiZ6j6uuPqEz1w/pQt04e1wVSPyfgTHGxKilGbu4a0YyyZuyOLZLM+4e3psuLepFO6yQsQRijDEhtmX3fh6ctZIPftjEwQ0rx+OqQCyBGGNMiOzPK2DqF2v5z8J08guV64d04ZpBnalTo3J+1FbOszLGmAhSVeat3MaEmStYtzOHob1bcscZvWjbpE60QwsrSyDGGFMB6dv2Mv6jFXy2ajtdW9Tj1b8dxbFdm0U7rIiwBGKMMeWQtT+PJ+el8fKXGdSuUZVxw3pxydHtqR6FmQGjxRKIMcaUQWGh8s6yjTw8J4Wd2Qc4/8i23HpKd5rWqxnt0CLOEogxxgRp2bpM7vkwmeUbd3NE+8a8fFl/DmndMNphRY0lEGOMKcW2rP08ODuF9777hZYNavLE+X0Z3qdVpWuWW1aWQIwxphi5+QW8tDiDp+ankVegXDuoM9cNrhy9yEPB/heMMSaABSlbGf/hCjJ25nBSz5bcOawn7ZvWjXZYMcUSiDHG+EjftocJM1eyKHU7nZrX5ZXL+3NCt+bRDismWQIxxhggM/sAT8xPY9qSddSpUZU7zujJqKM7UKNa4jTLLStLIMaYhJZXUMirS9YxaV4ae/bncdFR7bnppK4J2Sy3rCKeQESkLfA4cDIgwDzgJlVdH0TZdsC9wGCgGbAReAt4QFWzwxa0MabSUVUWpW7n3pkrWLM9m2O7NOPOYb3oflB8Tu4UDRFNICJSB1gA5AKjAQUmAAtF5LCSkoCI1MUlm+rAncB64EjgHqArcF54ozfGVBartrp6js9WbadTs7pMHZ3EkB4tEr5ZbllF+g7kSqAT0F1V0wFEZDmQBlwFPFZC2WNwiWKoqs71li0UkSbArSJSR1Vzwhe6MSbe7co+wKR5q3jt6/XUrVGVO4f14pIB7a2eo5winUCGA0uKkgeAqq4VkcXACEpOIDW8n1l+y38FquAehxljzJ8cyC9k2pJ1PDFvFdkHCrjoqHbcdFI3mtStUXphU6xIJ5DewPQAy5OBc0opOw93p/KQiFyDe4TVH7gReNbqQIwx/lSVBSnbuG/mStbsyOa4rq6eo1tLq+cIhUgnkCZAZoDlu4DGJRVU1f0icizwLi7hFHkB+Htx5URkDDAGoF27dmWN1xgTp1K37GHCzBV8nraDTs3r8tKlRzKoe3Or5wihaDTj1QDLSr2iIlILeBNoAVzC73cg44B84JqAB1OdAkwBSEpKCnRsY0wlsnNvLo/PW8XrX6+nfq3q3HVmLy4ekFjDrEdKpBNIJu4uxF9jAt+Z+PobMAjooqqrvWWfichuYIqIPKuqP4YsUmNMXDmQX8grX2bw5II0cg4UMOroDtx4YlcaWz1H2EQ6gSTj6kH89QJWlFL2UCDTJ3kU+cb72ROwBGJMglFV5q7YygMfryRjZw6DujfnjjN60qWF1XOEW6QTyAxgooh0UtU1ACLSAddE97ZSym4BGotIF99WXMBR3s9fQh2sMSa2Ld/4KxNmruSbtbvo3LwuL112JIO7t4h2WAlDVCNXLeB1BvwR2AfcgasPuReoDxymqnu97doDq4HxqjreW9YBWI5LJPfh6kCScJ0KVwH9VbWwpOMnJSXp0qVLQ39ixpiI2vTrPh6Zk8r73/9C07o1uOnkbpx/ZFur5wgDEVmmqkmB1kX0DkRVs0VkCG4ok2m4yvP5uKFM9vpsKkBVXP+OorIZIjIAuBvXe70ZsAFXQX5facnDGBP/9uzPY/Ki1Uz9Yi0KXDOoM9cM6kyDWtWjHVpCingrLG/Mq5GlbJNBgJZZqroCODc8kRljYlV+QSFvfLuBSfNWsWPvAf7StxW3Du1Om8Z1oh1aQrPReI0xMatowMP7P15J2ra99O/QhKmje9KnbaNoh2awBGKMiVErNmVx/8cr+SJ9Bx2a1uHZi49gaO+W1hEwhlgCMcbElK1Z+5k4J5V3vttIw9rVGTfMdQS0AQ9jjyUQY0xMyDmQz3OfrmHKZ2vILyzkimM78vfBXWlYxyrIY5UlEGNMVBUUKu8u28jEuals25PLGYcezL9O7UG7plZBHussgRhjouaLtB1MmLmClC17OLxdIyZf3I8j2gca7cjEolITiIhUB04Hlqvq2vCHZIyp7FZt3cP9H69kUep22japzdMXHs4Zhx5sFeRxptQEoqp5IvIWcCpgCcQYU27b97iRct/4Zj11a1bj9tN7MHpgB2pWqxrt0Ew5BPsIaw1uGHVjjCmzfQcKeHHxWiYvWs3+PDdS7g0ndrUZAeNcsAnkYWCsiCxQ1e3hDMgYU3kUVZA/+kkqW7NyOaVXS247rQedmteLdmgmBIJNIENw83isFZElwGb+ODGUquroUAdnjIlPRT3IH5yVQupWV0H+9IX9OLKDVZBXJsEmkGOBPGA70Nl7+bKZ/owxgBti/YGPU/hqzU46NK3D5Iv6ceohB1kFeSUUVAJR1Y7hDsQYE9827MrhkTmpzPhxE03r1mD8iN5c0L+dDbFeiVk/EGNMhWRmH+DphelM+2odVarA9UO6MOb4TtS3IdYrvaATiIjUAS4HTsDVh+wEFgEvq2pOWKIzxsSs/XkFvPxlBv9ZmE52bj7nJrXlppO6cVDDWtEOzURIUAlERA7CJYtuwDrcrICdgLOB60VkkKpuDVeQxpjYUVCofPD9Lzw6N5VNu/dzYo8W/Ou0HnRraXOQJ5qyNONtDBynqouLForIQOBd4CHg0pBHZ4yJKZ+t2s4Ds1JYuTmLw9o0ZOK5fRjYuVm0wzJREmwCOQ34l2/yAFDVL0XkDuDBkEdmjIkZyZt28+CsFD5P20HbJrV58oLDGXbowVSpYi2rElmwCaQesKmYdRu99caYSmZjZg6PzV3F+z/8QsPa1blzWC8uHtDOhh4xQPAJJBW4BJgdYN3FQErIIjLGRN3unDyeWZTOS19mIMDVJ3Tm6hM607C2tawyvws2gUwE/isiLYHXcT3RDwLOB07CJRdjTJzLzS9g2lfreGpBOln78xjZrw03n9yNVo1qRzs0E4OC7Uj4qteMdzzwgs+qrcDVqvp6OIIzxkRGYaHy4fJNPDInlY2Z+zihW3NuO60HPQ9uEO3QTAwLuh+Iqk4RkReA7rh+ILuAVFUtDFdwxpjw+zJ9B/fPWsnPv2TRu1UDHvzrYRzb1VpWmdIF2w9kHPCCqm4CVvqtOxi4UlXHhyE+Y0yYpGzJ4sFZKSxK3U7rRrWZdF5fhvdpZS2rTNCCvQO5C1eBHqglVitvvSUQY+LA5t37eGzuKt75biP1a1Zj7Ok9ueTo9tSqbi2rTNkEm0BK+krSGMgNQSzGmDDK2p/Hs4tWM/WLtajCFcd25LrBXWhUxyZ1MuVTbAIRkUG4eUCKXCUiw/w2qw2cASSHPDJjTEgcyC/kta/X8eT8NDJz8vhL31bcckp32japE+3QTJwr6Q7kBOAO798KXBZgmwPACuCGYA8oIm2Bx4GTcXc284CbVHV9kOV74h6XDQbqAuuBZ1T1iWBjMCYRqCozf9rMw7NTWb8rh2O6NOXfp/XkkNYNox2aqSSKTSCqeg9wD4CIFAJHq+rXFTmY1xR4Ae6R12hcYpoALBSRw1Q1u5TySV75RcAVwG6gK9YT3pg/+HrNTu7/eCU/btxNj4Pq88rl/Tm+azOb1MmEVLD9QEI1I8yVuFF8u6tqOoCILAfSgKuAx4orKCJVgFeA+ap6ls+qhSGKzZi4t2rrHh6alcL8lG0c3LAWE8/pw1mHt6aqtawyYRBsM97LgPaqeneAdXcDa1X1lSB2NRxYUpQ8AFR1rYgsBkZQQgIBBgG9gKuDidmYRLJl934e/2QVby/bQN0a1fjnqd25/JiO1rLKhFWwrbBuBKYWs24bcBPu7qA0vYHpAZYnA+eUUvZY72ctEVkCHAFkAm/gRgreF8TxjalUsvbn8dynrmVVQaFy2TEd+fvgLjSuay2rTPgFm0C6UHxLq5VA5yD30wT3oe9vF645cElaeT/fBJ4GbgOScBXqbYGzAhUSkTHAGIB27doFGaYxsc2/ZdWIvq241VpWmQgLNoHkA8WNbdC8jMfUAMuCeUBbVA/zqqqO8/69SESqAg+KSC9VXfGng6lOAaYAJCUlBTq2MXGjsNC1rHpkjmtZNbCza1l1aBtrWWUiL9gE8g2u7uGtAOuuBr4Ncj+ZuLsQf40JfGfia6f38xO/5XNxE1r1xTUpNqZS+nL1Dh6clcJyr2XVy5cdyQndmlvLKhM1wSaQ+4B5IvI1bjTeX4DWuKa0/XB9OoKRjKsH8deL0j/8ix6h+d9FFL17bFBHUymlbMnioVkpLEzdTitrWWViSLDNeD8VkbOBScBzPqsygJGquijI480AJopIJ1VdAyAiHYBjcHUaJZmF6z9yKvCRz/Kh3s+lQcZgTFzwHbOqXs1q/Pu0Howe2MFaVpmYUZbh3KcD00WkO9AU2KGqq8p4vOeBv3v7uQN3N3EvsAGfxCQi7YHVwPiiUX5VdaeIPADcKSJZuA6FScA44BXfpsHGxLPd+/J49tPVvGhjVpkYF3QCKaKqqeU9mKpmi8gQ3FAm03CPn+bjhjLZ67OpAFX5veK8yHhgD3AtcCtuZsRHcEnImLiWm1/Aq0vW89SCNH7NyeOsw1tz88ndrGWViVllSiAi0gc3oVQt/3Wq+t9g9uGNeTWylG0yCNAyS1UV19mwpA6HxsSVotkAJ85NZcOufRzbpRm3ndbDxqwyMS/YnuiNgJnAgKJF3k/fCu2gEogx5ndfpu/ggVkp/PTLbnoe3ID/Xn4ox3cra8t4Y6Ij2DuQ+3H1HscDn+M67e0GLgeOBs4PS3TGVFIrN7vZAD9d5WYDfOzcPvylb2ubDdDElWATyFDcyLxLvN83quoyXEe+ybihTkaFIT5jKpVNv+7j0bmreO97Nxvg7af3YNTR1rLKxKdgE8jBwBpVLRCR/UB9n3Xv4cajMsYUY/e+PJ5ZlM5LizMAuPK4Tlw7qLO1rDJxLdgEsgVo5P17He6x1SLv9y6hDcmYyiM3v4BpX63jqQXpZO3P46y+rbn5lG60aWwtq0z8CzaBfIFLGh/hmt/e5XUAzMdNDDUjLNEZE6cKC5UZP27ikTmp/PLrPo7v1pzbTu1Br1YNoh2aMSETbAK5h99Hw30EV6F+HlAHlzyuD31oxsSnL9J28MCslSRvyqJ3qwY8NPIwju1a3FikxsSvYhOIiIwHnlXVTaq6GlgtIlVUNQ+4xXsZYzzJm3bz4KwUPk/bQetGtZl0Xl+G92llLatMpVXSHchY3COrTQDesOkHRORIVf0uEsEZEw82797HxDmuZVWDWtW544yeXDygvbWsMpVeSQkk0Ncm+ypljGfPfjdm1QufuzGrxhzXiWsHdaFhnerRDs2YiCjzWFjGJLq8gkLe+HYDkz5Zxc7sAzYboElYlkCMCZKqMm/lNh6ctZLV27Pp37EJL57ekz5tG0U7NGOiorQEcqaIHOL9uwpu7KvhItLXf0NVfTHEsRkTM5Zv/JX7Zq7k67W76NS8Ls+PSuKkni1sNkCT0EpLIGMDLBsXYJkClkBMpbMxM4eJc1L54IdNNKlbg3tH9Ob8/u2oXtV/pgFjEk9JCaRjxKIwJsZk7c/jmYWreXHxWgS4dlBnrh7UmQa1rILcmCLFJhBVXRfJQIyJBXkFhby2ZB1PzE8jMyePv/Zrza2ndKdVo9rRDs2YmGOV6MbgKsjnJG/lodkprN2RzdGdmjL2jJ42qZMxJbAEYhLe9+szuf/jlXybkUmXFvV48dIkBne3CnJjSmMJxCSsDbtyeHhOKh/+uIlm9Wpw31mHcF5SW6pZBbkxQbEEYhLO7n15/GdhOi8vzqBKFbhhSBfGnNCZejXt7WBMWdg7xiSM/IJC/vfNeh6fl0ZmzgFG9mvDrad056CGtaIdmjFxqUwJRESqAL1ww7kvVdXssERlTIh9umo7Ez5aQdq2vQzo1IQ7zuhlFeTGVFDQCURErgPuwiUPgCOB70TkA2CBqj4Z+vCMqZj0bXuYMHMli1K3075pHZ675AhO6dXSKsiNCYGgEoiIXAk8gettPhd4y2f158BIwBKIiRmZ2QeYNG8Vr369njo1qjL29J6MGtiemtVsiHVjQiXYO5CbgUdV9V/evCC+UoD/C21YxpTPgfxC/vtVBk/OT2Nvbj4XHtWOf5zUjab1akY7NGMqnWATSEdgTjHrsoFGIYnGmHJSVT5ZsZUHZrmOgMd1bcadw3rRrWX9aIdmTKUVbALZAXQoZl134JeQRGNMOazYlMWEmSv4cvVOOjevy0uXHsmg7s2tnsOYMAu2x9SHwDgR6eSzTEWkGfAP4INgDygibUXkHRHZLSJZIvKeiLQLPuTf9vNvEVER+aKsZU3lsH1PLre9u5wznvqcFZuzuGd4b2bfdDyDe1gvcmMiIdg7kDuAIcDPwNe44dufBHoA24DxwexEROoAC4BcYLS3nwnAQhE5LNhmwV4iG+sd2ySY3PwCXvwig/8sTGd/XgGXDezIjSd2talkjYmwoBKIqu4UkSTgJmAosNor+zTwuKpmBXm8K4FOQHdVTQcQkeVAGnAV8FiQ+5kMvIZ7fGadIROEqjJ/5TbunbmCdTtzOKlnC24/vSedmteLdmjGJKSgP3xVdQ9wr/cqr+HAkqLk4e13rYgsBkYQRAIRkQuBfsAFwHsViMXEkfRtexn/0Qo+W7Wdzs3r8t/L+3N8t+bRDsuYhBbpb++9gekBlicD55RWWEQaA48D/1TVXfacu/LL2p/Hk/PSePnLDGpXr8qdw3ox6uj2NiOgMTGg2AQiIgvKsB9V1ROD2K4JkBlg+S6gcRDlHwFWAS8HG5iIjAHGALRrV+a6ehMlhYXKO8s28vCcFHZmH+C8pLbcOrQ7zaw/hzExo6Q7kCq4Su5glOVWINA+Sy0vIscBo4B+qhpsXKjqFGAKQFJSUtDlTPQsW5fJPR8ms3zjbvq1a8RLl/bn0DY2bpUxsaakKW0HheF4mbi7EH+NCXxn4us5YCqwUUQaecuqAVW93/epam6I4jRRsDVrPw/NSuG973+hRf2aTDqvLyP6trImucbEqEjXgSTj6kH89QJWlFK2p/e6OsC6TFx/lEkVCc5ER1Gz3KcXpJFXoFwzqDPXDe5i83MYE+PKOpx7Y6Ar8KcJFFT1syB2MQOYKCKdVHWNt88OwDHAbaWUHRxg2SSgKnA9kB5gvYlxi1K3cfeMZDK8Zrl3nNGLDs3qRjssY0wQgh2NtxZuJN5zKb6+IphhTp8H/g5MF5E7cPUh9wIbcI+oio7XHtfXZLyqjgdQ1UUB4voVqBZonYltv/y6j/EfJjMneSudmtXl5cuOZFD3FtEOyxhTBsHegdwJDML1Hp8GXAfsBy4FDgZuDGYnqpotIkNwTXGn4ZLRfOAmVd3rs6ngEpK11axkDuQX8vzna3hqQRoA/ze0O1cc19GGWTcmDkkwDZpEJAX3uOh5IA9IUtXvvHVvA5tUNagkEk1JSUm6dOnSaIeRsBan7+DO6T+zZns2Q3u35M5hvWjTuE60wzLGlEBElqlqUqB1wd6BtAOSVbVARPIA34fULwIvEeRdiEk8W3bv596ZK5i5fDPtm9bhpcuOZLA9rjIm7gWbQHYCRQMObQD64GYiBGgG1A5xXKYSyCso5OXFGUyat4r8QuUfJ3XjqhM6Uau6Pa4ypjIINoEsAQ4HZgHvAveKSH0gH7gFsCHVzR8sWbOTcdN/ZtXWvQzp0YK7z+xNu6b2uMqYyiTYBPIQ7jEWuOHXu+CGcK+KSy7XhD40E4+27dnP/TNX8sEPm2jdqDbPj0ri5F4tox2WMSYMgh3OfSmw1Pv3HmCkiNQEapZhKHdTiRUWKm98u4EHZ61kf14h1w/pwrWDulC7hj2uMqayKndXX2/YEBs6xJC6ZQ+3v/8Ty9ZlcnSnpkw46xA62xwdxlR6wXYkfBxopqqXBFg3DdiqqreGOjgT2/bnFfDUgjSe+3QN9WtVY+I5fRjZr7WNXWVMggj2DmQ4cHcx6+Z46yyBJJDP07Zzxwc/s25nDiP7tWHsGT1pUrdGtMMyxkRQsAmkNa75biAbvfUmAezYm8uEj1bwwQ+b6NisLq9feRQDOzeLdljGmCgINoFk4lpeLQqwrguwJ1QBmdikqry1dAP3f5xCzoF8bjixK9cO6mx9OoxJYMEmkHnAWBH5UFW3Fi0UkZbA7cAn4QjOxIb1O3O47b3lfLl6J/07NOH+vx5Clxb1ox2WMSbKyjKY4rdAmoh8xO+PrYbhWmLdEZ7wTDQVFCovf5nBxDmpVK0i3H/WoZx/ZFuqVLFKcmNM8P1AMkTkSFznwZOBpsAO4H3gLlVdF74QTTSkbd3DP99dzvfrf2VIjxbcd9YhHNzQRqwxxvwu6H4gqpqBm5PcVGJ5BYVMXrSapxekU7dmVZ44vy/D+9i0ssaYPytXR0IRaYibmXCLqm4MbUgmWn7auJv/e+dHUrbsYdhhB3P38N40q1cz2mEZY2JUsQlERIYCg1X1Nr/lY4FxRWVF5E1glKrmhzNQEz65+QVMmpfGlM/W0LRuDaZccgSn9D4o2mEZY2JcSXcgV+OmnP2NiJyMm4L2J+AFoCdwFbAMeDRMMZow+vmX3dz81g+s2rqXc5PaMPaMXjSsXT3aYRlj4kBJCeRwXLLwdRluKtuhqroFKHo2fiGWQOJKXkEh/1mYztML0mlSt4ZN8mSMKbOSEkgLYLXfspOBL4qSh2cm8KcxskzsWrV1Dze/9QM//5LFX/q24u7hvWlUx4YhMcaUTUkJZA8+U9eKSFdc890lfttl4eYFMTGuoFB54fM1PDp3FfVqVWPyRf047dCDox2WMSZOlZRAUoARuDsMvH8rMNdvu47AVkxMy9iRzS1v/8iydZkM7d2S+8461FpYGWMqpKQE8jjwnog0wSWIS3GV54v9tjsL+DEs0ZkKKxrD6u4ZK6heVZh0Xl9G9LV+HcaYiis2gajqByJyE27O8ya4R1dXq+pvLbNEpA0wGPhnmOM05ZCZfYB/v/cTs5O3MLBzUx49t4/1JjfGhEyJHQlV9UngyRLWbwQahTgmEwKL03dw81s/sCv7ALef3oMrju1kY1gZY0Kq3FPamtiUm1/Ao3NXMeWzNXRuXpepo4/kkNYNox2WMaYSsgRSiaRv28MN//uBFZuzuHhAO8ae3ovaNayBnDEmPCyBVAKqytvLNjJu+s/UqVGNF0YlcVKvltEOyxhTyVWJ9AFFpK2IvCMiu0UkS0TeE5F2QZRLEpEpIpIiIjkisl5EXhORjpGIO1Zl5+Zzy1s/8s93lnN428bMvvE4Sx7GmIiI6B2IiNQBFuAmoRqN61cyAVgoIoepanYJxc8HeuMq9ZNxE1rdCSwVkb6qWtyc7ZVWypYsrnvtO9bsyOamk7py/ZCuVLWKcmNMhET6EdaVQCegu6qmA4jIciANNyjjYyWUfUhVt/suEJHFwFpvv+PCEnEMUlXe/HYDd81IpkHt6rx2xVEM7Nws2mEZYxJMpBPIcGBJUfIAUNW1XiIYQQkJxD95eMvWich23N1IQsjOzef2939i+g+bOK5rMx47ty/N61uPcmNM5EU6gfQGpgdYngycU9adiUhP3KCPKysYV1xYuyObq6YtJX3bXm45uRvXDe5ifTuMMVET6QTSBMgMsHwX0LgsOxKRasCzwHZgagnbjQHGALRrV2pdfcyat2Ir/3jzB6pVFf57+VEc29UeWRljoisazXg1wLLyfI1+GhgInKGqgZKSO5jqFGAKQFJSUqBjx7TCQmXS/DSenJ/GIa0bMPmiI2jbpE60wzLGmIgnkEzcXYi/xgS+MwlIRB7A3VWMVlX/0YErjd05edz05vcsTN3OyH5tuO+sQ6hV3ToGGmNiQ6QTSDKuHsRfL2BFMDvw5mS/DbhBVaeFMLaYkrplD2OmLeWXzH3cO6I3Fw9obyPoGmNiSqQ7Es4ABohIp6IFItIBOMZbVyIRuQHXb2Ssqj4VriCjbWHKNkZO/pKcAwW8edUALjm6gyUPY0zMiXQCeR7IAKaLyAgRGY5rlbUBeK5oIxFpLyL5IjLOZ9n5wCRgNrBARAb4vHpF8iTCRVWZ+sVa/vbKt7RvWocZfz+GI9oHeuJnjDHRF9FHWKqaLSJDcJNVTcNVns8HblLVvT6bCm6aXN8Ed6q3/FTv5etTYFCYwo6IvIJC7pqRzOtfr2do75Y8fl5f6tSwocqMMbEr4p9QqroeGFnKNhn4tcxS1UtxsyJWOrtz8rj29WUsTt/JtYM6c+sp3a1/hzEm5tlX3CjbmJnD6Be/Yf2uHCae04ezj2gT7ZCMMSYolkCiaOXmLEa/+A378wp49W9HcVSnptEOyRhjgmYJJEq+Wr2TMf9dSr1a1XjnmoF0a1k/2iEZY0yZWAKJgpnLN/OPN3+gfdM6vHJ5f1o1qh3tkIwxpswsgUTY61+vZ+wHP3FEu8a8MDqJRnVqRDskY4wpF0sgEfTS4rXc8+EKBndvzuSLj7BhSYwxcc0SSIRMXrSah2anMLR3S566oB81qkV8NmFjjAkpSyBhpqo8MT+NSfPSOLNPKx47tw/Vq1ryMMbEP0sgYfb4J6t4ckE6Zx/RhodGHmZzlhtjKg37KhxGzyxK58kF6Zyb1IaHLXkYYyoZSyBh8vLitTw8O5XhfVrxwF8Ps6FJjDGVjiWQMHhr6Qbu/nAFJ/dqyaPn9rE7D2NMpWQJJMTmr9zKbe8u57iuzXj6wsOtwtwYU2nZp1sILd/4K39//Xt6t2rIsxcfQc1q1s/DGFN5WQIJkQ27crj85aU0qVuDqZcmUbemNXAzxlRu9ikXAntz87n85W85kF/AG2OOokX9WtEOyRhjws4SSAWpKre+9SNrdmQz7fL+dGlho+oaYxKDPcKqoMmfrmZ28hb+fVoPBnZpFu1wjDEmYiyBVMDi9B1MnJPKsMMO5m/Hdox2OMYYE1GWQMopM/sAN7/1Ax2a1eXhsw9DxPp6GGMSiyWQclBVbn//J3ZlH+DJ8w+nTg2rSjLGJB5LIOUw/YdNzPp5C7ec0p1DWjeMdjjGGBMVlkDK6NecA9z70Qr6tm3Elcd1inY4xhgTNfbspYwemp3Kr/vymHbWoTbGlTEmodkdSBms2JTF/75Zz6UDO9CrVYNoh2OMMVFlCaQMHp2bSoNa1bhhSNdoh2KMMVFnCSRI36/PZH7KNq46oTMN61SPdjjGGBN1EU8gItJWRN4Rkd0ikiUi74lIuyDL1hKRR0Rks4jsE5GvROT4cMcM8OLiDOrXqsalAztE4nDGGBPzIppARKQOsADoAYwGLgG6AgtFpG4Qu5gKXAmMA4YBm4E5ItI3LAF7tu/JZdZPmzk3qa2NsmuMMZ5IfxpeCXQCuqtqOoCILAfSgKuAx4orKCJ9gAuBy1X1JW/Zp0AyMB4YHq6g567YQn6hcm5S23Adwhhj4k6kH2ENB5YUJQ8AVV0LLAZGBFE2D3jTp2w+8AYwVERqhj5cZ8HKbbRvWoduLeuF6xDGGBN3Ip1AegM/B1ieDPQKouxaVc0JULYG0KXi4QX286bdHNG+sY13ZYwxPiKdQJoAmQGW7wIaV6Bs0fo/EZExIrJURJZu37496ECL5BUUcmyX5hzftXmZyxpjTGUWjRphDbAsmK/2Up6yqjoFmAKQlJQUqHyJqletwqPn9ilrMWOMqfQifQeSSeA7hcYEvrvwtauEskXrjTHGREikE0gyri7DXy9gRRBlO3pNgf3LHgDS/1zEGGNMuEQ6gcwABojIb8PYikgH4BhvXWllqwPn+JStBpwHzFXV3JBHa4wxpliRTiDPAxnAdBEZISLDgenABuC5oo1EpL2I5IvIuKJlqvoDrgnvJBG5QkROxDXh7QjcFblTMMYYAxFOIKqaDQwBVgHTgNeAtcAQVd3rs6kAVQPEdxnwEjABmAm0BU5V1e/CHLoxxhg/EW+FparrgZGlbJNBgNZVqroPuNl7GWOMiSIbjdcYY0y5WAIxxhhTLqJa5r51cUtEtgPrylm8GbAjhOHEskQ510Q5T7BzrYwidZ7tVTXgUBwJlUAqQkSWqmpStOOIhEQ510Q5T7BzrYxi4TztEZYxxphysQRijDGmXCyBBG9KtAOIoEQ510Q5T7BzrYyifp5WB2KMMaZc7A7EGGNMuVgCMcYYUy4JnUBEpK2IvCMiu0UkS0TeE5F2QZatJSKPiMhmEdknIl+JyPHhjrk8RORsEXlXRNZ5saaKyAMiUj+IslrMq28EQi8zERlUTLy/BlE2bq4pgIgsKuH6zC6lbMxeVxFpIyJPef//OV5cHQJs11hEXhCRHSKSLSLzROTQII9RRUT+LSIZIrJfRH4UkRKHWAq1YM5TRE4UkVdFZLX3N7laRCaLSIsgj5FRzHX+SyjOIRozEsYEb16RBUAuMBo32+EEYKGIHOYN/FiSqcAZwP8Ba4DrgDkicrQ3cnAsuRVYD9wObAQOB+4GBovIQFUtLKX8y/iMluxZFeIYQ+0G4Fuf3/ODKBNP1xTgWqCB37KjgccofXoEiN3r2gU4F1gGfA6c4r+BiAjuHDsC1+MmpPs37v3bV1U3lnKMe3Hvi7Hecc4H3haRYar6cahOpBSlnidwNVAP99m0BugK3AMM9T6n9gYo428O7v3uK7WcMf+RqibkC7gRKAC6+CzriPugubmUsn1wCecyn2XVvIsyI9rnFiDe5gGWjfLOYUgpZRWYEO1zKMO5DvJiPqmM5eLqmpZwHlNxX4qaxOt1Bar4/PsKL9YOftuM8JYP9lnWEDcz6ZOl7L+F9390j9/y+cDyGDvPQO/d471tLw/iGBnAq+E6h0R+hDUcWKKqv81kqKprgcW4P87Syubh5icpKpuPm59kqIjUDH245aeq2wMsLvp23jqSscSwuLqmgYhIbdyEax+qatxO8ayl3xGDu16bVHWhT7ndwIeU/v4dCtQAXvVb/ipwqIh0LEO45RbMecb6ezeRE0hv4OcAy5Nx0+SWVnatquYEKFsDd2sa607wfq4MYttrRCTXe067QESOC2dgIfKaiBSIyE4ReT2Iuq3KcE3/CtQHXgly+3i8rkVKev+2E5F6pZTN5c/TYCd7P0t7/0dbWd67AGd61zhXRJaEqv4DEjuBNME9N/W3C2hcgbJF62OWiLQGxgPzVHVpKZu/invWfhIwBmgKLBCRQeGMsQJ2A4/iHgkMwT3rPgn4qpSKx7i+pp5RwDZgVhDbxtt19Vfa9SrpPdwE+FW9ZzwBysbstfYavkzCJY8PgijyIa6OaChwEbAfeF9ELg5FPAlbie4J1IvyTxNZFbNNectGlffNbDquruey0rZX1Ut8fv1cRKbjvvlNAI4NS5AVoKrfA9/7LPpURD4DvsFVrN9RTNG4vaYAItIKlwye8B69lSjermsAFblecXmtRaQa8D/co6tjgrzO1/vt431gCfAAf36EV2aJfAeSSeBvGo0J/M3G164SyhatjzkiUgvXcqUTMFRLb6nyJ6q6Bzed8JEhDi9s1E15vIqSY47La+rjYtz7OdjHV38Qh9e1tOtV0nt4F9DYa8kVqGzMXWsRKbq2JwF/UdXl5dmPqhYAbwNtROTgisaVyAkkGfcs1F8vYEUQZTt6TYH9yx7gz89Wo05EqgPvAv2B01X1p4rsjsDf4GJZaTHH3TX1Mwr4UVV/rMA+4um6lvT+Xa8lN29NBmoCnQOUhdLf/9HwLHAecL6qzq/gvooSZ4WvdSInkBnAABHpVLTA68RzDKW3oZ8BVMe1eCkqWw13geeqam7Io60A79vLa8CJwAhVXVKBfTXA9ZX4OkThhZ2IJAHdKDnmuLqmvrzz60057z68fcTbdZ0BtBaRogrlonM4k9Lfv7NxXwou8lt+MfCz1xozZohIUZ3eZar6QQX3VQ33N75eVbdUOLhItXmOtRdQF/et8idcs7/hwI+4zjr1fLZrj6svGOdX/g3cbfIVuA/md3AVVP2ifW4BznUyv3eUHOD3alPceeI6Wj0PXIjrXzHa+/86ABwX7fMq5lxf887zr7hK9Ftws7atB5pVlmvqF/eTuCbILQOsi8vrCpztvYr+dq/xfj/BW18F+BLYgOsEOBRYhHv81NZvX/nAVL9lD3rX9mbv/2AyUAicGWPn+S9v+dQA793OJZ0ncIH3Nz0KGOz9P33u7e/8kMQf7T+UKP+RtsM91skC9uBaNXTw26aD9x9+t9/y2rgev1u8P8SvgUHRPqdizjPDO4dAr7uLO0/ct7nFuA/gPGAn7ttd/2ifUwnn+m9gOa41Vp73ATMFOLgyXVOfmKsD23F9PwKtj8vrWsLf6yKfbZoAL+KSRg6uI2CfYvb1st+yqrgGFetwTXqXA2fH2nnikmJx2/if0x+W4ZLMAmCrd513A/Nw9Z8hid+GczfGGFMuiVwHYowxpgIsgRhjjCkXSyDGGGPKxRKIMcaYcrEEYowxplwsgRhjjCkXSyAmYXnToaqIPBbFGC71nWo0wPp6InKbiHwrbtrlAyKyXkTeFpERAcZzqkgsF/tNe5rog62aUlgCMQnJZ/IlgIti4MPyr7jpaH/jDbv/Da438ie4YVVOBu7CjaTwAW5ss1CZ7cUwNYT7NJVYtN80xkTLWbj5xD8GTgdOBT6KYjzfq2qG37LXcNOvJqnqap/lnwIvichJuN7FFeINtJmvqjuAHSJyakX3aRKD3YGYRDUaN+7VpcA+3HhBfyIiF4hIiojsF5GfRGS4iCwSkUV+2zUTkcki8os381uKiIwpb3AiMgA389x9fsnjN6o6T1VTvO1v9Y7b3G8/IiJrROR/3u8dvMdT14rIwyKyCTeUR6PyxmoSlyUQk3B8Jl96U92c0x8Aw0Wksd92J+PuAlKAkcBE3Gxw3fy2a4AbW+oM4G7v54fAZBH5w4Q+ZXCi9zPYu6IXcYMB+k8SdgrQEXjOb/lY3HmMwd2N7S9fmCaR2SMsk4guwX15+q/3+yu4kUvPw827UOQe3NwQZ6k3aJyI/AQsw01QVeRG3Ki3h6pqmrdsnog0Au4SkckaxOxxftp4P9f7LvSG5vf94leoqoWquktE3gTGiMgj+vsgd1cBqaq6yG//W33Py5jysDsQk4hGAWmq+pX3+zxgEz6PsUSkKpAEvOv7IatudkP/+SJOxY3cu1ZEqhW9gDm4ucZ7UXbFta56BjeyatFrvN+6znh3L96Mc2fy57sPgA8seZiKsgRiEoqIHIn7QH9PRBp5dwn1gfeAo0Wk6PFUM9xQ6dsC7Gar3+8tgOP54wd7Hm7qUHBJpKw2eD/b+S2/Hzft7J+mnlXVb4ClwNXeoitwc0QEmmhqczliMuYP7BGWSTSjvZ//8l7+RuHmiSiaK6NFgG1a8sdHSztxiebGYo6ZWo44F3g/hwGPFy1U1fVFxy6mC8hk4DmvCfAVwNuqGmiOb7v7MBVmdyAmYYhIDdysbF/jZmjzf/0AXCIioqoFuG/zI30764nIEbhKaV+zgR64aUKXBnjtKWus3uO1z4GxIuI/d3dJ/oebHO113N3LsyVvbkz52R2ISSTDcI+TbglQqYyIPIf7Bj8IWIjrsDcXeF9EpuAea92Nm7Gw0Kfo47gK+M9F5HHcHUddXFI5TlVHlDPeC3Gz7C0VkWdwCSUbaI5rXQUuWfxGVfeJyMvAP4CfVPXLch7bmFLZHYhJJKNxH7hvF7P+f7g+IaMBVPUT4CKgJ/A+7pHXLbgE8lsHPlXdDQzEdUr8F67y/EVgBC4RlYuqbsRV5D+Cq6h/G5dQJuES4XBVfShA0aLzC1R5bkzI2JS2xpSBiLQB0nEd/O4Nwf4uBV4CugDrytHcN9A+78PVx7RS1awylBPcXOHjgDuB6qGIx1Re9gjLmGJ442U9hmvmuwPoBPwTyAFeCPHh0osOW94diMjhQHdc8phSluThuQiYVt7jm8RjdyDGFMOrdH8TGIB7ZJSNq4e4XVV/DtExmuJTKa+qSyuwrwxcC7E5wCVlrbwXkSa4JFnhWExisARijDGmXKwS3RhjTLlYAjHGGFMulkCMMcaUiyUQY4wx5WIJxBhjTLn8P4gLNtdqbUVcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# you can plot R(t) like this:\n",
    "logz=np.arange(3,-2,-0.01)\n",
    "z=10**logz\n",
    "R=universe.scale_factor(z)\n",
    "t=universe.age(z)\n",
    "plt.plot(t,R)\n",
    "plt.xlabel('Age [Gyr]')\n",
    "plt.ylabel('Scale Factor')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fe462b0-6c2a-4f54-9f1e-d20574b7b2ae",
   "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": 5
}