clear
clc

% resonance order
% Main resonance => Nr=0
Nr=00;

c=3e8;
lambda0=1030e-9;
dlambda0=2e-9;

FSR=33e6;
T1=115e-6;
F=30e3;
CEP=2*pi/10;

L0=c/FSR;
f0=c/lambda0;
df0=f0*dlambda0/lambda0;
n0=round(f0/FSR);
dn0=round(df0/FSR);
n=n0-3*dn0:n0+3*dn0;

fcav=n*FSR;

Pin=sech(1.77*(fcav-f0)/df0).^2;
Ein=sqrt(Pin);

figure(1)
clf
plot(c./fcav*1e9,Pin)
grid on
xlabel('wavelength (nm)')
ylabel('power (A.U.)')
title('normalized laser power spectrum')

df_FSR=linspace(-FSR/2,FSR/2,1e6);
Ec=Ffp(df_FSR,FSR,T1,F,1);

figure(2)
clf
plot(F*df_FSR/FSR,abs(Ec).^2,'r')
grid on
xlim([-5 5])
xlabel('normalized frequency (Hz / LW)')
ylabel('real part (A.U.)')
title('FP-cavity intra-cavity power gain')

dLmax=lambda0/1000;
Nk=1e2;
dL=linspace(-dLmax/2,dLmax/2,Nk+1)-L0*CEP/(2*pi)*FSR/f0;
for k=1:length(dL)
    df=FSR*(dL(k)+Nr*lambda0)/L0;
    Frep=FSR+df;
    flas=(n+CEP/2/pi)*Frep;
    
    [Ec,Er]=Ffp(flas,FSR,T1,F,Ein);
    Gc(k)=sum(abs(Ec).^2)/sum(Pin);
    Gr(k)=sum(abs(Er).^2)/sum(Pin);
end

figure(3)
clf
plot(dL/1e-6,Gc,'.')
grid on
xlabel('length scan (µm)')
ylabel('cavity power gain (A.U.)')

figure(4)
clf
plot(dL/1e-6,Gr,'.')
grid on
xlabel('length scan (µm)')
ylabel('cavity refleced power gain (A.U.)')

% cavity field & reflectivity of the FP cavity
function [Ec,Er] = Ffp(df,FSR,T1,F,Ein)
rho=1-pi/F;
t1=sqrt(T1);
r1=sqrt(1-t1^2);
r2=rho/r1;
if r2>1
    error('r2 > 1 !!!!')
end
Fc=1i*t1./(1-rho*exp(-1i*2*pi*df/FSR));
Ec=Fc.*Ein;
Fr=r1+1i*t1*r2*exp(-1i*2*pi*df/FSR).*Fc;
Er=Fr.*Ein;
end