#CalculateAverageChargeLoss.py
#
#PURPOSE:
#  To numerically integrate the charge loss and 
#  obtain the mean value mean value of the 
#  measured charge after IPCT
#

from scipy.integrate import quad
execfile('/afs/slac/u/ki/lances/python/python_HxRG/HxRG_Setup.py')

##PARAMETERS #########################################
Qdep    = 1620.
m       = 0.20
Qlost   = Qdep*m
tau     = 1/30.
t_init  = 0
t_final = 10.6
T       = t_final-t_init
y_init  = Qlost*exp(-(t_init*tau))
y_final = Qlost*exp(-(t_final*tau))

print 'Charge starts at :' + str(y_init)
print 'Charge ends   at :' + str(y_final)

##MEAN
MeanFunc = lambda t, tau, Q, t_final: Q*exp(-t*tau)/t_final

#Numerical
MeanInt = quad(MeanFunc, t_init, t_final, args=(tau, Qlost, t_final))
#Analytical 
MeanAna = -1/(tau*(t_final-t_init))*(y_final-y_init)
#Plot an analytical solution
t        = frange(50000)/100
MeanAnaR =  Qlost*(1/(tau*(t-t_init))*(exp(-t*tau)-exp(-t*t_init))+1)
mplot.figure(0)
mplot.clf()
mplot.plot(t, MeanAnaR)

print 'Average Charge Lost numerical:'
print MeanInt
print 'Average Charge Measured analytical: ' 
print MeanAna
print 'Average Charge Lost anaytical'
print Qlost-MeanAna

##MEAN
StdFunc = lambda y, tau, Q, t_init, t_final, y_init, y_final: \
                (-1/(tau*(t_final-t_init)))*(1/y)*\
                (y-(y_final-y_init)/(-tau*(t_final-t_init)))**2

#Probability Integral
ProbFunc = lambda y, tau, t_init, t_final: \
          1/(-tau*(t_final-t_init))*(1/y)
ProbInt  = quad(ProbFunc, y_init, y_final, args=(tau, t_init, t_final))
print 'Total Probability'
print ProbInt

#Numerical
StdInt = quad(StdFunc, y_init, y_final, \
              args=(tau, Qlost, t_init, t_final, y_init, y_final))
#Analytical
StdAna = -(y_final**2-y_init**2)/(2*tau*T) -\
          (y_final-y_init)**2/(tau**2*T**2)

print 'Variance in Charge Lost numerical:'
print StdInt
print 'Std in Charge Lost numerical:'
print sqrt(StdInt)
print 'Variance in Charge Lost analytical:'
print StdAna
print 'Std in Charge Lost Analytical:'
print sqrt(StdAna)