This notebook was put together by [Jake Vanderplas](http://www.vanderplas.com) for UW's [Astro 599](http://www.astro.washington.edu/users/vanderplas/Astr599/) course. Source and license info is on [GitHub](https://github.com/jakevdp/2013_fall_ASTR599/).

Introduction to SciPy

Scipy is the premiere package for scientific computing in Python. It has many, many useful tools, which are central to scientific computing tasks you'll come across.

This notebook will walk you through a number of exercises. It purposely does not give you everything you need: it tries to mimic reality, when you'll have to use the built-in documentation + Google searches to figure out how to proceed.

A Quick Survey of SciPy

Here we'll list and discuss some of the submodules available in SciPy:

• File input/output: scipy.io
• Special functions: scipy.special
• Linear algebra operations: scipy.linalg
• Fast Fourier transforms: scipy.fftpack
• Optimization and fit: scipy.optimize
• Statistics and random numbers: scipy.stats
• Interpolation: scipy.interpolate
• Numerical integration: scipy.integrate
• Signal processing: scipy.signal
• Image processing: scipy.ndimage
• Sparse Matrices: scipy.sparse
• Physical Constants: scipy.constants
• Spatial metrics: scipy.spatial

The various exercises below will give you a chance to explore these functions:

scipy.constants

Use values within scipy.constants to compute the Stefan-Boltzmann constant $\sigma$ via the following relation:

Confirm that your expression is correct by comparing the value you calculated to the Stefan-Boltzmann constant available within scipy.constants

scipy.special

Using functions within scipy.special as well as the matplotlib plotting functions we covered yesterday, plot the first six Bessell functions $J_n(x)$ of the first kind, for $0 \le x \le 20$.

Using functions within scipy.special, compute the value of $\log(1000!)$, where log is the natural logarithm, and ! indicates the factorial function (hint: what is the relationship between factorial and the gamma function?).

Bonus: how many digits does $1000!$ have when expressed as a base 10 number?

scipy.linalg

The scipy.linalg module offers a number of interesting routines, some of which overlap with the numpy.linalg module. Both are wrappers of the LAPACK and BLAS tools on your system, so using them will lead to executions as fast as any other languages.

The two linear algebra interfaces are very similar; scipy, thought, adds wrappers of more specialized routines, and often offers more options to tune the results of linear algebra operations.

• Use Scipy's linear algebra tools to compute the value of x, y, and z:

  $$2x + 3y - z = 5$$ $$3x - 2y + 4z = 6$$ $$-x + 2y -z = 12$$

  Hint: remember how a linear system of equations can be turned into a matrix expression? There are several ways of doing this within SciPy.

• Find the inverse of the following matrix:

   [[ 2   3  -1]
    [ 3  -2   4]
    [-1   2  -1]]
  
  

  Find the product of this inverse with the vector

   [5  6  12]

• How many nonzero eigenvalues does the following matrix have?

scipy.interpolate

1. Define an array x which contains 10 points from 0 to $2\pi$

2. Compute an array y which is the sine of the points in x

3. Define an array x_new which contains 1000 points from 0 to $2\pi$

4. Use functions in scipy.interpolate to interpolate x and y onto x_new, using cubic interpolation. Call the results y_new.

5. Use the following commands to plot the results:

    plt.plot(x, y, '.k')
    plt.plot(x_new, y_new, '-k')

scipy.integrate

Numerical integrations is central to many cosmological calculations. Here we'll calculate the angular diameter distance as a function of $z$. For simplicity here, we'll assume a flat universe ($\Omega_k=0$) with contributions only from $\Omega_M$ and $\Omega_\Lambda$.

The equations below refer to Hogg's 1999 paper.

1. computing $E(z)$

Write a function E(z, OmegaM) which computes $E(z)$, as defined by equation 14 of the Hogg paper. The function should take two arguments, z and OmegaM, in that order. Assume ($\Omega_M + \Omega_\Lambda = 1$, $\Omega_k = 0$), and use the default value $\Omega_M = 0.3$.

You should be able to call the function like this:

>>> E(1.0)
1.7606816861659007

Also define a function invE(z, OmegaM) which computes the inverse of E, above.

2. Comoving Distance

Using the tools in scipy.integrate, write a function D_C(z) that computes the comoving distance to $z$ in megaparsecs (eqn. 15 of Hogg). Use $D_H = 3000 Mpc$.

You should be able to call the function like this:

>>> D_C(1)
2314.281199283433

3. Angular Diameter Distance

write a function D_A which computes the angular diameter distance as a function of $z$ (Eqn. 18 in Hogg). Note that because we're assuming a flat cosmology, $D_M = D_C$.

4. Plot the angular diameter distance

Plot the angular diameter distance as a function of $z$, with $z$ ranging from 0 to 4.

Note: your function most likely expects a scalar argument for $z$, so passing an array will not work!

The angular diameter distance is essentially a plot of the ratio of the object's physical size to its apparent (angular) size. You should see a turnover in the plot around $z=1$: what does this mean for the angular size of an object as it gets further away?

scipy.optimize

The optimize submodule contains functions used for function minimization and optimization. There are lots of options in there for multi-dimensional function optimization, which we'll learn about later in the quarter.

• Algebraically determine the value of $x$ which minimizes this function:

  $$f(x) = x^2 + 2x - 4$$

• Numerically compute the value of $x$ which minimizes this function, using a suitable routine from scipy.optimize

• Using the same routine, find the value of $z$ for which $D_A(z)$, above, is maximized. You can reuse your function from the previous section.

Bonus #1

If you're finished and still have time, create a plot of the angular diameter distance turnover redshift as a function of $\Omega_M$, assuming a flat universe with only matter and dark energy ($\Omega_M + \Omega_\Lambda = 1$, $\Omega_k = 0$).

You'll need to combine integration and optimization together, being careful about how you pass-through the additional arguments.

Try to do this without using any loops, i.e. define suitable functions and use the list comprehensions we talked about yesterday.

Bonus #2

If you still have time, repeat the above exercises with a more general function that allows non-flat cosmologies (note that $D_C(z) \ne D_M(z)$ in this case, as indicated in eqn. 16 of the Hogg paper).

Create a contour plot showing the angular diameter distance turnaround redshift as a function of both $\Omega_M$ and $\Omega_\Lambda$.