Firstly you want to construct a master flat field image, which is done in a very similar way to constructing the master bias image:
(rm flat.lis) foreach file ([flat field images].sdf) echo $file:r >> flat.lis end
ls flat*sdf | sed -e 's/\.sdf//' > flat.lisand then
medsky flat.lis masterflat scaled=true
This creates a median-stacked master flat field image called masterflat.sdf. This image has the spatial response of the slit along its rows, the spectral response of the flat-fielding lamp along its columns, and pixel-to-pixel variations due to slight differences in sensitivity for each pixel.
We want to keep the slit response and the pixel-to-pixel variation in the master flat field, but the spectral response of the flat-fielding lamp must be removed. Firstly, collapse the master flat (by adding the counts in each row up) onto the dispersion axis (y-axis) using the starlink package figaro command ystract:
ystract masterflat [xhigh] [xlow] masterflat_spec
where [xhigh] and [xlow] are specified to avoid any useless areas at the sides of each of the spectra, and could be specified as max min if you want.
There are often many more counts towards the red than the blue in flat fields, because the lamps used to illuminate these are much cooler than the stars you have observed, so we will want to work with the logarithm of the spectral response. Take the logarithm of the masterflat_spec.sdf file using the starlink package kappa command log10:
log10 masterflat_spec masterflat_spec_log
Now we have isolated the spectral response of the master flat field, we want to fit it with a polynomial. The pamela command polfit fits an n-order polynomial to one-dimensional data:
polfit masterflat_spec_log masterflat_spec_log_fit [npoly] -3 3 50 TRUE TRUE
where npoly is the order of the polynomial and should be no more than 5 unless you have a good reason. This creates the file masterflat_spec_log_fit.sdf, which is the evaluation of the fitted polynomial at each point in the masterflat_spec_log.sdf.
Important: if there is a lot of structure in the flat field due to the response of a dichroic element in the spectrograph (e.g. the ISIS instrument on the WHT) then do not try to remove it. Leaving it in the flat field means that it will be removed from the science images when they are divided by the master flat.
Masking: if there are dodgy bits at each end of masterflat_spec_log which are causing the fit to be bad, then these areas must be excluded from the fit. We can use the starlink package kappa command ardmask. Create a file mask.ard which will mask out the bits between [pixel-1] to [pixel-2] and [pixel-3] to [pixel-4]. Its syntax should look like:
DIMENSION(1) LINE([pixel-1],[pixel-2]) LINE([pixel-3],[pixel-4])
Once you have a mask.ard file you can apply it:
ardmask masterflat_spec_log mask.ard masterflat_spec_log_mask
and now you use masterflat_spec_log_mask for the fitting with polfit.
Now we have fitted and evaluated a polynomial to the logarithm of the one-dimensional spectral response of the flat-field lamp, we anti-log it and divide the masterflat image by it. The commands are
exp10 masterflat_spec_log_fit masterflat_spec_fit
isydiv masterflat masterflat_spec_fit masterflat_fit
We now have a master flat field image which still contains the effect of the illumination of the slit and the pixel-to-pixel variation (both of which affect the science spectra) but does not contain the spectral response of the flat-field lamp (which does not affect the science spectra). However, its average value will be very small as ystract added pixels together rather than averaging them. Therefore we will use the starlink package figaro command istat to find the mean value of the corrected masterflat:
istat masterflat_fit min max min max
This gives the mean value of the image brightness. Be careful that this mean is not strongly affected by dodgy values in a few bad pixels - if it is then calculate it using a part of the image which is not affected by these. Divide the masterflat by its mean to transform its average value to unity:
icdiv masterflat_fit [mean value] masterflat_fit_mean
Now pamela uses the inverse of the master flat (called the balance) rather than the master flat itself, so we must calculate the inverse of masterflat_fit_mean. Do this by dividing the unit image (the one full of 1's) by the flat field to make the balance image:
idiv unit masterflat_fit_mean balance
Now we have the means by which to flat-field all the science images. The pamela take balance.sdf directly, so we don't need to apply it to the science images ourselves.