PSF_GAUSSIAN

[PL]

"""
Simpler version of `PSF_GAUSSIAN` in IDL.

:Authors: Pey Lian Lim (Python)

:Organization: Space Telescope Science Institute

:History:
* 2010/08/17 PLL converted from IDL to Python.
"""

# External modules
import numpy

#-----------
def GaussPsf2D(npix, fwhm, normalize=True):
"""
Parameters
----------
npix: int
Number of pixels for each dimension.
Just one number to make all sizes equal.

fwhm: float
FWHM (pixels) in each dimension.
Single number to make all the same.

normalize: bool, optional
Normalized so total PSF is 1.

Returns
-------
psf: array_like
Gaussian point spread function.
"""

# Initialize PSF params
cntrd = (npix - 1.0) * 0.5
st_dev = 0.5 * fwhm / numpy.sqrt( 2.0 * numpy.log(2) )

# Make PSF
i = range(npix)
psf = numpy.array( [numpy.exp(-(((cntrd-x)/st_dev)**2+((cntrd-y)/
st_dev)**2)/2) for x in i for y in i] )
psf = psf.reshape(npix, npix)

# Normalize
if normalize: psf /= psf.sum()

return psf

[Matt Davis]

We may want to set up a repository for this kind of thing, and a wiki
for describing it, sooner than later.

- Matt

[Rene Breton]

Hey PL,

That's a useful one.

May I suggest you use numpy.indices or something of that kind instead of two for loops, because for loops are awfully slow in python?

x, y = numpy.indices([npix,npix]) - (npix-1)*0.5

psf = numpy.exp( -0.5 * ((x**2 + y**2)/st_dev**2) )

Cheers,

Rene

[Erik Tollerud]

I think this is a great example of the sort of thing that should go in
the IDL affiliated package that's being discussed in the "astropy.idl"
group discussion...

--
Erik Tollerud