% This is a self-consistent mode model for the IMC IM and PRC mirrors
% small fuge-factors have been applied to the IM3 RoC and IM3-IM4 distance
% to match the eigen-modes of IMC and PRC
%
% The model was extended to include the SRC on 20171207.
% additional useful reference: E1400205 (not used here)
%
% based on T1200470
addpath(genpath('../lib'));
set(0,'DefaultAxesFontSize',16);
set(0,'DefaultAxesFontWeight','normal');

% fudge factors
% these fudghe factors were good to match RITM=1939
%n=1.44963;
%epsilon3=-0.012483303719256;% change to IM3 RoC to match PRC eigenmode
%epsilon4= 0.031219289189112;% change to IM3 - IM4 length to match PRC eigenmode
% new fudge factors to match mean ITM radius for O3 
%  they are calculated by find_epsilon3_epsilon4.m
n=1.44963;
epsilon3=-0.016544302558746;% change to IM3 RoC to match PRC eigenmode
epsilon4=0.021401970648811;% change to IM3 - IM4 length to match PRC eigenmode

% estimate: guess - need more precise data!
PeriscopeToViewport = 1.6+.82+.29;

%parameters
f1=9099451;
f2=5*f1;
c=299792458;
lambda=1.064e-6;

RPRM=-11;
RPR2=-4.555;
RPR3=36.0;
LPRM=16.6128;         %PRM to PR2
LPR2=16.1626-0.0075;  %PR2 to PR3
LPR3=24.5;            %PR3 to ITM lens

LPRC=3.5*c/2/f1;  % nominal PRC length, given the frequency f1
LSRC=17*c/2/f2;   % nominal SRC length - note that the last 5mm SR2 move did not attempt to preserve that!
LPR3=LPRC-LPR2-LPRM;

RSR3=36.013; % from nebula page, L -  Is that the right value?
RSR2=-6.424;
RSRM=-5.678;

% old reference numbers from OMCpath.m, not used
%loc.ITMY=[-200.0,4983.1,-80.0];
%loc.SR3=[-174.2,-19615.9,-94.5];
%loc.SR2=[-594.1,-4178.1,-104.4];
%loc.SRM=[305.4,-19908.6,-113.2];
%d.ITMYSR3=norm(loc.SR3-loc.ITMY)/1000;
%d.SR3SR2=norm(loc.SR2-loc.SR3)/1000;
%d.SR2SRM=norm(loc.SRM-loc.SR2)/1000;
%d.SRC=d.ITMYSR3+d.SR3SR2+d.SR2SRM;
%%%%%




LARM=3994.5;
LIM4=0.413;  % IM4 to PRM HR
LIM3=1.17;




% X arm
XARM.lambda=lambda;
XARM.q=1;
XARM.N=2;
XARM.dist=[3994.5,3994.5,0];
XARM.ifoc=[2/2244.2,2/1940.3]; % https://galaxy.ligo.caltech.edu/optics/ on 20171208, ETM13, ITM7 (new after replacement in fall 2017)
%XARM.ifoc=[2/2241.54,2/1939.9]; %old: ETM8 ITM3
XARM.label={'ETMX','ITMX'};
XARM=setTelescopeEigenBeam(XARM);

SNPX.lambda=lambda;
SNPX.q=getTelescopeBeam(XARM,'ITMX',1,-1);
SNPX.N=1;
SNPX.dist=[0,0.08]; % 8cm Schnupp Asy, x is longer, from Thesis Evan Hall P1600295
SNPX.ifoc=[-(n-1)/2*XARM.ifoc(2)];
SNPX.label={'ITMXlensr'};

% Y arm
YARM.lambda=lambda;
YARM.q=1;
YARM.N=2;
YARM.dist=[3994.5,3994.5,0];
YARM.ifoc=[2/2246.9,2/1939.2]; % galaxy.ligo.caltech.edu on 20171208, ETM16, ITM11 (new after replacement in fall 2017)
%YARM.ifoc=[2/2238.9,2/1939.2]; %old: ETM12 ITM11
YARM.label={'ETMY','ITMY'};
YARM=setTelescopeEigenBeam(YARM);

SNPY.lambda=lambda;
SNPY.q=getTelescopeBeam(YARM,'ITMY',1,-1);
SNPY.N=1;
SNPY.dist=[0,0.00]; % 8cm Schnupp Asy, x is longer, i.e. y=0, from Thesis Evan Hall P1600295
SNPY.ifoc=[-(n-1)/2*YARM.ifoc(2)];
SNPY.label={'ITMYlensr'};

% construct mean ITM and mean Schnupp
% Note that this is not quite accurate:
% The Michelson x arm sees 85mm more glass from the beamsplitter
% but close enough for now....
qx=getTelescopeBeam(SNPX,'ITMXlensr',1,1);
qy=getTelescopeBeam(SNPY,'ITMYlensr',1,1);
ex=1/qx;
ey=1/qy;
em=(ex+ey)/2; % average the inverse q, since A=1/q*exp(i-kr^2/(2q)) 
qm=1/em;
msnp=(SNPX.dist(2)+SNPY.dist(2))/2; %mean schnupp is 4cm
qm=qm-msnp;
pp=printMode(qm);
RITM=pp.R*n; % calculate the mean ITM radius based on the mean beam curvature







fITM=-RITM/(n-1);
fPRM=-RPRM/(n-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Input Mode Cleaner
IMC.lambda=lambda;
IMC.q=0; %0.2325+13.383212207463499i; % at MC1 HR (is eigenmode of OMC)
IMC.N=4;
IMC.dist=[16.2406,16.2406,0.4650/2,0.4650/2,0];
IMC.ifoc=[2/27.346,0,0,0];
IMC.label={'MC2','MC3','Waist','MC1'};
IMC.Pin=50;
IMC=setTelescopeEigenBeam(IMC); % sets the input beam

%Input Mode Cleaner, from input to output only
IMCp=IMC;
IMCp.N=2;
IMCp.dist=IMC.dist([1,2,5]);
IMCp.ifoc=IMC.ifoc([1,2]);
IMCp.label=IMC.label([1,2]);

% from IMC to PRM
IM.lambda=lambda;
IM.q=getTelescopeBeam(IMCp,'MC3',true,1);
IM.N=6;
IM.dist=[0.4650/2,0.4282+0.0845*n,1.8445-0.4650/2-(0.4282+0.0845*n), 1.1703, 1.5889-LIM4+epsilon4,LIM4,0];
IM.ifoc=[0,0,2/12.8,2/(-6.24+epsilon3),0,1/fPRM];
IM.label={'Waist','IM1','IM2','IM3','IM4','PRMlens'};
IM.Pin=IMCp.Pin;

% PRC forward path
PRC.lambda=lambda;
PRC.q=getTelescopeBeam(IM,'PRMlens',1,1);
PRC.N=3;
PRC.dist=[LPRM,LPR2,LPR3,0];
PRC.ifoc=[2/RPR2,2/RPR3,1/fITM];
PRC.label={'PR2','PR3','ITMlens'};
PRC.Pin=IM.Pin;
% PRC return path
PRCr=PRC;
PRCr.N=5;
PRCr.dist=[0,0,LPR3,LPR2,LPRM,0];
PRCr.ifoc=[-2/RITM,1/fITM,2/RPR3,2/RPR2,2/RPRM];
PRCr.label={'ITMback','ITMlensr','PR3r','PR2r','PRM'};

% PRC round trip
PRCrt=concatTelescope(PRC,PRCr);
PRCrt=setTelescopeEigenBeam(PRCrt);

% input beam from PSL
INP.lambda=lambda;
INP.q=0;
INP.N=6;
INP.dist=[0,PeriscopeToViewport*1000,4018,211.8,140,256.3,167+84.5*n]/1000;
INP.ifoc=[0,0,0,0,0,0];
INP.label={'PSLperiscope','INPUTViewPort','PeriscopeTop','PeriscopeBottom','Steer1','Steer2'};
INP.Pin=IMC.Pin;
INP.q=IMC.q-sum(INP.dist);

%SRC mirror positions from E1400205
BS.x=-203.1;
ITMX.x=5012.6;
ITMY.x=-200.8;
SRM.x=303.8;
SR2.x=-595.4;
SR3.x=-174.2;

BS.y=-183.7;
ITMX.y=-199.8;
ITMY.y=4982.7;
SRM.y=-19874.5;
SR2.y=-4161.1;
SR3.y=-19616.1;

BS.z=-82.5;
ITMX.z=-83.1;
ITMY.z=-79.2;
SRM.z=-113.6;
SR2.z=-104.3;
SR3.z=-94.1;

% relevant optic thicklnesses
BS.thickness=60.41; %mm, from https://galaxy.ligo.caltech.edu/optics/components/BS_06.txt
ITMX.thickness=199.98;  %mm, from https://galaxy.ligo.caltech.edu/optics/components/ITM_07.txt
ITMY.thickness=199.904; %mm, from https://galaxy.ligo.caltech.edu/optics/components/ITM_11.txt
CPX.thickness=99.82;  %mm, from https://galaxy.ligo.caltech.edu/optics/components/CP_01.txt
CPY.thickness=100.24; %mm, from https://galaxy.ligo.caltech.edu/optics/components/CP_05.txt

%SR2 move in fall 2017, 5cm towards BS
SR2.y=SR2.y+5; % in mm

BS.p=[BS.x;BS.y;BS.z]/1e3;
ITMX.p=[ITMX.x;ITMX.y;ITMX.z]/1e3;
ITMY.p=[ITMY.x;ITMY.y;ITMY.z]/1e3;
SRM.p=[SRM.x;SRM.y;SRM.z]/1e3;
SR2.p=[SR2.x;SR2.y;SR2.z]/1e3;
SR3.p=[SR3.x;SR3.y;SR3.z]/1e3;

SRC_dxGlass=(CPX.thickness+ITMX.thickness+BS.thickness*sqrt(2)*2)/1e3; % distance x SRC spends in glass
SRC_dyGlass=(CPY.thickness+ITMY.thickness+BS.thickness*sqrt(2)*1)/1e3; % distance x SRC spends in glass

e14_LSR2=norm(SR2.p-SR3.p);  %SR3 to SR2, from E1400205, before SR2 move
e14_LSRM=norm(SRM.p-SR2.p);  %SR2 to SRM, from E1400205, before SR2 move
e14_L_ITMX_to_SR3_NGC=norm(SR3.p-BS.p)+norm(BS.p-ITMX.p); % no glass correction yet
%e14_L_ITMY_to_SR3_NGC___=norm(SR3.p-BS.p)+norm(BS.p-ITMY.p); % no glass correction yet
e14_L_ITMY_to_SR3_NGC=norm(SR3.p-ITMY.p); % no glass correction yet

e14_L_ITMX_to_SR3_MODE=e14_L_ITMX_to_SR3_NGC -SRC_dxGlass +SRC_dxGlass/n;
e14_L_ITMY_to_SR3_MODE=e14_L_ITMY_to_SR3_NGC -SRC_dyGlass +SRC_dyGlass/n;
e14_L_ITMX_to_SR3_LONG=e14_L_ITMX_to_SR3_NGC -SRC_dxGlass +SRC_dxGlass*n;
e14_L_ITMY_to_SR3_LONG=e14_L_ITMY_to_SR3_NGC -SRC_dyGlass +SRC_dyGlass*n;

e14_SCHNUPP=e14_L_ITMX_to_SR3_LONG-e14_L_ITMY_to_SR3_LONG;
e14_L_ITM_to_SR3_MODE=(e14_L_ITMX_to_SR3_MODE+e14_L_ITMY_to_SR3_MODE)/2;
e14_L_ITM_to_SR3_LONG=(e14_L_ITMX_to_SR3_LONG+e14_L_ITMY_to_SR3_LONG)/2;
e14_LSRC=e14_L_ITM_to_SR3_LONG+e14_LSR2+e14_LSRM;

%
%e12_LSR2=15.4601;  %SR3 to SR2, from E1200616, before SR2 move
%e12_LSRM=15.7400;  %SR2 to SRM, from E1200616, before SR2 move




%SRC single pass
SRC.lambda=lambda;
SRC.q=getTelescopeBeam(PRC,'ITMlens',1,1);
SRC.N=4;
SRC.dist=[0,0,e14_L_ITM_to_SR3_MODE,e14_LSR2,e14_LSRM];
SRC.ifoc=[-2/RITM,1/fITM,2/RSR3,2/RSR2];
SRC.label={'ITMback','ITMlensr','SR3','SR2'};

%SRC return pass
SRCr.lambda=lambda;
SRCr.q=getTelescopeBeam(SRC,'SR2',1,1);
SRCr.N=4;
SRCr.dist=[0,e14_LSRM,e14_LSR2,e14_L_ITM_to_SR3_MODE,0];
SRCr.ifoc=[2/RSRM,2/RSR2,2/RSR3,1/fITM];
SRCr.label={'SRM','SR2r','SR3r','ITMlens'};

SRCrt=concatTelescope(SRC,SRCr);
SRCrt=setTelescopeEigenBeam(SRCrt);



%disp('MC3')
%printMode(getTelescopeBeam(IM,0));
%disp('IM2')
%printMode(getTelescopeBeam(IM,'IM2'));
%disp('IM3')
%printMode(getTelescopeBeam(IM,'IM3'));
%disp('IM4')
%printMode(getTelescopeBeam(IM,'IM4'));
%plotTelescope(concatTelescope(IMC,IMCp,IM))
%plotTelescope(concatTelescope(INP,IMCp,IM,PRC))
%plotTelescope(PRCrt)
%overlap=abs(exactGaussianBeamOverlap(getTelescopeInputBeam(PRC),getTelescopeInputBeam(PRCrt)))^2
figure(1)
plotTelescope(INP);
subplot(211)
title('Input beam: From PSL periscope to MC3')
figure(2)
plotTelescope(IMC);
subplot(211)
title('Input Mode Cleaner')
figure(3)
plotTelescope(IM);
subplot(211)
title('Input Mirros: from IMC to PRM')
figure(4)
plotTelescope(PRC);
subplot(211)
title('Power Recycling Cavity')
figure(5)
plotTelescope(SRC);
subplot(211)
title('Signal Recycling Cavity')
figure(6)
plotTelescope(XARM);
subplot(211)
title('X arm')
figure(7)
plotTelescope(YARM);
subplot(211)
title('Y arm')

display('Mode on ITMX HR, returning from arm:');
printMode(SNPX.q);
display('Mode on ITMY HR, returning from arm:');
printMode(SNPY.q);

qp=getTelescopeBeam(PRCrt,'ITMlensr',1,0);
[amplitudeOverlap,powerOverlap,powerLoss] = gaussianBeamOverlapAmplitude( qm,qp);
display(['Mode-matching of PRC to average arm: ',num2str(powerOverlap*100),' %']);

[amplitudeOverlap,powerOverlap,powerLoss] = gaussianBeamOverlapAmplitude( XARM.q,getTelescopeBeam(PRC,'ITMlens',1,1));
display(['Mode-matching of PRC to X arm: ',num2str(powerOverlap*100),' %']);

[amplitudeOverlap,powerOverlap,powerLoss] = gaussianBeamOverlapAmplitude( YARM.q,getTelescopeBeam(PRC,'ITMlens',1,1));
display(['Mode-matching of PRC to Y arm: ',num2str(powerOverlap*100),' %']);


[amplitudeOverlap,powerOverlap,powerLoss] = gaussianBeamOverlapAmplitude( SRCrt.q,getTelescopeBeam(PRC,'ITMlens',1,1));
display(['Mode-matching of SRC to PRC: ',num2str(powerOverlap*100),' %']);

[amplitudeOverlap,powerOverlap,powerLoss] = gaussianBeamOverlapAmplitude( XARM.q,SRCrt.q);
display(['Mode-matching of SRC to X arm: ',num2str(powerOverlap*100),' %']);

[amplitudeOverlap,powerOverlap,powerLoss] = gaussianBeamOverlapAmplitude( YARM.q,SRCrt.q);
display(['Mode-matching of SRC to Y arm: ',num2str(powerOverlap*100),' %']);

CD=gaussianBeamContrastDefect(qx,qy);
display(['Contrast defect between x and y arm, incl Schnupp asymmetry, but without reflectivity imbalance: ',num2str(CD)]);


