##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
#! /usr/bin/env python
# ConvolveGaussian.py
#
# PURPOSE:
# Generate an image of Gaussians on a white noise background and 
# convolve them with an IPC mask.
#
# CALLING SEQUENCE:
# run ~/python/IPCTest/ConvolveGaussian.py
#
import sys,os
import pyfits
from numpy import *
from pylab import *

import numpy.random as ran
import Gaussian
import convolve

#Get some of the keywords from the keywords file 
if os.path.isfile('/afs/slac/u/ki/lances/python/H4RGAnly/keywords.py'):
  execfile('/afs/slac/u/ki/lances/python/H4RGAnly/keywords.py');

PlotDir = RootDir+'IPC/'

#Generate the original Gaussian
FWHMX  = 10
FWHMY  = 10
Height = 40
FullX  = 51
FullY  = 51
Gauss = Gaussian.gaussian2d(FullX,FullY, Height,FullX/2,FullY/2, FWHMX, FWHMY)
print 'Gauss FWHMX : ' +str(FWHMX)
print 'Gauss FWHMY : ' +str(FWHMY)

#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

#The criteria for an image that contains a series of psf after 
#convultion
MeanBackGround = 2.3
StdDBackGround = .38
PSFReg	       = 51.
NumPSFPerRow   = 20.
ImDim          = PSFReg*NumPSFPerRow
#Generate the Image and add White noise
Im    = zeros((ImDim,ImDim),dtype=float32)
Noise = ran.normal(MeanBackGround,StdDBackGround,size=(ImDim,ImDim))

#Start out with alpha of 0 - No IPC
Alpha = 0.20
for Row in range(NumPSFPerRow):
  Height = Height + 20
  Gauss = Gaussian.gaussian2d(FullX,FullY, Height,FullX/2,FullY/2, FWHMX, FWHMY)
  for Col in range(NumPSFPerRow):
    #Generate the IPC Mask
    Alpha = Alpha-0.0005
    IPCMask[IPCCN] = 1-4*Alpha
    IPCMask[IPCNN] = Alpha
    IPCMask[IPCFN] = 0.0

    #Convolve the Gaussian with this IPC mask
    IPCGauss = convolve.convolve2d(Gauss,IPCMask)
    NoiseIPCGauss = ran.poisson(IPCGauss)
    Im[Col*PSFReg:(Col+1)*PSFReg,Row*PSFReg:(Row+1)*PSFReg] = \
      Noise[Col*PSFReg:(Col+1)*PSFReg,Row*PSFReg:(Row+1)*PSFReg]+IPCGauss[:,:]

#Write the Image to a fits file for examination with pyraf
PSFImName = PlotDir+'IPCTest.fits'
if os.path.exists(PSFImName) : os.remove(PSFImName)
hdu = pyfits.PrimaryHDU(Im)
hdu.writeto(PSFImName)
 
matplotlib.pylab.clf()
matplotlib.pylab.imshow(Im)