%% compute frequency dependent TF of cavity / double cavity
% based on André Thüring PhD Thesis
% https://pure.mpg.de/pubman/faces/ViewItemOverviewPage.jsp?itemId=item_2320864_1
% section 2.2.1, in particular equation (2.52)

% frequency vector
f = 0:1e-1:1e4;

% cavity parameters
L = 3e3; % [m] arm cavity length
c = 3e8; % [m/s] speed of light
l = 12; % [m] recycling cavity length

% cavity phases
SRdetuned = 40; % [nm] SR cavity detuning
Phi1_0 = SRdetuned*1e-9/1064e-9*2*pi;
Phi2_0 = 0.0;
Phi1 = Phi1_0 + 2*pi*f*l/c; % SR cavity phase
Phi2p = Phi2_0 + 2*pi*f*L/c; % arm cavity phase
Phi2m = Phi2_0 - 2*pi*f*L/c; % arm cavity phase

% mirror reflectivities/transmission
rho2 = (1 - pi/450); % IM
rho3 = sqrt(1 - 4.5e-6); % EM
tau2 = sqrt(1 - rho2^2); % IM
tau1 = sqrt(0.4); % SR
SRCloss = 20e-2;
rho1 = sqrt(1 - 0.4 - SRCloss); % SR

% optial response
rho23p = (rho2 - rho3.*exp(2j*Phi2p))./(1 - rho2*rho3.*exp(2j*Phi2p));
rho23m = (rho2 - rho3.*exp(2j*Phi2m))./(1 - rho2*rho3.*exp(2j*Phi2m));
CF1p = 1j*tau1./(1 - rho1.*rho23p.*exp(2j*Phi1)); % carrier field in SR cavity
CF2p = -1j*tau2./(1 - rho2.*rho3.*exp(2j*Phi2p)).*CF1p; % carrier field in arm cavity
tau3MCp = -1j*tau2.*exp(1j*Phi2p)./(1 - rho2*rho3*exp(2j*Phi2p)).*CF1p.*exp(1j*Phi1);
CF1m = 1j*tau1./(1 - rho1.*rho23m.*exp(2j*Phi1)); % carrier field in SR cavity
CF2m = -1j*tau2./(1 - rho2.*rho3.*exp(2j*Phi2m)).*CF1m; % carrier field in arm cavity
tau3MCm = -1j*tau2.*exp(1j*Phi2m)./(1 - rho2*rho3*exp(2j*Phi2m)).*CF1m.*exp(1j*Phi1);
tau3MCp = conj(tau3MCp);

tau3MC = (tau3MCp + tau3MCm)/2;


% plot response
figure;
subplot(2,1, 1)
loglog(f, abs(tau3MC))
grid;
subplot(2, 1, 2)
semilogx(f, 180/pi*angle(tau3MC));
grid;

title(['SR response, detuned ' num2str(SRdetuned) ' nm, SRC loss ' num2str(SRCloss) ])
