##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
#! /usr/bin/env python
##
##PURPOSE: To make a Gaussian and convolve it with an IPC mask
import sys,os
import pyfits
from pylab import *
from numpy import *
import Gaussian
import convolve

#Generate the original Gaussian
FWHMX = 20
FWHMY = 20
FullX = 200	#Even Size of Gaussian in X
FullY = 200	#Even Size of Gaussian in Y
Gauss = Gaussian.gaussian2d(FullX,FullY, 400, FullX/2, FullY/2, FWHMX, FWHMY)
print 'Gauss FWHMX : ' +str(FWHMX)
print 'Gauss FWHMY : ' +str(FWHMY)
#Show a contour plot of the gaussian 
#matshow(Gauss, cmap = cm.gist_earth_r)

#IPC Mask
IPCMask = ones((3,3), dtype = float64)
IPCCN   = [1,1]			#Index of Center
IPCNN   = [[1,0,2,1],[0,1,1,2]] #Indices of Nearest neighbors
IPCFN   = [[0,2,0,2],[0,0,2,2]] #Indices of Corner neighbors

#Use the assumption in Moore (p.117) that Center Pixel = 1-4a
Alpha = .20
IPCMask[IPCCN] = 1-4*Alpha
IPCMask[IPCNN] = Alpha
IPCMask[IPCFN] = 0.0

#Convolve the Gaussian with the IPC Mask
IPCGauss = convolve.convolve2d(Gauss,IPCMask)
#Fit it to find the new FWHM
Height, MuX, MuY, FWHMX, FWHMY = Gaussian.moments(IPCGauss)
print 'IPC FWHMX   : ' +str(FWHMX)
print 'IPC FWHMY   : ' +str(FWHMY)

#Difference of the two
Diff = Gauss - IPCGauss
matplotlib.pylab.plot(Diff[FullX/2,:],'g')

#Do the plotting
matplotlib.pylab.clf()
subplot(221)
matplotlib.pylab.plot(IPCGauss[FullX/2,:],'b')
matplotlib.pylab.plot(Gauss[FullX/2,:],'r')
subplot(222)
matplotlib.pylab.plot(Diff[FullX/2,:],'g')

#Show the image
subplot(223)
matplotlib.pylab.imshow(Diff)