# Extracting PSF Images

## Generic Code to Reconstruct the PSF

The Point-Spread Function (PSF) reconstruction can be performed without any specialized tools.

To read the PSF info from a psField file for the r band, using IDL, read extension 3:

IDL> pstruct = mrdfits(psfield_file, 3)

(recall u,g,r,i,z == 0,1,2,3,4 so just add one to the index).

The resulting structure can then be used to reconstruct the image. The following IDL code, taken from the SDSSIDL routine "sdss_psfrec.pro", demonstrates the algorithm to reconstruct the PSF at location (row,col) in the field.

nrow_b=(pstruct.nrow_b)
ncol_b=(pstruct.ncol_b)
;assumes they are the same for each eigen so only use the 0 one
rnrow=(pstruct.rnrow)
rncol=(pstruct.rncol)
nb=nrow_b*ncol_b
coeffs=fltarr(nb)
ecoeff=fltarr(3)
cmat=pstruct.c
rcs=0.001
FOR i=0L, nb-1L DO coeffs[i]=(row*rcs)^(i mod nrow_b) * (col*rcs)^(i/nrow_b)
FOR j=0,2 DO BEGIN
FOR i=0L, nb-1L DO BEGIN
ecoeff[j]=ecoeff[j]+cmat(i/nrow_b,i mod nrow_b,j)*coeffs[i]
ENDFOR
ENDFOR
psf = (pstruct.rrows)[*,0]*ecoeff+$(pstruct.rrows)[*,1]*ecoeff+$
(pstruct.rrows)[*,2]*ecoeff


The above code can be easily rewritten in other programming languages. For example, the following function allows to reconstruct the PSF in Python.

from astropy.io import fits
import numpy as np

def reconstructPSF(psFieldFilename, filter, row, col):

filterIdx = 'ugriz'.index(filter) + 1
psField = fits.open(psFieldFilename)
pStruct = psField[filterIdx].data

nrow_b = pStruct['nrow_b']
ncol_b = pStruct['ncol_b']

rnrow = pStruct['rnrow']
rncol = pStruct['rncol']

nb = nrow_b * ncol_b
coeffs = np.zeros(nb.size, float)
ecoeff = np.zeros(3, float)
cmat = pStruct['c']

rcs = 0.001
for ii in range(0, nb.size):
coeffs[ii] = (row * rcs)**(ii % nrow_b) * (col * rcs)**(ii / nrow_b)

for jj in range(0, 3):
for ii in range(0, nb.size):
ecoeff[jj] = ecoeff[jj] + cmat[int(ii / nrow_b), ii % nrow_b, jj] * coeffs[ii]

psf = pStruct['rrows'] * ecoeff + \
pStruct['rrows'] * ecoeff + \
pStruct['rrows'] * ecoeff

psf = np.reshape(psf, (rnrow, rncol))
# psf = psf[10:40, 10:40]  # Trim non-zero regions.

return psf


## Stand Alone Code

We have also created a standalone code that serves to both read them and as a template library for inclusion in other codes. The code is available as: readAtlasImages-v5_4_11.tar.gz.

### Compiling

1. make clean
2. make

If you are on a big-endian machine, remove -DSDSS_LITTLE_ENDIAN from CFLAGS in the Makefile.

### Using

% read_PSF -h
Usage: read_PSF [options] input-file hdu output-file
-?      This message
-h      This message
-i      Print an ID string and exit
-v      Turn up verbosity (repeat flag for more chatter)


To reconstruct the z PSF (i.e., the 5th HDU) at the position (row, col) = (500, 600) from run 1336, column 2, field 51 you'd say:

%  read_PSF psField-001336-2-0051.fit 5 500.0 600.0 foo.fit

The desired PSF would appear as an unsigned short FITS file in foo.fit; the background level is set to the standard 'soft bias' of 1000. If you want a floating image, change a line in the read_PSF.c; look for

        /* create a float region */