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(-2*pi,2*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
hold on
plot(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')

% 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