import common import matplotlib.pyplot as plt import numpy as np import problem11 import problem12 import problem13 import problem14 import setup_matplotlib from problem9 import problem9 if __name__ == "__main__": z_start = 2e7 z_end = 0 N_start = common.N_of_a(common.a_of_z(z_start)) N_end = common.N_of_a(common.a_of_z(z_end)) N = np.linspace(N_start, N_end, 300) # problem 9 ------------------------------------------------------------------------------------ Omega_m_ipl, Omega_r_ipl, Omega_phi_ipl, omega_phi_ipl = problem9(N, 'power law') Omega_m_exp, Omega_r_exp, Omega_phi_exp, omega_phi_exp = problem9(N, 'exponential') # plotting results for problem 9 _, ax = plt.subplots() z = common.z_of_a(common.a_of_N(N)) ax.plot(z, Omega_m_ipl, label=r'$\Omega_m$') ax.plot(z, Omega_r_ipl, label=r'$\Omega_r$') ax.plot(z, Omega_phi_ipl, label=r'$\Omega_\phi$') ax.plot(z, Omega_m_exp, '--', color='C0') ax.plot(z, Omega_r_exp, '--', color='C1') ax.plot(z, Omega_phi_exp, '--', color='C2') # plotting settings ax.set_xscale('log') ax.set_xlabel('Redshift z') ax.set_ylabel(r'Density Parameters $\Omega_i$') ax.xaxis.set_inverted(True) # legend stuff ax.plot(np.nan, np.nan, label='Power Law Potential', color='k') ax.plot(np.nan, np.nan, label='Exponential Potential', color='k', linestyle='--') # get handles and labels handles, labels = ax.get_legend_handles_labels() # specify order of items in legend order = [1, 3, 0, 4, 2] # add legend to plot ax.legend([handles[idx] for idx in order],[labels[idx] for idx in order], loc='upper center', bbox_to_anchor=(0.5, 1.15), ncol=3) plt.savefig('../output/density_parameters.png', bbox_inches='tight', ) _, ax = plt.subplots() ax.plot(z, omega_phi_ipl, label=r'Power Law Potential', color='k') ax.plot(z, omega_phi_exp, '--', label=r'Exponential Potential', color='k') ax.set_xscale('log') ax.set_xlabel('Redshift z') ax.set_ylabel(r'Equation of State $\omega_\phi$') ax.xaxis.set_inverted(True) ax.legend() plt.savefig('../output/equation_of_state.png', bbox_inches='tight') # problem 10 ----------------------------------------------------------------------------------- H_ipl = common.H_over_H0(N[::-1], Omega_m_ipl[-1], 'quintessence', Omega_r_ipl[-1], Omega_phi_ipl[-1], omega_phi_ipl[::-1])[::-1] H_exp = common.H_over_H0(N[::-1], Omega_m_exp[-1], 'quintessence', Omega_r_exp[-1], Omega_phi_exp[-1], omega_phi_exp[::-1])[::-1] Omega_m0_lcdm = 0.3 H_lcdm = common.H_over_H0(N, Omega_m0_lcdm, 'lcdm') fig, ax = plt.subplots() ax.plot(z, H_ipl, label=r'Power Law Potential', color='k') ax.plot(z, H_exp, '--', label=r'Exponential Potential', color='k') ax.plot(z, H_lcdm, '-.',label=r'$\Lambda\mathrm{CDM}$', color='k') ax.set_xlabel(r'$z$') ax.set_ylabel(r'$H(z)/H_0$') ax.set_xscale('log') ax.set_yscale('log') ax.xaxis.set_inverted(True) ax.legend() fig.savefig('../output/hubble_parameter.png', bbox_inches='tight') # problem 11 ----------------------------------------------------------------------------------- t0_ipl = problem11.t0H0(H_ipl, N) t0_exp = problem11.t0H0(H_exp, N) t0_lcdm = problem11.t0H0_lcdm(Omega_m0_lcdm) print(f'Dimensionless age of the universe according to:') print(f' Quintessence model with a power law potential: {(t0_ipl):.3f}') print(f' Quintessence model with an exponential potential: {(t0_exp):.3f}') print(f' \N{GREEK CAPITAL LETTER LAMDA}CDM model: {(t0_lcdm):.3f}') # problem 12 ----------------------------------------------------------------------------------- # Use subset of original solution from z=2 to z=0 mask = z <= 2.0 N_recent = N[mask] z_recent = z[mask] Omega_m_ipl_recent = Omega_m_ipl[mask] Omega_r_ipl_recent = Omega_r_ipl[mask] Omega_phi_ipl_recent = Omega_phi_ipl[mask] omega_phi_ipl_recent = omega_phi_ipl[mask] Omega_m_exp_recent = Omega_m_exp[mask] Omega_r_exp_recent = Omega_r_exp[mask] Omega_phi_exp_recent = Omega_phi_exp[mask] omega_phi_exp_recent = omega_phi_exp[mask] d_L_ipl = problem12.d_L(N_recent, Omega_m_ipl_recent[-1], 'quintessence', Omega_r_ipl_recent[-1], Omega_phi_ipl_recent[-1], omega_phi_ipl_recent) d_L_exp = problem12.d_L(N_recent, Omega_m_exp_recent[-1], 'quintessence', Omega_r_exp_recent[-1], Omega_phi_exp_recent[-1], omega_phi_exp_recent) fig, ax = plt.subplots() ax.plot(z_recent, d_L_ipl, label=r'Power Law Potential', color='k') ax.plot(z_recent, d_L_exp, '--', label=r'Exponential Potential', color='k') ax.set_xlabel(r'$z$') ax.set_ylabel(r'$H_0d_L(z)/c$') ax.xaxis.set_inverted(True) fig.savefig('../output/luminosity_distance.png', bbox_inches='tight') plt.close() # problem 13 ----------------------------------------------------------------------------------- z_obs, d_L_obs, sigma_d_L = problem13.get_data() fig, ax = plt.subplots() ax.errorbar(z_obs, d_L_obs, yerr=sigma_d_L, fmt='o', label='Observations', color='k') ax.plot(z_recent, d_L_ipl / common.h * common.c, label=r'Power Law Potential', color='k') ax.plot(z_recent, d_L_exp / common.h * common.c, '--', label=r'Exponential Potential', color='k') ax.set_xlabel(r'$z$') ax.set_ylabel(r'$d_L(z)$ [Gpc]') ax.xaxis.set_inverted(True) ax.legend() fig.savefig('../output/luminosity_distance_with_data.png', bbox_inches='tight') plt.close() # Convert to Gpc and compute chi-squared chi2_ipl = problem13.chi2(z_recent, d_L_ipl / common.h * common.c, z_obs, d_L_obs, sigma_d_L) chi2_exp = problem13.chi2(z_recent, d_L_exp / common.h * common.c, z_obs, d_L_obs, sigma_d_L) print(f'\nChi-squared values:') print(f' Power law potential: {chi2_ipl:.2f}') print(f' Exponential potential: {chi2_exp:.2f}') # problem 14 ----------------------------------------------------------------------------------- best_Omega_m0 = problem14.best_fit_Omega_m0(N_recent, z_recent, z_obs, d_L_obs, sigma_d_L) print(f'\nBest-fit Omega_m0 for LCDM: {best_Omega_m0[0]:.3f}')