# Python program to load two FITS images,
# take their difference,
# copy the header and add WCS keywords,
# then write out the difference image and header 
# as a new FITS file

# import needed extensions
from numpy import *
#import matplotlib.pyplot as plt # plotting package
#import matplotlib.cm as cm # colormaps
import pyfits


# file names, change as appropriate
skyfile = 'm39-b-10s.fits'
biasfile = 'bias-medavg.fits'
outputfile = 'm39wcs.fits'

# read in the files
sky, header = pyfits.getdata(skyfile, 0, header=True)
bias = pyfits.getdata(biasfile)

# find the difference between the images
reduced = sky-bias

# add keywords to the header
# need to calculate values for CRPIX1/2, CRVAL1/2, CD's
header['WCSAXES'] = 2
header['WCSNAME'] = 'DSS'
header['RADESYS'] = 'ICRS'
header['CTYPE1'] = 'RA---TAN'
header['CUNIT1'] = 'deg' 
header['CRPIX1'] = (1024.0, 'X reference pixel')
header['CRVAL1'] = ( 320.0, 'RA of reference pixel')
header['CTYPE2'] = 'DEC--TAN'
header['CUNIT2'] = 'deg'
header['CRPIX2'] = (1024.0, 'Y reference pixel')
header['CRVAL2'] = (  45.0, 'Dec of reference pixel')
header['CD1_1'] =  -0.0002802103
header['CD1_2'] =  -0.0000006392
header['CD2_1'] =  -0.0000006192
header['CD2_2'] =  0.0002797130

# write out the new FITS file
pyfits.writeto(outputfile, reduced, header, clobber=True)