[AstroPy] Fitting in 3 dimensions

[Daniel Evans]

Dear all,

I have a cube of images derived from an integral field spectrograph, with
axes (x, y, wavelen), with a size of something like (32, 32, 39). In these
images is a star which I would like to use astropy.modelling to fit an
AiryDisk2D to, doing so simultaneously in my N wavelength slices. I'd like
to do it simultaneously in order to better constrain the PSF position and
size using data from all wavelength slices simultaneously - I don't have
any obvious "correct" answer to fix the inputs to.

I envision doing this by fitting N Airy disk models, with x_0 and y_0 fixed
between them, the radius dependent on the wavelength in slice N, and the
amplitude being free to vary. A semi-psudocode version is of what I imagine
the custom model's evaluate() function to be is:

####

def evaluate(x, y, wavelength, x_cen, y_cen, radius_0, amplitude_0...):
    z = numpy.empty_like(x)
    global wavelengthValues

    for N in range(len(wavelengthValues)):
        radius = radius_0 * wavelengthValues[N] / wavelengthValues[0]
        amplitude = amplitude_N

        wavelengthMask = wavelength == N

        z[wavelengthMask] = models.AiryDisk2D.evaluate(
            x[wavelengthMask], y[wavelengthMask], amplitude,
            x_cen, y_cen, radius)

    return z

####

1. This requires fitting in 3D, which I can't see any simple way to do in
astropy - everything seems to support x,y as input and z as output. Have I
missed something, or would I have to really go back to basics?

2. My suggested model is also relatively clunky due to having to define N
amplitude parameters. I glossed over this in the code above, but in reality
the only way I've thought of to do it is to define N input parameters, make
a list from those N input parameters, and then select the correct one from
the list. That's a lot of typing! Is there a more elegant way to do this?

3. Alternatively, can I attack this problem in 2D? Could I easily tile my N
wavelength slices together into a single image, and then called a series of
AiryDisk2D() models that generate data only for the relevant (x, y) of each
tile?

Regards,
Daniel
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/astropy/attachments/20171018/7046c56b/attachment.html>