clear all
close all

P0A_V = [9.34 9.25 9 8.5 7.3 5.2 3.5 1.9 0.7 0.2];
Step0A_V = [0:0.1:0.9];
P0A_H = [9.15 9.1 9 8.75 8.4 7.9 6.6 4.45 2 0.6 0.15];
Step0A_H = [0:0.1:1];

P2A_V = [116 113 111 99 80 53 35 19 6 1.8];
Step2A_V = [0:0.1:0.9];
P2A_H = [113 112.7 112.5 112 110.5 107 101 93.5 78 53 24 7.8 1.9];
Step2A_H = [0:0.1:1.2];

P4A_V = [325 311 288 240 178 116 73 40 13];
Step4A_V = [0:0.1:0.8];
P4A_H = [323 320 318 310 290 254 224 172 104 45 12];
Step4A_H = [0:0.1:1];

subplot(3,1,1)
xlabel('Déplacement des D-shape, o = vertical et * = horizontal')
ylabel('Puissance stockée')
title('Mesures à 0A')
hold on
plot (Step0A_V,P0A_V/9.34,'o-')
plot (Step0A_H,P0A_H/9.15,'*-')
hold off

subplot(3,1,2)
xlabel('Déplacement des D-shape, o = vertical et * = horizontal')
ylabel('Puissance stockée')
title('Mesures à 2A')
hold on
plot (Step2A_V,P2A_V/116,'o-')
plot (Step2A_H,P2A_H/113,'*-')
hold off

subplot(3,1,3)
xlabel('Déplacement des D-shape, o = vertical et * = horizontal')
ylabel('Puissance stockée')
title('Mesures à 4A')
hold on
plot (Step4A_V,P4A_V/325,'o-')
plot (Step4A_H,P4A_H/323,'*-')
hold off

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])

function [TAS,r,r_prod] = GetCavity(cav_name,varargin)
% Return T and r coefficient of a given cavity
% TAS vector contains the 4 T coeffs, then 4 A coeffs, then 4 S coeffs

if strcmp(cav_name,'SBOX_ULE')==1
    TAS(1) = 180e-6;   % T
    TAS(2) = 2e-6;
    TAS(3) = 2e-6;
    TAS(4) = 2e-6;
    TAS(5) = 1.15e-6;  % A
    TAS(6) = 1.27e-6;
    TAS(7) = 1.2e-6;
    TAS(8) = 1e-6;
    TAS(9) =  7e-6;     % S
    TAS(10) = 4.5e-6;
    TAS(11) = 3.6e-6;
    TAS(12) = 9e-6;
%     TAS(1) = 180e-6;   % T
%     TAS(2) = 3.2e-6;
%     TAS(3) = 2.8e-6;
%     TAS(4) = 2.85e-6;
%     TAS(5) = 30e-6;  % A
%     TAS(6) = 30e-6;
%     TAS(7) = 30e-6;
%     TAS(8) = 30e-6;
%     TAS(9) = 20e-6;     % S
%     TAS(10) = 20e-6;
%     TAS(11) = 20e-6;
%     TAS(12) = 20e-6;
    
elseif strcmp(cav_name,'ThomX')==1
    TAS(1) = 120e-6;   % T
    TAS(2) = 1.5e-6;
    TAS(3) = 1.5e-6;
    TAS(4) = 1.5e-6;
    TAS(5) = 0.4e-6;  % A
    TAS(6) = 0.24e-6;
    TAS(7) = 0.24e-6;
    TAS(8) = 0.27e-6;
    TAS(9) = 4e-6;     % S
    TAS(10) = 4.5e-6;
    TAS(11) = 10e-6;
    TAS(12) = 4.5e-6;
    
elseif strcmp(cav_name,'MIGHTY_low')==1
    TAS(1) = 1060e-6;
    TAS(2) = 330e-6;
    TAS(3) = 330e-6;
    TAS(4) = 330e-6;
    TAS(5:12) = 0;

elseif strcmp(cav_name,'Fab_cav')==1
    TAS(1) = 100e-6;
    TAS(2) = 10e-6;
    TAS(3) = 10e-6;
    TAS(4) = 10e-6;
    TAS(5:12) = 0;
end

switch nargin
    case 2
        TAS = repmat(TAS,numel(varargin{1}),1);
        TAS(:,1) = varargin{1};
    case 3
        TAS = repmat(TAS,numel(varargin{1}),1);
        TAS(:,1) = varargin{1};
        TAS(:,2) = varargin{2};
end
    
% Field reflection coeffs
rr = @(TAS) (1-sum(TAS,2)).^(1/2);
for ii = 1:4
    r(:,ii) = rr(TAS(:,ii:4:12));
end
r_prod = prod(r,2);
end