##^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
##DecayJacobian
##
##PURPOSE:
##  This is the Jacobian used to try and fit the rise/decay function
##  in Persistence_Examine_Exposure
##
##INPUTS
##  x - The readtime array (Read1, Read2, ...., HxRG.NAxis3-1)
##  y - The startime array (Exp1Time, Exp2Time, .... #Exposures)
##  p - A 5 element array
##      p[0] = N1 number of traps
##      p[1] = Time constant Tau1
##      p[2] = N2 number of traps
##      p[3] = Time constant Tau2
##      p[4] = Estimate of sky flux
##OPTIONAL INPUTS:
##  Col - 0 - Do the col_deriv=0
##        1 - Do the col deriv=1, apparently this works alot better
#         
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')
import pdb

def DecayJacobian(p, y, x, z, Col=1):
  #'The Jacobian of the errfunc is -Jacobian of the func'
  N1, Tau1, N2, Tau2, Sk= p

  if Col == 1:
    J = zeros([5, len(x)], dtype=float32 )

    ##dS/dN1
    J[0,:] = -exp(-y/Tau1)*(1-exp(-x/Tau1))
    ##dS/dTau1
    J[1,:] = -N1*(\
               (y/Tau1**2)*exp(-y/Tau1)*(1-exp(-x/Tau1)) + \
               exp(-y/Tau1)*(-x/Tau1**2)*exp(-x/Tau1))
    ##dS/dN2
    J[2,:] = -exp(-y/Tau2)*(1-exp(-x/Tau2))
    ##dS/dTau2
    J[3,:] = -N2*(\
               (y/Tau2**2)*exp(-y/Tau2)*(1-exp(-x/Tau2)) + \
               exp(-y/Tau2)*(-x/Tau2**2)*exp(-x/Tau2))
    ##dS/dSk 
    J[4,:] = -x
    return J
  else:
    J = zeros([len(x),5], dtype=float32 )

    ##dS/dN1
    J[:,0] = -exp(-y/Tau1)*(1-exp(-x/Tau1))
    ##dS/dTau1
    J[:,1] = -N1*(\
               (y/Tau1**2)*exp(-y/Tau1)*(1-exp(-x/Tau1)) + \
               exp(-y/Tau1)*(-x/Tau1**2)*exp(-x/Tau1))
    ##dS/dN2
    J[:,2] = -exp(-y/Tau2)*(1-exp(-x/Tau2))
    ##dS/dTau2
    J[:,3] = -N2*(\
               (y/Tau2**2)*exp(-y/Tau2)*(1-exp(-x/Tau2)) + \
               exp(-y/Tau2)*(-x/Tau2**2)*exp(-x/Tau2))
    ##dS/dSk 
    J[:,4] = -x
    return J


##SimpleExpRiseJacobian
##
##PURPOSE:
##  This is the Jacobian used to try and fit an exponential rise
##
##INPUTS
##  x - The readtime array (Read1, Read2, ...., HxRG.NAxis3-1)
##  y - The values array (ADU1, ADU2, .... #ADUN)
##  p - A 5 element array
##      p[0] = N1 number of traps
##      p[1] = Time constant Tau11
##      p[2] = First ADU value
##OPTIONAL INPUTS:
##  Col - 0 - Do the col_deriv=0
##        1 - Do the col deriv=1, apparently this works alot better
#         
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')
import pdb

def SimpleExpRiseJacobian(p, x, y, Col=1):
  #'The Jacobian of the errfunc is -Jacobian of the func'
  N1, Tau1, Off = p

  if Col == 1:
    J = zeros([3, len(x)], dtype=float32 )

    ##dS/dN1
    J[0,:] = -(1-exp(-x/Tau1))
    ##dS/dTau1
    J[1,:] = N1*(x/(Tau1**2))*exp(-x/Tau1) 
    ##dS/dOffset
    J[2,:] = -1 
    return J
  else:
    J = zeros([len(x),3], dtype=float32 )

    ##dS/dN1
    J[:,0] = -(1-exp(-x/Tau1))
    ##dS/dTau1
    J[:,1] = N1*(x/(Tau1**2))*exp(-x/Tau1)               
    ##dS/dOffset
    J[:,2] = -1         
    return J

##SimpleDecayJacobian
##
##PURPOSE:
##  This is the Jacobian used to try and fit an exponential decay
##      p[0]*exp(-t/p[1])+p[2]
##INPUTS
##  x - The readtime array (Read1, Read2, ...., HxRG.NAxis3-1)
##  y - The values array (ADU1, ADU2, .... #ADUN) or 
##                       (Slope1, Slope2, ...)
##  p - A 3 element array
##      p[0] = Initial val
##      p[1] = Time constant Tau1
##      p[2] = Offset from 0
##OPTIONAL INPUTS:
##  Col - 0 - Do the col_deriv=0
##        1 - Do the col deriv=1, apparently this works alot better
#         
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')
import pdb

def SimpleExpDecayJacobian(p, x, y, Col=1):
  #'The Jacobian of the errfunc is -Jacobian of the func'
  N1, Tau1, Off = p

  if Col == 1:
    J = zeros([3, len(x)], dtype=float32 )

    ##dS/dN1
    J[0,:] = -exp(-x/Tau1)
    ##dS/dTau1
    J[1,:] = -N1*(x/(Tau1**2))*exp(-x/Tau1)
    ##dS/dOffset
    J[2,:] = -1
    return J
  else:
    J = zeros([len(x),3], dtype=float32 )

    ##dS/dN1
    J[:,0] = -exp(-x/Tau1)
    ##dS/dTau1
    J[:,1] = -N1*(x/(Tau1**2))*exp(-x/Tau1)
    ##dS/dOffset
    J[:,2] = -1
    return J

##SimpleDecayJacobian2
##
##PURPOSE:
##  This is the Jacobian used to try and fit an exponential decay
##      p[0]*exp(-t/p[1])+p[2]*exp(-t/p[3])+p[4]
##INPUTS
##  x - The readtime array (Read1, Read2, ...., HxRG.NAxis3-1)
##  y - The values array (ADU1, ADU2, .... #ADUN) or 
##                       (Slope1, Slope2, ...)
##  p - A 5 element array
##      p[0] = Initial val
##      p[1] = Time constant Tau1
##      p[2] = Initial val
##      p[3] = Time constant Tau2 
##      p[4] = Offset from 0
##OPTIONAL INPUTS:
##  Col - 0 - Do the col_deriv=0
##        1 - Do the col deriv=1, apparently this works alot better
#         
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')
import pdb
 
def SimpleExpDecayJacobian2(p, x, y, Col=1):
  #'The Jacobian of the errfunc is -Jacobian of the func'
  N1, Tau1, N2, Tau2, Off = p
  
  if Col == 1:
    J = zeros([5, len(x)], dtype=float32 )
 
    ##dS/dN1
    J[0,:] = -exp(-x/Tau1)
    ##dS/dTau1
    J[1,:] = -N1*(x/(Tau1**2))*exp(-x/Tau1)
    ##dS/dN1
    J[2,:] = -exp(-x/Tau2)
    ##dS/dTau1
    J[3,:] = -N2*(x/(Tau2**2))*exp(-x/Tau2)
    ##dS/dOffset
    J[4,:] = -1
    return J
  else:
    J = zeros([len(x),5], dtype=float32 )
 
    ##dS/dN1
    J[:,0] = -exp(-x/Tau1)
    ##dS/dTau1
    J[:,1] = -N1*(x/(Tau1**2))*exp(-x/Tau1)
    ##dS/dN1
    J[:,2] = -exp(-x/Tau2)
    ##dS/dTau1 
    J[:,3] = -N2*(x/(Tau2**2))*exp(-x/Tau2)
    ##dS/dOffset
    J[:,4] = -1
    return J

##SimpleDecayJacobian3
##
##PURPOSE:
##  This is the Jacobian used to try and fit an exponential decay with
##      p[0]*exp(-t/p[1])+p[2]*exp(-t/p[2])
##INPUTS
##  x - The readtime array (Read1, Read2, ...., HxRG.NAxis3-1)
##  y - The values array (ADU1, ADU2, .... #ADUN) or 
##                       (Slope1, Slope2, ...)
##  p - A 4 element array
##      p[0] = Initial val
##      p[1] = Time constant Tau1
##      p[2] = Initial val
##      p[3] = Time constant Tau2 
##OPTIONAL INPUTS:
##  Col - 0 - Do the col_deriv=0
##        1 - Do the col deriv=1, apparently this works alot better
#         
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')
import pdb

def SimpleExpDecayJacobian3(p, x, y, Off, Col=1):
  #'The Jacobian of the errfunc is -Jacobian of the func'
  N1, Tau1, N2, Tau2 = p

  if Col == 1:
    J = zeros([4, len(x)], dtype=float32 )

    ##dS/dN1
    J[0,:] = -exp(-x/Tau1)
    ##dS/dTau1
    J[1,:] = -N1*(x/(Tau1**2))*exp(-x/Tau1)
    ##dS/dN1
    J[2,:] = -exp(-x/Tau2)
    ##dS/dTau1
    J[3,:] = -N2*(x/(Tau2**2))*exp(-x/Tau2)
    return J
  else:
    J = zeros([len(x),4], dtype=float32 )

    ##dS/dN1
    J[:,0] = -exp(-x/Tau1)
    ##dS/dTau1
    J[:,1] = -N1*(x/(Tau1**2))*exp(-x/Tau1)
    ##dS/dN1
    J[:,2] = -exp(-x/Tau2)
    ##dS/dTau1 
    J[:,3] = -N2*(x/(Tau2**2))*exp(-x/Tau2)
    return J


##DecayJacobian2
##
##PURPOSE:
##  This is the Jacobian used to try and fit the rise/decay function
##  in Persistence_Examine_Exposure
##
##INPUTS
##  x - The readtime array (Read1, Read2, ...., HxRG.NAxis3-1)
##  y - The startime array (Exp1Time, Exp2Time, .... #Exposures)
##  p - A 5 element array
##      p[0] = N1 number of traps
##      p[1] = Time constant Tau11
##      p[2] = Time constant Tau12
##      p[3] = N2 number of traps
##      p[4] = Time constant Tau21
##      p[5] = Time constant Tau22
##      p[6] = Estimate of sky flux
##OPTIONAL INPUTS:
##  Col - 0 - Do the col_deriv=0
##        1 - Do the col deriv=1, apparently this works alot better
#         
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')
import pdb

def DecayJacobian2(p, y, x, z, Col=1):
  #'The Jacobian of the errfunc is -Jacobian of the func'
  N1, Tau11, Tau12,  N2, Tau21, Tau22, Sk= p

  if Col == 1:
    J = zeros([7, len(x)], dtype=float32 )

    ##dS/dN1
    J[0,:] = -exp(-y/Tau11)*(1-exp(-x/Tau12))
    ##dS/dTau11
    J[1,:] = -N1*(y/Tau11**2)*exp(-y/Tau11)*(1-exp(-x/Tau12)) 
    ##dS/dTau11
    J[2,:] = -N1*(x/Tau12**2)*exp(-y/Tau11)*exp(-x/Tau12) 
    ##dS/dN2
    J[3,:] = -exp(-y/Tau21)*(1-exp(-x/Tau22))
    ##dS/dTau11
    J[4,:] = -N1*(y/Tau21**2)*exp(-y/Tau21)*(1-exp(-x/Tau22)) 
    ##dS/dTau11
    J[5,:] = -N1*(x/Tau22**2)*exp(-y/Tau21)*exp(-x/Tau22)    
    ##dS/dSk 
    J[6,:] = -x
    return J
  else:
    J = zeros([len(x),7], dtype=float32 )

    ##dS/dN1
    J[:,0] = -exp(-y/Tau11)*(1-exp(-x/Tau12))
    ##dS/dTau11
    J[:,1] = -N1*(y/Tau11**2)*exp(-y/Tau11)*(1-exp(-x/Tau12)) 
    ##dS/dTau11
    J[:,2] = -N1*(x/Tau12**2)*exp(-y/Tau11)*exp(-x/Tau12)    
    ##dS/dN2
    J[:,3] = -exp(-y/Tau21)*(1-exp(-x/Tau22))
    ##dS/dTau11
    J[:,4] = -N1*(y/Tau21**2)*exp(-y/Tau21)*(1-exp(-x/Tau22))
    ##dS/dTau11
    J[:,5] = -N1*(x/Tau22**2)*exp(-y/Tau21)*exp(-x/Tau22)
    ##dS/dSk 
    J[:,6] = -x 
    return J


##SpatialSpreadJacobian
##
##PURPOSE:
##  This is the Jacobian used to try and fit the spreading function, f,
##  in the form cos(r^(f))
##  in Persistence_Examine_Exposure
##
##INPUTS
##  x - The readtime array (Read1, Read2, ...., HxRG.NAxis3-1)
##  p - A 4 element array
##      p[0]*exp(-(x+p[1])/p[2])+p[3]
##
##OPTIONAL INPUTS:
##  Col - 0 - Do the col_deriv=0
##        1 - Do the col deriv=1, apparently this works alot better
#         
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')
import pdb

def SpatialSpreadJacobian4(p, x, y, Col=0):
  #'The Jacobian of the errfunc is -Jacobian of the func'
  A, TOff, Sig, YOff = p

  if Col == 1:
    J = zeros([4, len(x)], dtype=float32 )

    ##dS/A
    J[0,:] = -exp(-(x+TOff)/Sig)
    ##dS/dTOff
    J[1,:] = -(A*(-1/Sig)*exp(-(x+TOff)/Sig))
    ##dS/dSig
    J[2,:] = -(A*((x+TOff)/(Sig**2))*exp(-(x+TOff)/Sig))
    ##dS/dYOff
    J[3,:] = -1
    pdb.set_trace()
    return J

  else:
    J = zeros([len(x),4], dtype=float32 )

    ##dS/A
    J[:,0] = -exp(-(x+TOff)/Sig)
    ##dS/dTOff
    J[:,1] = -(A*(-1/Sig)*exp(-(x+TOff)/Sig))
    ##dS/dSig
    J[:,2] = -(A*((x+TOff)/(Sig**2))*exp(-(x+TOff)/Sig))
    ##dS/dYOff
    J[:,3] = -1

    return J

def SpatialSpreadJacobian3(p, x, y, Col=0):
  #'The Jacobian of the errfunc is -Jacobian of the func'
  A, Sig, YOff = p

  if Col == 1:
    J = zeros([3, len(x)], dtype=float32 )

    ##dS/A
    J[0,:] = -exp(-(x)/Sig)
    ##dS/dSig
    J[1,:] = -(A*((x)/(Sig**2))*exp(-(x)/Sig))
    ##dS/dYOff
    J[2,:] = -1
    pdb.set_trace()
    return J

  else:
    J = zeros([len(x),3], dtype=float32 )

    ##dS/A
    J[:,0] = -exp(-(x)/Sig)
    ##dS/dSig
    J[:,1] = -(A*((x)/(Sig**2))*exp(-(x)/Sig))
    ##dS/dYOff
    J[:,2] = -1

    return J