##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
#Functions that will return gaussians
#
#1D Gaussian
#
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')
from numpy import *
import pdb

#1D Gaussian Generation
#INPUTS:   Height  = Amplitude
#	   XCenter = Mean x
#	   FWHM    = The Full Width at Half Max = 2.3548*Sigma
#
#RETURNS:  Data    = The Gaussian
def gaussian1d(Height,XCenter,FWHM):
    Sigma  = float(FWHM)/2.3548
    return lambda x: Height*exp(-(XCenter-x)**2/(2*Sigma**2))

#1D Gaussian Fit
#INPUTS:   Data    = The Gaussian Data
#          X       = The x values for the data
#
#RETURNS:  Fit	   = The fit to the data
#          Mu	   = The mean
#	   Height  = The amplitude
#	   FWHM    = The Full Width at Half Max
def gaussian1dfit(data, X=0) :

  if len(X) == 1:                   
    X = arange(data.size)                         # Form an array of x values

  #Calculate moments and such
  x = sum(X*data)/sum(data)                       # First Moment
  Var   = abs(sum((X-x)**2*data)/sum(data))       # Second Moment
  Std   = sqrt(Var)                               # Std. Dev
  Max   = data.max()                              # Get the maximum data value
  Gauss  = lambda t : Max*exp(-(t-x)**2/(2*Std**2))
  Fit    = Gauss(X)

  #Return Values
  Height = Max                                    #Height
  Mu     = x                                      #Mean
  FWHM   = 2.3548*Std                             #FWHM

  return Fit, Height, Mu, Std, FWHM

#2D Gaussian
#INPUTS:   Height  = Amplitude
#          XCenter = Mean x
#          FWHM    = The Full Width at Half Max = 2.3548*Sigma
#
#RETURNS:  Data    = The Gaussian
def gaussian2d(XEls, YEls, Height, XCenter, YCenter, FWHMX, FWHMY):
    X, Y    = mgrid[0:XEls,0:YEls]
    SigmaX  = float(FWHMX)/2.3548
    SigmaY  = float(FWHMY)/2.3548
    Gauss2d = lambda x, y: Height*exp(\
	   -(((XCenter-x)/SigmaX)**2+((YCenter-y)/SigmaY)**2)/2)
    Data = Gauss2d(X,Y)
    return Data

#moments
#INPUTS:    Data    = 2-d Gaussian data
#KEYWORDS   SubSize = int, a value to extract a subportion of the array
#RETURNS:   Height  = Amplitude
#	    MuX     = Mean X
#	    MuY     = Mean Y
#           FWHMX   = Full Width at Half Max in X direction
#           FWHMY   = Full Width at Half Max in X direction
def moments(Data, SubSize=1):
 
  if SubSize != 1:
    MidData = (Data.shape)[0]
    Data    = Data[MidData-SubSize:MidData+SubSize, \
                   MidData-SubSize:MidData+SubSize] 

  Total  = Data.sum()
  X, Y   = indices(Data.shape)
  MuData = Data.min()
  MuX    = (X*(Data-MuData)).sum()/Total
  MuY    = (Y*(Data-MuData)).sum()/Total
  col    = Data[:, int(MuY)]-MuData
  YWidth = sqrt(abs((arange(col.size)-MuY)**2*col).sum()/col.sum())
  row    = Data[int(MuX), :]-MuData
  XWidth = sqrt(abs((arange(row.size)-MuX)**2*row).sum()/row.sum())
  Height = Data.max()
  FWHMX  = 2.3548*XWidth
  FWHMY  = 2.3548*YWidth
  FWHM   = (FWHMX+FWHMY)/2

  #Get Ellipticity
  MidX   = X.mean()  ;   MidY = Y.mean()
  X      = X-MidX    ;   Y    = Y-MidY
  Ixx    = sum((Data-MuData)*X**2)/Total
  Ixy    = sum((Data-MuData)*X*Y )/Total
  Iyy    = sum((Data-MuData)*Y**2)/Total
  E1     = (Ixx-Iyy)/(Ixx+Iyy)
  E2     = 2*Ixy/(Ixx+Iyy)
  E      = (E1**2+E2**2)**.5
  EA     = atan2(E2, E1)/2.

  return Height, MuX, MuY, FWHMX, FWHMY, FWHM, E, EA