import scipy as sp
import matplotlib.pyplot as plt
from scipy import optimize, special
import datetime
#
####Read data
#data1 = sp.loadtxt("T0071CH1.txt")
data1 = sp.loadtxt('T0071CH1.csv',delimiter=',', skiprows=16, usecols=(0,1))
x1 = data1.T[0]
y1 = data1.T[1]
#
def exp_decay(parameter, t_):
	S1 = parameter[0]
	N1 = parameter[1]
	t0 = parameter[2]
	decay = sp.select([t_ < t0, t_ >= t0], [N1, N1*sp.exp(-(t_ - t0)*S1)])
	return decay
#
def erf_decay(parameter, t_):
	S2 = parameter[0]
	N2 = parameter[1]
	T0 = parameter[2]
	decay = N2 - N2*special.erf((t_ - T0)*S2)
	return decay
#
def fit_func_exp(parameter, t_, f):
	residual = f - exp_decay(parameter, t_)
	return residual
#
def fit_func_erf(parameter, t_, f):
    residual = f - erf_decay(parameter, t_)
    return residual
#
####Fitting
init_parameter_exp = [0.0, 0.0, 0.0]
result_exp = optimize.leastsq(fit_func_exp,init_parameter_exp,args=(x1,y1))
print(result_exp)
S1_fit=result_exp[0][0]
N1_fit=result_exp[0][1]
t0_fit=result_exp[0][2]
#
init_parameter_erf = [0.0, 0.0, 0.0]
result_erf = optimize.leastsq(fit_func_erf,init_parameter_erf,args=(x1,y1))
print(result_erf)
S2_fit=result_erf[0][0]
N2_fit=result_erf[0][1]
T0_fit=result_erf[0][2]
#
#
####PLOT
plt.rc('text', usetex=True)
time = sp.linspace(-0.005, 0.005, 1000)
fig, ax = plt.subplots(1, 2, figsize=(13, 6))
ax[0].plot(x1, y1, 'o', color='darkgray', markersize=0.1, alpha=0.8, label='raw data')
ax[0].plot(time, exp_decay(result_exp[0], time), 'r-', alpha=0.8, label='exp fit')
ax[0].plot(time, erf_decay(result_erf[0], time), 'b-', alpha=0.8, label='erf fit')
ax[0].set_ylim(0.0, 0.14)
ax[0].set_xlim(-0.005, 0.005)
ax[0].grid(linestyle='--', alpha=0.8)
ax[0].set_xlabel('time [s]')
ax[0].set_ylabel('Intensity [V]')
ax[0].legend(loc='upper right')
ax[0].set_title('Laser intensity decay by cutting beam')
#
plt.rc('text', usetex=True)
plt.rc('font', family='Serif')
s1 = '$ g(t) = N exp (-S(t-T)), (t > T) $'
s2 = '$ f(t) = N - N \int_{-\infty}^t exp \{-S^2(t-T)^2\} dt $'
ax[1].text(0.0, 0.9, s1, ha="left", fontsize=15)
ax[1].text(0.1, 0.8, 'S = {0}'.format(result_exp[0][0]), ha="left")
ax[1].text(0.1, 0.7, 'N = {0}'.format(result_exp[0][1]), ha="left")
ax[1].text(0.1, 0.6, 'T = {0}'.format(result_exp[0][2]), ha="left")
ax[1].text(0.0, 0.4, s2, ha="left", fontsize=15)
ax[1].text(0.1, 0.3, 'S = {0}'.format(result_erf[0][0]), ha="left")
ax[1].text(0.1, 0.2, 'N = {0}'.format(result_erf[0][1]), ha="left")
ax[1].text(0.1, 0.1, 'T = {0}'.format(result_erf[0][2]), ha="left")
ax[1].set_xticks([])
ax[1].set_yticks([])
for side in ['top','bottom', 'left', 'right']:
	ax[1].spines[side].set_visible(False)
#
plt.savefig("{0}_HandDecayFit.jpg".format(datetime.date.today()), dpi=216)
plt.show()