## Interpolation

Interpolation is a used for many astronomical applications. Interpolation is required to combine sub-pixel dithered images or spectroscopy, sample grids of stellar evolution or stellar atmosphere models, calculate extinction from observed extinction curves, and many many more applications. The scipy.interpolate package in python has some nice built-in interpolation functions and I have gathered a few links describing the capabilities (in addition to the documentation). I would recommend using splrep/splev over interp1d for speed.

- Scipy cookbook on interpolation
- A slideshow outlining the interpolation capabilities in Scipy
- A thread explaining terms in spline interpolation/smoothing

## Integration

Within Scipy, there is an integrate package with several different functions that perform definite or indefinite integrals. I have only played with these briefly; however, I tested both the `scipy.integrate.quad()` and `scipy.integrate.romberg()` functions. These all work in roughly the same way by taking a user-defined function, and the upper and lower boundaries of the integral. Here is a brief example:

from scipy import integrate def myfunc(x, a, b): return (x**b) + a # These are the arguments that will be passed as a and b to myfunc() args = (1.0, -2.0) # Integrate myfunc() from 0.5 to 1.5 results = integrate.quad(myfunc, 0.5, 1.5, args) print 'Integral = ', results[0], ' with error = ', results[1] |

The well-established Fortran QUADPACK integration code lies underneath `quad()` while `romberg()` appears to be a python only implementation (??). While I would much rather use `quad()` for its robustness, it only supports user-functions that accept single value floats; while `romberg()` will pass in vectors of values. For some of my applications, `myfunc()` is computationally intensive, so there is an order of magnitude speed gain by having a properly vectorized `myfunc()` used in vector-form by the integration code. Only with `romberg()` have I matched the speed of IDL’s qpint1d (from Markwardt library, also based on QUADPACK, but vectorized) for identical `myfunc()`.

{ 3 comments… read them below or add one }

You might find that the interpolation package (splrep and splint) is faster and more robust than the integration package (quad or romberg) for doing integration, too! If you have a reasonably smooth function (such that the piece-wise polynomial spline is a good model) you can evaluate it at a reasonably small number of points and do the (very fast) spline integration. This can improve speed for a single integration because the number of function evaluations needed to get an accurate spline model may be small compared to the number of evaluations for quad(), but it is particularly useful if you will be integrating the same curve lots of times (i.e., if only the limits of the integration change) or for surface and volume integrals.

One of our readers ran into an issue with passing in vector quantities into the integration function when using romberg(). You can find the full discussion of the issue and a suggested solution on the scipy-users mailing list:

http://old.nabble.com/Speedup-with-integrate.romberg%28%29-and-vector-arguments-td31485594.html

How do I integrate an interpolated function in Python. I am interpolate Spectrum data and I want to integrate the approximative function from interpolation. Thank you