clear
clc

global datacut tcut Trt

[h,data]=myisfread('tek0005CH1.isf');
dt=h.Xincr;
Nt=h.Npoints;
t=0:dt:(Nt-1)*dt;
t=t';

% offset removing
offset=mean(data(0.9*Nt:Nt));
data=data-offset;

% magnitude calculation
datamax=mean(data(1:0.1*Nt));

% cuting data to remove constant part at the beginning of the signal
% and where log(data) becomes complex (signal close to zero is sometime negative)
id1=find(data<0.967*datamax);
id2=find(imag(log(data))~=0);
tcut=t(id1(1):id2(1)-1);
datacut=data(id1(1):id2(1)-1);

%fit
L=9;
c=3e8;
Trt=L/c;
rho=1-1.5e-4;
k=1-15e-4;
pc0=0.14;
t0=2e-3;
x0=[rho k pc0 t0];

options = optimset('Display','iter','PlotFcns',@optimplotfval,'TolX',1e-7,'TolFun',1e-7);
x=fminsearch('myerror',x0,options);

rho=x(1);
k=x(2);
pc0=x(3);
t0=x(4);
datafit=pc0/(rho-k)^2*((1-k)*rho.^((tcut-t0)/Trt)-(1-rho)*k.^((tcut-t0)/Trt)).^2;

%%
figure(1)
clf
plot(t*1e3,data)
grid on
xlabel('Time (ms)')
ylabel('Transmited power (A.U)')

figure(2)
clf
plot(tcut*1e3,log(datacut))
hold on
plot(tcut*1e3,log(datafit),'.')
grid on
xlabel('Time (ms)')
ylabel('Log(Transmited power) (A.U)')
legend('data','fit')

Tpin=Trt/2/abs(log(k));
Tcav=Trt/2/abs(log(rho));
F=2*pi*Tcav/Trt;

disp(['Finesse                 : ' num2str(F)])
disp(['Cavity decay time : ' num2str(Tcav*1e6) ' µs'])
disp(['Pin sw decay time : ' num2str(Tpin*1e6) ' µs'])