In [8]:
%autosave 10
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize, integrate
import astropy.units as u
import astropy.constants as c
from astropy import cosmology
from astroquery.vizier import Vizier
In [7]:
vizier = Vizier(row_limit=1000)
tables = vizier.get_catalogs('J/A+A/568/A22/tablef3')
table = tables[0]
table
Out[7]:
In [27]:
z = table['zcmb']
mb = table['mb']
x1 = table['x1']
color = table['c']
def mu(mb, x1, color, Mb, alpha, beta):
return mb - Mb + alpha*x1 - beta*color
def distance(z, H0, o_l):
"""Luminosity distance
Paramters
---------
z: float
Redshift
H0: float
Hubble constant, km/s/Mpc
o_l: float
Omega_Lambda
Returns
-------
float, pc
"""
if not isinstance(z, np.ndarray):
z = np.array([z])
dist = np.empty(shape=z.shape)
for i, zi in enumerate(z):
d = integrate.quad(
lambda z: 1. / np.sqrt(o_l + (1-o_l)*(1+z)**3),
0, zi
)[0] * (1+zi) * c.c.to_value(u.km/u.s) / H0
dist[i] = d
return 1e6 * dist
def mu_th(z, H0, o_l):
return 5 * np.log10(distance(z, H0, o_l)) - 5
def test_distance():
planck = cosmology.Planck15
z = 1
d1 = distance(
z,
planck.H0.to_value(u.km/u.s/u.Mpc),
planck.Ode0
)
d2 = planck.luminosity_distance(z).to_value(u.pc)
eps = 1e-3
assert abs(d1-d2) < d1*eps
test_distance()
def chi2(H0, mb, x1, color, z):
mu1 = mu(mb, x1, color, -19, 0.14, 3.1)
mu2 = mu_th(z, H0, 0.7)
return np.sum(np.square(mu1-mu2))
optimize.minimize(
chi2,
[50.0],
args=(mb, x1, color, z),
)
Out[27]:
In [ ]: