<jupyter_start><jupyter_text>Week 4 practical - due Friday 9/27 at 11pm* Make sure to put your name in the body of the notebook, and in the filename.* You may discuss these problems with other students in class, but you must turn in your own notebook.* List anyone with whom you discussed the work.* Please document your code, label your plots, and use markdown cells to add explanatory comments to the notebook.* Save your notebook as a pdf file and submit it in Canvas before the deadline.You may not distribute this notebook beyond our class, or post it anywhere online, without permission from the instructor.<jupyter_code>import numpy as np
import sympy as sp
from scipy import integrate
from matplotlib import rc
rc('text', usetex=True)
from sympy import init_printing
init_printing(use_latex=True)
from astropy.cosmology import WMAP9 as cosmo
from astropy import constants as const

import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib
matplotlib.rcParams['mathtext.fontset'] = 'stix'
matplotlib.rcParams['font.family'] = 'serif'
matplotlib.rcParams['font.sans-serif'] = 'stix'
matplotlib.rcParams['text.usetex'] = True
import matplotlib as mpl<jupyter_output><empty_output><jupyter_text>Hubble diagramType Ia supernovae can be used to measure the **Hubble diagram**, which plots the luminosity distance as a function of redshift. Here you will compute the Hubble diagram for different cosmologies.Consider a universe with mass density parameter $\Omega_M$ and dark energy characterized by $\Omega_L$. Let $d_H = c/H_0$ be the Hubble distance, where the Hubble constant is $H_0 = 100 h$ km/s/Mpc. The comoving distance to redshift $z$, in units of $d_H$, is$$  {\hat d}_C(z) = \int_{0}^{z} \frac{dz}{\sqrt{\Omega_M(1+z)^3 + \Omega_\Lambda + \Omega_k (1+z)^2}}$$where$$  \Omega_k = 1 - \Omega_M - \Omega_\Lambda$$The luminosity distance is then$$  d_L(z) = (1+z) d_H \times \begin{cases}   \Omega_k^{-1/2} \sinh\left(\Omega_k^{1/2} {\hat d}_C(z)\right) & \Omega_k > 0 \\   {\hat d}_C(z) & \Omega_k = 0 \\   |\Omega_k|^{-1/2} \sin\left(|\Omega_k|^{1/2} {\hat d}_C(z)\right) & \Omega_k < 0 \end{cases}$$Plot Hubble diagrams out to $z=2$ for the following cases:* flat cosmologies with $(\Omega_M,\Omega_\Lambda) = (0,1)$, $(0.3,0.7)$, $(0.7,0.3)$, and $(1,0)$* matter-only cosmologies with $\Omega_\Lambda = 0$ and $\Omega_M = 0.1$, $0.3$, and $0.5$ Note: `astropy` has functions to compute cosmological distances. In case you are interested, here is sample code that you could modify to test your results. **But you must write your own code here.**<jupyter_code>ztry = np.arange(0.1,1.1,0.1)
Otry = 0.3
Ltry = 0.7

from astropy.cosmology import LambdaCDM
cosmo = LambdaCDM(H0=100.0,Om0=Otry,Ode0=Ltry,Tcmb0=0)
# note that since we set H0=100 in the previous command, we get distances in units of h^{-1} Mpc
cosmo.luminosity_distance(ztry).value
def f1(Ωm, ΩΛ):
    Ωk = 1 - Ωm - ΩΛ
    f = lambda z: (Ωm * (1+z)**3 + ΩΛ + Ωk * (1+z)**2)**(-1/2)
    return f

def d_C(Ωm, ΩΛ):
    ## Defines and evaluates the dc function and returns the value  
    Ωk = 1 - Ωm - ΩΛ
    return integrate.quad(f1(Ωm, ΩΛ), 0 , 2)[0]

def partFunc(Ωm, ΩΛ):
    Ωk = 1 - Ωm - ΩΛ
    display(d_C(Ωm,ΩΛ))
    ### Need to define H with correct units
    d_H = (const.c/cosmo.H(0)).value
    
    if Ωk > 0:
        d_L = lambda z: (1+z) * d_H * (Ωk)**(-1/2) * np.sinh((Ωk**(1/2)) * d_C(Ωm, ΩΛ))
        
    elif Ωk == 0:
        d_L = lambda z: (1+z) * d_H * d_C(Ωm, ΩΛ)
        
    elif Ωk < 0:
        d_L = lambda z: (1+z) * d_H * np.absolute(Ωk**(-1/2))*np.sin(np.absolute(Ωk**(1/2))*d_C(Ωm, Ωk)) 
    
    return d_L

def hPlot(Ωm, ΩΛ):
    z = np.linspace(0, 2, 100)
    display(z)
    f = partFunc(Ωm, ΩΛ)
    display(f)
    plt.clf()
    plt.plot(z, f(z), 'r' )
    plt.show()
        
hPlot(0,1)
#print((const.c/cosmo.H(0)).value)<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup