% OMC parameters;           AJW, 1/05

% output mode cleaner:
% PZT-actuated-length, 4-mirror closed path, symmetric.
% T1 and T2 are power transmissivity of the input and output couplers,
% R3 and R4 are the power reflectivity of the internal reflectors.
% R3 is flat, R4 is curved (ROC = 1 meter) and PZT actuated.
% LMC is the HALF length of the full round-trip closed path,
% HALF of the length from from 1->2->3->4->1.
% Since R4 is curved, the waist is between mirrors 2 and 3, closer to 3.

%                            V
%                            |
%                          1 /-\
%                            |  \--------------------\
% 3 |------------------------|------------------------)  4 PZT
%    \--w-----------------\  |
%                          \-/ 2
%                            |
%                            V

clear all

c = 2.9979e8;
Lambda = 1.064e-6;

% output mode cleaner:
% PZT-actuated-length, 4-mirror
LOMC = 0.238;
Tmirr = 1.55e-2;
T1 = Tmirr;
T2 = T1;
% losses to to transmission, scattering, absorption.
L1 = 50e-6;
L2 = L1;
R1 = 1-T1-L1;
R2 = R1;
r1 = sqrt(R1);  r2 = sqrt(R2);
t1 = sqrt(T1);  t2 = sqrt(T2);
% two more HR mirrors, same losses
R3 = 1-L1; R4 = R3;
r3 = sqrt(R3); r4 = r3;
rprod = r1*r2*r3*r4;

% cavity parameters:
FSR = c/(2*LOMC)
fin = pi*sqrt(rprod)/(1-rprod)
BW = FSR/fin;
tstor = fin/(2*pi*FSR)
FWHM = FSR/fin
fpole = FSR*(1-rprod)/(rprod)/(2*pi)
T0 = abs(t1*t2./(1-rprod)).^2
Gain = (t1/(1-rprod))^2

% choose a g-factor:
RCM = 1.0; % radius of curved mirror in meters; assumed R=1.
gfac = (1-2*LOMC/RCM)^2;
ZOMC = LOMC*sqrt(RCM/LOMC-1)
% waist of beam in MC:
wsqr = (Lambda/pi)*ZOMC;
wOMC  = sqrt(wsqr)
zOMC  = pi*wOMC^2/Lambda
div = wOMC/zOMC
% Guoy phase shift
phshift = atan(LOMC/zOMC);
% 1ppm diameter of beam at flat and curved mirrors
dFM1ppm = 5.257 * wOMC*sqrt(1+(0.05/zOMC)^2)
dCM1ppm = 5.257 * wOMC*sqrt(1+(LOMC/zOMC)^2)

% transmission vs frequency shift from resonance
fs = [-1:0.02:1]*4e6;
dL = Lambda/FSR.*fs;
dphi = 2*pi*LOMC/c.*fs;
% cavity transmission
Tcav = abs(t1*t2./(1-rprod.*exp(i*2.*dphi))).^2;
TcavMax = max(Tcav)
% cavity reflectivity
rcav = (r1-r2*(1-L1).*exp(i*2.*dphi))./(1-rprod.*exp(i*2.*dphi));
Rcav = abs(rcav).^2;
% frequency shift of HOMs
dfhom = phshift/pi*FSR;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% summarize all parameters in figure
figure(1)
hold off;
plot(fs,Tcav,'b',fs,Rcav,'g',fs,imag(rcav),'r');
legend('OMC Transmission','OMC Reflection','OMC Imag(ref)',4)
xlabel('\Delta f, Hz')
title('40m Output mode cleaner')
text(fs(10),0.90,['L_{OMC} = ',num2str(LOMC),' mm'])
text(fs(10),0.80,['T_{mirr} = ',num2str(T1*1e6),' ppm'])
text(fs(10),0.70,['L_{back} = ',num2str(L1*1e6),' ppm'])
text(fs(10),0.60,['g_{cav} = ',num2str(gfac)])
text(fs(10),0.50,['waist = ',num2str(wOMC*1e6,3),' \mu m'])
text(fs(10),0.40,['div = ',num2str(div*1e6,3),' \mu rad'])
text(fs(10),0.30,['d(1ppm) CM = ',num2str(dCM1ppm*1e3,3),' mm'])
text(fs(10),0.20,['d(1ppm) FM = ',num2str(dFM1ppm*1e3,3),' mm'])
text(fs(60),0.90,['FSR = ',num2str(FSR*1e-6,3),' MHz'])
text(fs(60),0.80,['finesse = ',num2str(fin,3),''])
text(fs(60),0.70,['FWHM = ',num2str(FWHM*1e-3,4),' kHz'])
text(fs(60),0.60,['f_{pole} = ',num2str(fpole*1e-3,4),' kHz'])
text(fs(60),0.50,['Gain = ',num2str(Gain,3),''])
text(fs(60),0.40,['\tau_{stor} = ',num2str(tstor*1e9,3),' ns'])
text(fs(60),0.30,['\Delta f, HOM = ',num2str(dfhom*1e-6,3),' MHz'])
text(fs(60),0.20,['OMC transmission = ',num2str(TcavMax,3),''])
%print -dpng OMC_params.png

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HOM and RF sideband transmission 

fs = [-0.5:0.0001:0.5]*FSR;
fsM = fs*1e-6;
dL = Lambda/FSR.*fs;
dphi = 2*pi*LOMC/c.*fs;
TRmodeMax = zeros(1,length(fs));
% look at first 8 HOMs
hom = [0:1:8];
for j=1:length(hom),
  phmode = hom(j)*phshift;
  TRm(j) = abs(t1*t2.*exp(i*phmode)./(1-rprod.*exp(i*2.*phmode))).^2;
  TRmode = abs(t1*t2.*exp(i*(dphi+phmode))./(1-rprod.*exp(i*2.*(dphi+phmode)))).^2;
  TRmodeMax = max(TRmodeMax,TRmode);
end
% RF sidebands in MHz:
f11 = 33.207;
f12 = -f11;
f21 = 5*f11;
f22 = 5*f12;
af11 = find(fsM>f11);
Tf11 = TRmodeMax(af11(1));
af12 = find(fsM>f12);
Tf12 = TRmodeMax(af12(1));
af21 = find(fsM>f21);
Tf21 = TRmodeMax(af21(1));
af22 = find(fsM>f22);
Tf22 = TRmodeMax(af22(1));

% plot RF and HOM transmission
figure(2);
plot(fsM,TRmodeMax)
text(fsM(1),0.96,['f11, T11 = ',num2str(f11,3),', ',num2str(Tf11,3)])
text(fsM(1),0.90,['f12, T12 = ',num2str(f12,3),', ',num2str(Tf12,3)])
text(fsM(1),0.84,['f21, T21 = ',num2str(f21,3),', ',num2str(Tf21,3)])
text(fsM(1),0.78,['f22, T22 = ',num2str(f22,3),', ',num2str(Tf22,3)])
axis([fsM(1) fsM(end) 0 1])
ylabel('HOM Transmission')
xlabel('\Delta f, MHz')
title('40m Output mode cleaner')
%print -dpng OMC_RFT.png

% plot vs mode number
figure(3);
semilogy(hom,TRm,'ro')
xlabel('HOM, n+m')
ylabel('Transmission')
title('40m Output mode cleaner')
%print -dpng OMC_HOMT.png

if (0)
figure(33);
hom = [0:1:16];
fhom = dfhom.*hom;
plot(hom,fhom*1e-6,'ro')
xlabel('HOM, n+m')
ylabel('\Delta f, MHz')
title('40m Output mode cleaner')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% explore different g-factors 
% (for fixed RCM, this means different lengths)
frf = 33.207e6;
%frf=0;
gmax = 11;
% set mirror3=1 for 3-mirror cavity, 0 for 4-mirror cavity
mirror3 = 0;

% choose a range of g-factors:
g = 0.001:0.001:1;
mode = 0:gmax;
[Mode,G] = meshgrid(mode,g);

LMC = RCM*(1-G);
Dphi1 = 2*2*pi*LMC/c.*frf;
Dphi2 = 2*2*pi*LMC/c.*5*frf;
Guoy  = 2*Mode.*acos(sqrt(G));
% transmission of carrier and RF sidebands 
T1 = abs(t1*t2./(1-rprod*exp(i*Guoy))).^2;
T2 = abs(t1*t2./(1-rprod*exp(i*Guoy+i*Dphi1))).^2;
T3 = abs(t1*t2./(1-rprod*exp(i*Guoy-i*Dphi1))).^2;
T4 = abs(t1*t2./(1-rprod*exp(i*Guoy+i*Dphi2))).^2;
T5 = abs(t1*t2./(1-rprod*exp(i*Guoy-i*Dphi2))).^2;
if (mirror3)
    T1p = abs(t1*t2./(1-rprod*exp(i*Guoy+i*pi))).^2;
    T2p = abs(t1*t2./(1-rprod*exp(i*Guoy+i*Dphi1+i*pi))).^2;
    T3p = abs(t1*t2./(1-rprod*exp(i*Guoy-i*Dphi1+i*pi))).^2;
    T4p = abs(t1*t2./(1-rprod*exp(i*Guoy+i*Dphi2+i*pi))).^2;
    T5p = abs(t1*t2./(1-rprod*exp(i*Guoy-i*Dphi2+i*pi))).^2;
    T1 = max(T1,T1p);
    T2 = max(T2,T2p);
    T3 = max(T3,T3p);
    T4 = max(T4,T4p);
    T5 = max(T5,T5p);
end

% plot transmission vs half-length, g-factor, or RCM
% X = G;
X = LMC;   % to plot versus half-length for fixed RCM 
% X = 0.238./(1-G);  % to plot versus RCM for fixed half-length
figure(4)
semilogy(X,T1);
if (1)
hold on
semilogy(X,T2);
semilogy(X,T3);
semilogy(X,T4);
semilogy(X,T5);
hold off
end
axis([0.2 0.28 1e-4 1])
%xlabel('g = 1-L_{1/2} / R')
xlabel(' Half-length L_{1/2} (m)')
ylabel('Mode power transmission')
title('40m OMC')
grid on
%print -dpng OMC_gfac.png

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now look at tilt-locking in reflection
TRmode = abs(t1*t2.*exp(i*(dphi+phshift))./(1-rprod.*exp(i*2.*(dphi+phshift)))).^2;
r00 = ((1-L1)*r2*exp(i*2*(dphi))-r1)./(1-rprod.*exp(i*2.*(dphi)));
r01 = ((1-L1)*r2*exp(i*2*(dphi+phshift))-r1)./(1-rprod.*exp(i*2.*(dphi+phshift)));

% the tiltlock signal is given by the difference in power
% on a split photodiode, so long as the tem00 is centered on the PD.
% See, eg, Slagmolen, Shaddock, Gray & McClelland, 
% IEEE J. Quantum Electronics 38, 1521 (2002), eqn 31.
tilt = 0.1;   % radians
power = 1.;   % watts
tiltsig = power*cos(tilt)*sin(tilt)*2*sqrt(2/pi)*real(i*r00.*conj(r01));

figure(5)
plot(dphi,tiltsig,'b')
%plot(dphi,tiltsig,'b',dphi,real(r00),'r',dphi,real(r01),'g',dphi,imag(r00),'r:',dphi,real(r01),'g:')
ylabel('tiltlock signal')
xlabel('\Delta\phi, radians')
title('40m Output mode cleaner')
%print -dpng OMC_tiltsignal.png

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dither lock
% as long as the dither frequency is << FWHM of the cavity,
% it is amplified in the cavity as much as the carrier is.
fdith = 14000;  % Hz
ddith = 0.001;

% now the dither locking in transmission
tcarr = t1*t2./(1-rprod.*exp(i*2.*(dphi)));
tdiff = [tcarr(2:end) tcarr(1)] - [tcarr(end) tcarr(1:(end-1))] ;
tdiff = tdiff/(dphi(3)-dphi(1));
Tdithsig = power*ddith*real(tcarr.*conj(tdiff));
Tran = power*real(tcarr.*conj(tcarr));
% now the dither locking in reflection
rcarr = ((1-L1)*r2.*exp(i*2.*(dphi))-r1)./(1-rprod.*exp(i*2.*(dphi)));
rdiff = [rcarr(2:end) rcarr(1)] - [rcarr(end) rcarr(1:(end-1))] ;
rdiff = rdiff/(dphi(3)-dphi(1));
Rdithsig = power*ddith*real(rcarr.*conj(rdiff));
Refl = power*real(rcarr.*conj(rcarr));

figure(6)
plot(dphi,Tdithsig,'b',dphi,Rdithsig,'r')
ylabel('dither signal')
xlabel('\Delta\phi, radians')
legend('Transmission','Reflection')
title('40m Output mode cleaner')
%print -dpng OMC_dithersignal.png

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% look at dependence of finesse and BW on losses:
LOMC = 0.238;
FSR = c/(2*LOMC);
Loss = 0:0.0001:0.001;
r3 = sqrt(1-Loss);
r4 = sqrt(1-Loss);
Tmirr = 0.0155;
r1 = sqrt(1-Loss-Tmirr);
r2 = sqrt(1-Loss-Tmirr);
fin1 = pi*sqrt(r1.*r2.*r3.*r4)./(1-r1.*r2.*r3.*r4);
T1 = abs(sqrt(Tmirr)*sqrt(Tmirr)./(1-r1.*r2.*r3.*r4)).^2;
Tmirr = 0.0157;
r1 = sqrt(1-Loss-Tmirr);
r2 = sqrt(1-Loss-Tmirr);
fin2 = pi*sqrt(r1.*r2.*r3.*r4)./(1-r1.*r2.*r3.*r4);
T2 = abs(sqrt(Tmirr)*sqrt(Tmirr)./(1-r1.*r2.*r3.*r4)).^2;
Tmirr = 0.01;
r1 = sqrt(1-Loss-Tmirr);
r2 = sqrt(1-Loss-Tmirr);
fin3 = pi*sqrt(r1.*r2.*r3.*r4)./(1-r1.*r2.*r3.*r4);
T3 = abs(sqrt(Tmirr)*sqrt(Tmirr)./(1-r1.*r2.*r3.*r4)).^2;
Tmirr = 0.014;
r1 = sqrt(1-Loss-Tmirr);
r2 = sqrt(1-Loss-Tmirr);
fin4 = pi*sqrt(r1.*r2.*r3.*r4)./(1-r1.*r2.*r3.*r4);
T4 = abs(sqrt(Tmirr)*sqrt(Tmirr)./(1-r1.*r2.*r3.*r4)).^2;

% plot finesse vs loss per mirror
figure(7)
%plot(Loss*100,fin,'r',Loss*100,BW*1e-3,'b')
plot(Loss*100,fin1,'b',Loss*100,fin2,'r',Loss*100,fin3,'g',Loss*100,fin4,'c')
xlabel('Loss per mirror (and there are 4) (%)')
ylabel('Finesse')
title('40m Output mode cleaner')
%legend('T=0.257%','T=0.5%','T=1.0%','T=1.4%')
legend('T=1.55%','T=1.57%','T=1.0%','T=1.4%')
hold on
plot([0.03 0.03],[0 1200],'y')
hold off
grid on
%print -dpng OMC_finesseloss.png

% plot transmission vs loss per mirror
figure(8)
plot(Loss*100,T1,'b',Loss*100,T2,'r',Loss*100,T3,'g',Loss*100,T4,'c')
xlabel('Loss per mirror (and there are 4) (%)')
ylabel('Power transmission')
title('40m Output mode cleaner')
legend('T=1.55%','T=1.57%','T=1.0%','T=1.4%')
hold on
plot([0.03 0.03],[0 1],'y')
hold off
axis([0 .1 0 1])
grid on
%print -dpng OMC_transloss.png

return;