CODEWITHSUNDEEP
pythonscipyquad
user7248647
user7248647

Reputation: 67

Python: Plotting an Integral

I know there are some similar questions around here but none of them seems to really get to my problem.

My code looks like this:

import numpy 

import matplotlib.pyplot as plt

from scipy import integrate as 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,))
    return I


def d_L(z, omega_m , H_0=70):

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

    return distance

The functions do work, my problem: How can I plot d_L versus z ? Like it's obviously a problem that I have this integral function in the definition of my d_L and it depends on z and some args=(omega_m, ).

Upvotes: 3

Views: 9569

Answers (2)

Julien Marrec
Julien Marrec

Reputation: 11895

To build upon @eyllanesc's solution, here's how you'd plot several values of omega:

import numpy

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

z0 = -1.8
zf = 10
zs = numpy.linspace(z0, zf, 1000)

fig, ax = plt.subplots(nrows=1,ncols=1, figsize=(16,9))

for omega_m in np.linspace(0, 1, 10):
    d_Ls = numpy.linspace(z0, zf, 1000)
    for index in range(zs.size):
        d_Ls[index] = d_L(zs[index], omega_m=omega_m)
    ax.plot(zs,d_Ls, label='$\Omega$ = {:.2f}'.format(omega_m))
ax.legend(loc='best')
plt.show()

Plot with multiple lines

Upvotes: 2

eyllanesc
eyllanesc

Reputation: 243907

try with my solution:

import numpy

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

z0 = -1.8
zf = 10
zs = numpy.linspace(z0, zf, 1000)
d_Ls = numpy.linspace(z0, zf, 1000)
omega_m = 0.2

for index in range(zs.size):
     d_Ls[index] = d_L(zs[index], omega_m=omega_m)
plt.plot(zs,d_Ls)
plt.show()

Output:

enter image description here

Upvotes: 1

Related Questions