gistfile1.py

import numpy as numpy
from astropy.io import fits
from scipy.ndimage.filters import uniform_filter

#Directory: /Users/UCL_Astronomy/Documents/UCL/PHASG199/M33_UVOT_sum/UVOTIMSUM/M33_sum_epoch1_um2_norm.img
with fits.open('...') as ima_norm_um2:
    #Open UVOTIMSUM file once and close it after extracting the relevant values:
    ima_norm_um2_hdr  = ima_norm_um2[0].header
    ima_norm_um2_data = ima_norm_um2[0].data
    #Individual dimensions for number of x pixels and number of y pixels:
    nxpix_um2_ext1 = ima_norm_um2_hdr['NAXIS1']
    nypix_um2_ext1 = ima_norm_um2_hdr['NAXIS2']
    
#Compute the size of the images (you can also do this manually rather than calling these keywords from the header):
#Call the header and data from the UVOTIMSUM file with the relevant keyword extensions:
corrfact_um2_ext1 = numpy.zeros(ima_norm_um2_data.shape)
coincorr_um2_ext1 = numpy.zeros(ima_norm_um2_data.shape)

#Check that the dimensions are all the same:
#print(corrfact_um2_ext1.shape)
#print(coincorr_um2_ext1.shape)
#print(ima_norm_um2_data.shape)

alpha =  0.9842000
ft    =  0.0110329
a1    =  0.0658568
a2    = -0.0907142
a3    =  0.0285951
a4    =  0.0308063

ima_UVM2sum = uniform_filter(ima_norm_um2_data, size=9)
ima_UVM2sum_valid = ima_UVM2sum[4:-4,4:-4]

xvec_UVM2 = ft*ima_UVM2sum_valid
fxvec_UVM2 = 1 + (a1*xvec_UVM2) + (a2*xvec_UVM2**2) + (a3*xvec_UVM2**3) + (a4*xvec_UVM2**4)
Ctheory_UVM2 = - np.log(1-(alpha*ima_UVM2sum_valid*ft)) / (alpha*ft)

corrfact_um2_ext1[4:-4,4:-4] = Ctheory_UVM2*(fxvec_UVM2/ima_UVM2sum_valid)
coincorr_um2_ext1[4:-4,4:-4]  = corrfact_um2_ext1[4:-4,4:-4] * ima_sk_um2_valid

# Make a new image file to save the correction factors:
hdu_corrfact = fits.PrimaryHDU(corrfact_um2_ext1, header=ima_norm_um2_hdr)
fits.HDUList([hdu_corrfact]).writeto('.../M33_sum_epoch1_um2_corrfact.img')

# Make a new image file to save the corrected image to:
hdu_coincorr = fits.PrimaryHDU(coincorr_um2_ext1, header=ima_norm_um2_hdr)
fits.HDUList([hdu_coincorr]).writeto('.../M33_sum_epoch1_um2_coincorr.img')

import numpy as numpy
from astropy.io import fits

#Directory: /Users/UCL_Astronomy/Documents/UCL/PHASG199/M33_UVOT_sum/UVOTIMSUM/M33_sum_epoch1_um2_norm.img
with fits.open('...') as ima_norm_um2:
    #Open UVOTIMSUM file once and close it after extracting the relevant values:
    ima_norm_um2_hdr  = ima_norm_um2[0].header
    ima_norm_um2_data = ima_norm_um2[0].data
    #Individual dimensions for number of x pixels and number of y pixels:
    nxpix_um2_ext1 = ima_norm_um2_hdr['NAXIS1']
    nypix_um2_ext1 = ima_norm_um2_hdr['NAXIS2']
    #Compute the size of the images (you can also do this manually rather than calling these keywords from the header):
    #Call the header and data from the UVOTIMSUM file with the relevant keyword extensions:
    corrfact_um2_ext1 = numpy.zeros((ima_norm_um2_hdr['NAXIS2'], ima_norm_um2_hdr['NAXIS1']))
    coincorr_um2_ext1 = numpy.zeros((ima_norm_um2_hdr['NAXIS2'], ima_norm_um2_hdr['NAXIS1']))

#Check that the dimensions are all the same:
print(corrfact_um2_ext1.shape)
print(coincorr_um2_ext1.shape)
print(ima_norm_um2_data.shape)

# ------------------------------------------------------------------------------------------------------------------------------------------------------------------ #

# Carl: From here on in we just need to work on the for/else loops in order to update the the new correction factors image file and correction image file defined above.

# Define the variables from Poole et al. (2008) "Photometric calibration of the Swift ultraviolet/optical telescope":

alpha =  0.9842000
ft    =  0.0110329
a1    =  0.0658568
a2    = -0.0907142
a3    =  0.0285951
a4    =  0.0308063

ima_UVM2sum = uniform_filter(ima_norm_um2_data, size=9)

ima_UVM2sum_valid = ima_UVM2sum[4:-4,4:-4]

xvec_UVM2 = ft*ima_UVM2sum_valid
fxvec_UVM2 = 1 + (a1*xvec_UVM2) + (a2*xvec_UVM2*xvec_UVM2) + (a3*xvec_UVM2*xvec_UVM2*xvec_UVM2) + (a4*xvec_UVM2*xvec_UVM2*xvec_UVM2*xvec_UVM2)
Ctheory_UVM2 = -(numpy.log(1-(alpha*ima_UVM2sum_valid*ft)))/(alpha*ft)

corrfact_um2_ext1[4:-4,4:-4] = Ctheory_UVM2*(fxvec_UVM2/ima_UVM2sum_valid)
coincorr_um2_ext1[4:-4,4:-4] = corrfact_um2_ext1[4:-4,4:-4]*ima_UVM2sum_valid

# Make a new image file to save the correction factors:
hdu_corrfact = fits.PrimaryHDU(corrfact_um2_ext1, header=ima_norm_um2_hdr)
fits.HDUList([hdu_corrfact]).writeto('.../M33_sum_epoch1_um2_corrfact.img', clobber=True)

# Make a new image file to save the corrected image to:
hdu_coincorr = fits.PrimaryHDU(coincorr_um2_ext1, header=ima_norm_um2_hdr)
fits.HDUList([hdu_coincorr]).writeto('.../M33_sum_epoch1_um2_coincorr.img', clobber=True)