clear
clc

c=3e8;                      % speed of light

% laser parameters
lambda0=1031.6e-9;          % central wavelength (m)
dlambda0=6.2e-9;            % spectral LW (m)
Frep0=160.3e6;              % laser repetition rate (Hz)
CEP0=0;                     % arbitrary CEP value (rad)

% CVBG parameters
CVBG=3;                      % choose the version of the CVBG
switch CVBG
    case 1
        % N40-05
        lambda1=1031.61e-9;   % central wavelength (m)
        dlambda1=2.2e-9;      % spectral LW (m)
    case 2
        % N40-01
        lambda1=1031.55e-9;   % central wavelength (m)
        dlambda1=1.92e-9;     % spectral LW (m)
    case 3
        %N40-20
        lambda1=1031.64e-9;   % central wavelength (m)
        dlambda1=2.49e-9;     % spectral LW (m)
end
lambda_min=lambda1-dlambda1/2; % minimum wavelength limit of the CVBG
lambda_max=lambda1+dlambda1/2; % maximum wavelength limit of the CVBG

% wavelength vector
lambda=linspace(lambda0-5*dlambda0,lambda0+5*dlambda0,1e5);
% laser power vs wavelentgth
Plas=Plaser(lambda,lambda0,dlambda0,0,1);
% laser power after CVBG vs wavelentgth
Pcvbg=Plaser(lambda,lambda0,dlambda0,lambda_min,lambda_max);

figure(1)
clf
plot(lambda*1e9,Plas)
hold on
plot(lambda*1e9,Pcvbg)
grid on
xlabel('wavelength (nm)')
ylabel('laser power (a.u)')
title('laser spectral power before and after CVBG')
legend('before CVBG','after CVBG')

nmin=floor(c/lambda_max/Frep0);   % minimum laser resonance index 
nmax=ceil(c/lambda_min/Frep0);    % maximum laser resonance index 
nmean=(nmin+nmax)/2;              % average laser resonance index 
nv=nmin:nmax;                     % vector of resonance indexes
flas=(nv+CEP0/2/pi)*Frep0;        % vector of laser frequencies 
lambda=c./flas;                   % new vector of wavelength for the laser
% laser power after CVBG vs wavelentgth
Pcvbg=Plaser(lambda,lambda0,dlambda0,lambda_min,lambda_max);

figure(2)
clf
plot(lambda*1e9,Pcvbg)
grid on
xlabel('wavelength (nm)')
ylabel('laser power (a.u)')
title('laser spectral power after CVBG')

% FP-cavity description
FSR=Frep0;          % Free Spectral Range of the FP-cavity
F=23000;            % Finesse of the FP-cavity
LW=FSR/F;           % FP-cavity linewidth definition

N=1e3;              % Nb of CEP simulation steps
cepv=linspace(-4*pi,5*pi,N);    % CEP vector
Gcav=zeros(1,N);                % FP-cavity gain vector initialization
for k=1:N
    dfrep=-cepv(k)/2/pi/(nmean+cepv(k)/2/pi)*FSR; % dfrep = frep - FSR
    df=(nv-nmean).*dfrep;                      % df = flas(n) - n*FSR
    T=Airy(df,LW);                            % power FP-cavity gain vs df
    Gcav(k)=sum(T.*Pcvbg)/sum(Pcvbg);         % total FP-cavity gain
end

figure(3)
clf
semilogy(cepv/pi,Gcav)
grid on
xlabel('CEP/pi (rad/rad)')
ylabel('Relative cavity gain (a.u)')
title('Relative cavity gain vs CEP')
%legend('cvbg N40-05','cvbg N40-01','cvbg N40-20')
hold on

%% measurements
% dfrep = +28Hz (frep@950mA - frep@850mA) => dCEP = -2rad

% measurements on 22/02/2024 Vinput=30mV
Idm=[0 1 2 3];
Vm=[63 164 185 64]/3100;
semilogy((Idm-1.55)*2,Vm,'*') 

% measurements on 23/02/2024 Vinput=30mV
V0=2700;
NIdm=[0 1 2 3 4];
V850m=[45 90 303 66 37]/V0;
V900m=[41 126 202 61 36]/V0;
V950m=[58 164 147 55 34]/V0;
semilogy(NIdm*2-3.5,V850m,'o')
semilogy(NIdm*2-3.25,V900m,'o')
semilogy(NIdm*2-3,V950m,'o')

% measurements on 23/02/2024 Vinput=94mV
V1=650;
Vn850m=663/V1;
Vn900m=114/V1;
Vn950m=74/V1;
semilogy(0,Vn850m,'g+','linewidth',2)
semilogy(0.25,Vn900m,'+','linewidth',2)
semilogy(0.5,Vn950m,'+','linewidth',2)

legend('theoretical curve with F=23000','faom=250MHz I=950mA 22/02/2024', ...
'faom=250MHz I=850mA 23/02/2024','faom=250MHz I=900mA 23/02/2024','faom=250MHz I=950mA 23/02/2024', ...
'faom=210MHz I=850mA 23/02/2024','faom=210MHz I=900mA 23/02/2024','faom=210MHz I=950mA 23/02/2024')

% Laser power after CVBG function
function Pcvbg=Plaser(lambda,lambda0,dlambda0,lambda_min,lambda_max)
Plas=sech(1.7625*(lambda-lambda0)/dlambda0).^2;
Tcvbg=lambda>=lambda_min & lambda<=lambda_max;
Pcvbg=Plas.*Tcvbg;
end

% FP-cavity Airy function
function T=Airy(df,LW)
T=1./(1+(2*df/LW).^2);
end