import scipy.integrate as integrate
import numpy as np
import matplotlib.pyplot as plt

def integrand(x, a, b):
c = 1. - a - b
return 1/(x*(a*x**(-3)+c*x**(-2)+b)**0.5)


# the five models we wish to test
models = [ (0.,0.), (0.,1.), (1.,0.), (2.5,0.), (0.3,0.7) ]

# calculate t(a) for each model
for (om,ol) in models :
print("Model: om=%.2f, ol=%.2f"%(om,ol))

# we need to calculate t(a) for each a, that we want to
# plot
t_array = np.array([])
a_array = np.linspace(0.001,5.0,100)
ae = 1.e-5
for a1 in a_array :
    t = integrate.quad(integrand,ae, a1, args=(om,ol))
    t_array = np.append(t_array,t)