clear all; close all;
tic
addpath(genpath('C:\Users\amoudry\desktop\Fichiers Labo\Fichiers Pierre\Simulation\Personal codes\Various'))
[TAS,~,r] = GetCavity('SBOX_ULE');
[F,~] = Get_info(TAS(1:4),TAS(5:8),TAS(9:12));
lambda = 1030e-9;
c = 299792458;
FSR = 133.33e6;
w0 = 2*pi*c/lambda;
tau = 1e-12; % FWHM duration
a = 4*log(2)/tau^2;
E0 = 1;%(pi/2/a)^(1/2); % Energy to normalize gaussian spectrum (Input beam power = 1)
DeltaPhiCE = 0; % CEP
N = 1e5;

dk = 0:(N-1);
Aa = (r.^dk-r.^(2.*N-dk));
Bb = E0*TAS(1)./(1-r.^2);
Cc = (1-r.^(2.*N));

Nn = 5e2; % Increase Nn <-> increase resolution
dtt = -Nn:Nn;
dtt = dtt*lambda/c/(0.1*Nn); %1e6
Ecn = zeros(numel(dtt),1);
for ii = 1:numel(dtt)
    for ll = 0:3
%         ll = 0;
        dt = dtt(ii)+ll*lambda/c;
        Phid = DeltaPhiCE + w0.*dt;
        temp_vect = Aa.*cos(dk.*Phid).*exp(-a.*dk.^2.*dt.^2./2);
        Ecn(ii,ll+1) = Bb.*(2*sum(temp_vect)-Cc);
        disp([num2str(ii)]);
    end
end
toc
%% Time plots
% figure
% semilogy(dtt/lambda*c,Ecn/max(Ecn),'LineWidth',2)
% hold on
% semilogy(dtt/lambda*c,Ecn2/max(Ecn),'LineWidth',2)
% hold on
% semilogy(dtt/lambda*c,Ecn3/max(Ecn),'LineWidth',2)
% hold on
% semilogy(dtt/lambda*c,Ecn4/max(Ecn),'LineWidth',2)
% set(gca,'FontSize',15)
% xlabel('\DeltaT (\lambda/c)')
% ylabel('Energy (A.U.)')
% grid on
% legend('\DeltaT = 0','\DeltaT = \lambda/c','\DeltaT = 2\lambda/c','\DeltaT = 3\lambda/c')
% axis square

% figure
% semilogy(dtt/lambda*c,Ecn/max(Ecn),'LineWidth',2)
% grid on
% set(gca,'FontSize',25)
% % set(gca,'YLim',[1e-9 1e0])
% xlabel('\DeltaT (\lambda_0/c)')
% ylabel('log(Energie (u.a.))')

%% Frequency
nu0 = w0/2/pi;
frep = (1/FSR-nu0/FSR*dtt).^(-1); % Infinity in dtt = 1/nu0

fprintf('\nFinesse : %g\n\n',F);
% Get linewidth
figure
for jj = 1:4
%     Find the 2 minimas of Ecn_half. Take the corresponding frep and
%     substract them
    Ecn_half = abs(Ecn(:,jj)-max(Ecn(:,jj))/2);
    Ecn_half2 = sort(Ecn_half);
    [row1,~] = find(Ecn_half==Ecn_half2(1),1);
    [row2,~] = find(Ecn_half==Ecn_half2(2),2);
    if numel(row2)>1  % Sometimes row can be a vector
        row2 = row2(2);
    end
    dnu = abs(frep(row2)-frep(row1));
    fprintf('RES %g\nMax gain : %g. Linewidth : %g kHz\n\n',jj-1,max(Ecn(:,jj)),dnu/1e3);
    
%     plot((frep-FSR)/FSR,Ecn(:,jj)/max(Ecn(:,1)),'LineWidth',2)
%     hold on
    semilogy((frep-FSR)/FSR,Ecn(:,jj)/max(Ecn(:,1)),'LineWidth',2)
    xlim([-0.05 0.05])
    hold on
end
set(gca,'FontSize',15)
xlabel('(f_r_e_p-FSR)/FSR')
% ylabel('Energy (A.U.)')
ylabel('log(Energie (u.a.))')
grid on
legend('\DeltaT = 0','\DeltaT = \lambda/c','\DeltaT = 2\lambda/c','\DeltaT = 3\lambda/c')
axis square
% axis([-0.01 0.01 10^-6 1])