clear
clc

c=3e8;
lambda=1030e-9;
Pin=35;
Gcav=20e3;

%% 2M-cavity geometrical setup definition
FSR=216.67e6;           % Free Spectral Range of the FP-cavity
Lrt=c/FSR;              % round trip distance in the FP-cavity
L=Lrt/2;                % distance between mirrors
iR10=0;                 % cold ROC of M1
iR20=1/2.241;           % cold ROC of M2

%% 2M-cavity thermal setup definition
A_Coating=0.6e-6;       % absorption in the coatings
% Heraeus Suprasil 3001 parameters
n_Sup3001=1.45;         % refractive index for Fused Silica
A_Sup3001=0.3e-6;       % 0.3+/-0.2 ppm/cm @ 1064nm
kappa_Sup3001=1.38;     % 1.38 W/m/K @ 20°C    / 1.46W/m/K @ 100°C
alpha_Sup3001=0.6e-6;   % 0.51ppm/K @ 0-100°C / 0.59ppm/K @ 0-300°C
beta_Sup3001=8e-6;      % cf Suprasil 3001 documentation
% Corning 7972 ULE parameters
n_ULE=1.45;             % refractive index for ULE
kappa_ULE=1.31;         % 1.31 W/m/K @ 25°C
alpha_ULE=10e-9;        % Premium grade < 10ppb/K
beta_ULE=11e-6;         % 11.24ppm/K @ 40-60°C / 10.68ppm/K @ 20-40°C

%% 2M-cavity cold mode definition
zw0=iR10*L*(1-iR20*L)/(iR10+iR20-2*L*iR10*iR20);
zr0=sqrt(L*(1-L*iR10)*(1-L*iR20)*(iR10+iR20-L*iR10*iR20))/(iR10+iR20-2*L*iR10*iR20);

% complex radius at z=0 (M1)
q0=-zw0+1i*zr0;
% beam size at z=0 (M1)
wm10=sqrt(lambda/pi/zr0)*abs(q0);
% complex radius at z=L (M2)
qL=L-zw0+1i*zr0;
% beam size at z=L (M2)
wm20=sqrt(lambda/pi/zr0)*abs(qL);

% beam profiler position
Lb=0.67;
% definition of the z-axis
Nz=1e3;
z=linspace(-zr0,L+Lb,Nz);
idm=z<=0;
idp=z>=0 & z<=L;
idb=z>=L;

%% telescope definition
% cold optimization
zwT=zw0;
zrT=zr0;
% hot optimization
zwT=-0.38;
zrT=0.155;
qT=z-zwT+1i*zrT;
wT=sqrt(lambda/pi/zrT)*abs(qT);

%% geometrical coupling definition
C0=4*zrT*zr0/((zwT-zw0)^2+(zrT+zr0)^2);

Nk=200;
Pcav=zeros(1,Nk);
C=C0*ones(1,Nk);
wm1=wm10*ones(1,Nk);
wm2=wm20*ones(1,Nk);
kt=0.05;
iR1th_f=0;
iR2th_f=0;
iR1thl_f=0;
iR2thl_f=0;

figure(1)
clf
hold on
grid on
xlabel('z position (m)')
ylabel('beam size (µm)')
plot(z(idm),wT(idm)*1e6,'r')
ylim([0 900])

for k=1:Nk

    % cavity power
    Pcav(k)=Gcav*C(k)*Pin*(k>1);
    % absorbed power in coatings
    Pa=A_Coating*Pcav(k);
    % thermal ROC for M1 and M2
    iR1th_i=-alpha_Sup3001/(2*pi*kappa_Sup3001*wm1(k)^2)*Pa;
    iR2th_i=-alpha_ULE/(2*pi*kappa_ULE*wm2(k)^2)*Pa;
    % thermal lens for M1 and M2
    iR1thl_i=beta_Sup3001/(2*pi*kappa_Sup3001*wm1(k)^2)*Pa;
    iR2thl_i=-0*beta_ULE/(2*pi*kappa_ULE*wm2(k)^2)*Pa;

    % slow thermal effect simulation
    iR1th_f=iR1th_f+kt*(iR1th_i-iR1th_f);
    iR2th_f=iR2th_f+kt*(iR2th_i-iR2th_f);
    iR1thl_f=iR1thl_f+kt*(iR1thl_i-iR1thl_f);
    iR2thl_f=iR2thl_f+kt*(iR2thl_i-iR2thl_f);

    % total ROC for M1 and M2
    iR1=iR10+iR1th_f;
    iR2=iR20+iR2th_f;
    % total ROC in tranmission for M1 and M2
    iR1t=iR10+iR1thl_f;
    iR2t=iR20+iR2thl_f;

    % cavity mode parameters
    zw=iR1*L*(1-iR2*L)/(iR1+iR2-2*L*iR1*iR2);
    zr=sqrt(L*(1-L*iR1)*(1-L*iR2)*(iR1+iR2-L*iR1*iR2))/(iR1+iR2-2*L*iR1*iR2);
    q0=-zw+1i*zr;
    qL=L-zw+1i*zr;
    q=z-zw+1i*zr;
    w=sqrt(lambda/pi/zr)*abs(q);

    % beam from telescope to cavity
    zwTA=(zwT+2*iR1t*(zrT^2+zwT^2))/(1+4*iR1t*zwT+4*iR1t^2*(zwT^2+zrT^2));
    zrTA=zrT/(1+4*iR1t*zwT+4*iR1t^2*(zwT^2+zrT^2));
    qTA=z-zwTA+1i*zrTA;
    wTA=sqrt(lambda/pi/zrTA)*abs(qTA);

    % beam after the cavity
    zwOUT=(zw+2*iR2t*(zr^2+zw^2))/(1+4*iR2t*zw+4*iR2t^2*(zw^2+zr^2));
    zrOUT=zr/(1+4*iR2t*zw+4*iR2t^2*(zw^2+zr^2));
    qOUT=z-zwOUT+1i*zrOUT;
    wOUT=sqrt(lambda/pi/zrOUT)*abs(qOUT);

    % plots
    plot(z(idp),w(idp)*1e6,'k')
    plot(z(idp),wTA(idp)*1e6,'r')
    %plot(z(idb),wOUT(idb)*1e6,'b')

    % coupling calculation
    if k<Nk
        wm1(k+1)=sqrt(lambda/pi/zr)*abs(q0);
        wm2(k+1)=sqrt(lambda/pi/zr)*abs(qL);
        C(k+1)=4*zrTA*zr/((zwTA-zw)^2+(zrTA+zr)^2);
    end

end

figure(2)
clf
plot(Pcav/1e3)
grid on
ylim([0 max(Pcav/1e3)])
ylabel('cavity power (kW)')

figure(3)
clf
plot(C)
grid on
ylim([0 1])
ylabel('coupling (A.U)')

zwT=(zw-2*iR1t*(zr^2+zw^2))/(1-4*iR1t*zw+4*iR1t^2*(zw^2+zr^2));
zrT=zr/(1-4*iR1t*zw+4*iR1t^2*(zw^2+zr^2));

disp(['waist position for telescope from hot cavity  : ' num2str(zwT) ' m'])
disp(['Rayleigh length for telescope from hot cavity : ' num2str(zrT) ' m'])