dove
dove

Reputation: 37

Python plotting an integral of a function

I found a code that calculates:

d_{L}(z) = (1+z) \int_{0}^{z} \frac{cdz}{H(z)}

where

H(z) = H_{0} \sqrt{(1+z)^3\Omega _{m0}+\Omega _{d0}}

But instead of changing \Omega, I want to change z as z = [0.0, 0.1, 0.5, 2.0, 4.0, 5.0]. Here is the code and my modification on z at the end:

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


def H(z, omega_m, H_0=70):
    omega_lambda = 1 - omega_m
    z_prime = ((1 + z) ** 3)
    wurzel = numpy.sqrt(omega_m * z_prime + omega_lambda)

    return H_0 * wurzel


def H_inv(z, omega_m, H_0=70):
    return 1 / (H(z, omega_m, H_0=70))


def integral(z, omega_m, H_0=70):
    I = integrate.quad(H_inv, 0, z, args=(omega_m,))[0]
    return I


def d_L(z, omega_m, H_0=70):
    distance = (2.99 * (10 ** 8)) * (1 + z) * integral(z, omega_m, H_0)
    return distance

z = [0.0, 0.1, 0.5, 2.0, 4.0, 5.0]
omega_m = 0.27


plt.figure()

for i in range(len(z)):
    d_Ls = d_L(z[i], omega_m=omega_m)
    plt.plot(d_Ls, '-', label = 'z = %.1f' % (z[i]))
plt.legend();

But the plot shows nothing. empty_plot

I don't understand what is the problem!

Upvotes: 0

Views: 83

Answers (1)

Stefan
Stefan

Reputation: 957

You are trying to plot single points as a line, which is not possible (as of course you need at least 2 points). You can easily modify that by using scatter instead of plot. You will also need to define the other coordinate of the point, which is z at index i:

for i in range(len(z)):
    d_Ls = d_L(z[i], omega_m=omega_m)
    plt.scatter(z[i], d_Ls, label = 'z = %.1f' % (z[i]))

Upvotes: 1

Related Questions