% H1 AS beam modeling script

% This script takes the input beam from the MC waist and propagates using
% the (hopefully) realistic as-built parameters of the PRC, ITMs, SRC and
% predicts various beam parameters at the AS port and in HAM6.

% It also takes measured beam properties from ISCT6 and propagates the beam 
% upstream to the SRM and downstream to the OMC, and measures the overlap at
% the OMC and the Guoy phase and spot size at the WFS

% Thing to do - vary assumptions of optics like OMs and ISCT6 lenses and see how
% sensitive measurements on ISCT6 will be to errors in these parameters

%% 
clear all
close all
%% Initialization

lambda = 1.064e-6;

% Get as-built optical parameters from H1parameters file
H1parametersNEW_ITMs;

% Build the beam from the IMC to the BS
% Lots of this is taken from Lisas script in LHO:13815

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

%% Cross check of PRC, SRC and Schnupp Asy lengths

% Optical length PR3 --> ITMX AR
dPR3_ITXAR = dPR3_BSHR + dBSHR_BSAR * 1.45 + dBSAR_CPX1 + dCPX1_CPX2 *1.45 + dCPX2_ITXAR;
dPRM_ITXAR = dPRM_PR2 + dPR2_PR3 + dPR3_ITXAR;
dPRM_BSHR = dPRM_PR2 + dPR2_PR3 + dPR3_BSHR;

% Optical length BS --> ITMY HR / ITX HR
dBSHR_ITYHR = dBSHR_CPY1 + dCPY1_CPY2 * 1.45 + dCPY2_ITYAR + dTYAR_ITYHR * 1.45;
dBSHR_ITXHR =  dBSHR_BSAR * 1.45 + dBSAR_CPX1 + dCPX1_CPX2 * 1.45 + dCPX2_ITXAR + dITXAR_ITXHR * 1.45;

% Optical length SR3 --> ITMY AR
dSR3_BSHR = dSR3_BSAR + dBSAR_BSHR * 1.45;
dSR3_ITYAR = dSR3_BSHR + dBSHR_CPY1 + dCPY1_CPY2 *1.45 + dCPY2_ITYAR;


% Cross check optical length of signal recycling cavities
PRC = dPRM_PR2 + dPR2_PR3 + dPR3_BSHR + (dBSHR_ITXHR + dBSHR_ITYHR) / 2;
SRC = dSRM_SR2 + dSR2_SR3 + dSR3_BSHR + (dBSHR_ITXHR + dBSHR_ITYHR) / 2;

Schnupp = dBSHR_ITXHR - dBSHR_ITYHR;

% Build beam path to AS port
asX = beamPath;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Beam propagation from the IMC, based on Pauls measured IO parameters

% Input waist in IMC is in between MC1 and MC3

% Seed for X path
asX.seedq = beamq(-0.2325 + 13.3377i); % q at MC3 HR - path for single bounce X
asX.seedz = 0;
% Seed for Y path - do not think this is necessary?
%asY.seedq = beamq(-0.2325 + 13.3377i); % q at MC3 HR - path for single bounce Y
%asY.seedz = 0;


% From IMC to PRM
% MC3 substrate path length = 84.5mm; 122.6 MC3-AR surface to IM1 = 428.2mm
dWaist_SM1 = (122.6 + 428.2)*1e-3;
asX.addComponent(component.flatMirror(dWaist_SM1, 'SM1f'));
% IM1->IM2 = 1293.8mm, IM2 Rc = 12800mm, IM2 AOI = 7deg
asX.addComponent(component.curvedMirror(ROC_PMMT1*cos(7*pi/180), asX.component('SM1f').z + dSM1_PMMT1_eff, 'PMMT1f'));
%IM2->IM3 = 1170.4mm, IM3 Rc = -6240mm
asX.addComponent(component.curvedMirror(ROC_PMMT2*cos(7*pi/180), asX.component('PMMT1f').z + dPMMT2_PMMT1_eff, 'PMMT2f'));
%IM3 AOI = 7.1deg, IM3->IM4 = 1174.5mm
asX.addComponent(component.flatMirror(asX.component('PMMT2f').z + dSM2_PMMT2_eff, 'SM2f'));
% PRM lens, IM4->PRM-AR surface = 413.5mm PRM substrate path length = 73.7mm
asX.addComponent(component.dielectric(inf, R_PRM, tPRM, nCPX, asX.component('SM2f').z + dPRM_SM2_eff + tPRM/idxPRM, 'PRMHR'));

% Path from PRM to ITX HR
asX.addComponent(component.curvedMirror(R_PR2,  asX.component('PRMHR').z + dPRM_PR2, 'PR2'));
asX.addComponent(component.curvedMirror(R_PR3, asX.component('PR2').z + dPR2_PR3, 'PR3'));
asX.addComponent(component.dielectric(R_BS, inf, dBSHR_BSAR, nSRM, asX.component('PR3').z + dPR3_BSHR, 'BSHR'));

% Build Y Path
asY = asX.duplicate;  

% X - Path from BSHR to ITXHR
asX.addComponent(component.dielectric(R_BS, inf, tCPX, nCPX, asX.component('BSHR').z + dBSAR_CPX1, 'CPXf'));
asX.addComponent(component.curvedMirror(-R_ITMX/nITMX, asX.component('CPXf').z + dCPX2_ITXAR + tITMX/nITMX, 'ITXHR'));
% X - Path back from ITXHR to BS HR
asX.addComponent(component.dielectric(inf, inf, 100e-3, nCPX, asX.component('ITXHR').z + dCPX2_ITXAR, 'CPXb'));
asX.addComponent(component.dielectric(inf, -R_BS, 2*dBSHR_BSAR, nSRM, asX.component('CPXb').z + dBSAR_CPX1, 'BSAR'));
% X - Path back from BS AR to SRM
asX.addComponent(component.curvedMirror(R_SR3, asX.component('BSAR').z + dSR3_BSAR, 'SR3'));
asX.addComponent(component.curvedMirror(R_SR2, asX.component('SR3').z + dSR2_SR3 , 'SR2'));
asX.addComponent(component.dielectric(-R_SRM, inf, tSRM, nSRM, asX.component('SR2').z + dSRM_SR2, 'SRM'));


% Y - Path from BSHR to ITYHR   
asY.addComponent(component.dielectric(R_CPY, inf, tCPY, nCPY, asY.component('BSHR').z + dBSHR_CPY1, 'CPYf'));
asY.addComponent(component.curvedMirror(-R_ITMY/nITMY, asY.component('CPYf').z + dCPY2_ITYAR + tITMY/nITMY, 'ITYHR'));
% Y - Path back from ITYHR to BS AR
asY.addComponent(component.dielectric(inf, -R_CPY, 100e-3, nCPX, asY.component('ITYHR').z + dCPY2_ITYAR, 'CPYb'));
asY.addComponent(component.dielectric(inf, -R_BS, dBSHR_BSAR, nSRM, asY.component('CPYb').z + dBSHR_CPY1, 'BSAR')); 
% Y - Path back from BS AR to SRM
asY.addComponent(component.curvedMirror(R_SR3, asY.component('BSAR').z + dSR3_BSAR, 'SR3'));
asY.addComponent(component.curvedMirror(R_SR2, asY.component('SR3').z + dSR2_SR3 , 'SR2'));
asY.addComponent(component.dielectric(-R_SRM, inf, tSRM, nSRM, asY.component('SR2').z + dSRM_SR2, 'SRM'));



% Overlap between beams after SRM
qX = asX.qPropagate([asX.component('SRM').z + 100e-2]).q ;
qY = asY.qPropagate([asY.component('SRM').z + 100e-2]).q ;

overlapXY = overlap(beamq(qX), beamq(qY));
disp(['overlap of Qs after SRM = ',num2str(overlapXY)])


% Save a copy for the ASAIR Path
nominalASAIR = asX.duplicate;

% OutputTelescope - X path
asX.addComponent(component.curvedMirror(R_OM1, asX.component('SRM').z + dSRM_OM1-tSRM, 'OM1'));
asX.addComponent(component.curvedMirror(R_OM2, asX.component('OM1').z + dOM1_OM2, 'OM2'));

% Save a copy for the AS WFS Path
nominalWFSX = asX.duplicate;

asX.addComponent(component.flatMirror(asX.component('OM2').z + dOM2_OM3, 'OM3'));
asX.addComponent(component.flatMirror(asX.component('OM3').z + dOM3_OMC, 'OMC waist'));
asX.addComponent(component.flatMirror(asX.component('OMC waist').z + 0.0434, 'OMC A'));
asX.addComponent(component.flatMirror(asX.component('OMC waist').z + 0.4845, 'OMC B'));


% OutputTelescope - Y path
asY.addComponent(component.curvedMirror(R_OM1, asY.component('SRM').z + dSRM_OM1-tSRM, 'OM1'));
asY.addComponent(component.curvedMirror(R_OM2, asY.component('OM1').z + dOM1_OM2, 'OM2'));

nominalWFSY = asY.duplicate;

asY.addComponent(component.flatMirror(asY.component('OM2').z + dOM2_OM3, 'OM3'));
asY.addComponent(component.flatMirror(asY.component('OM3').z + dOM3_OMC, 'OMC waist'));
asY.addComponent(component.flatMirror(asY.component('OMC waist').z + 0.0434, 'OMC A'));
asY.addComponent(component.flatMirror(asY.component('OMC waist').z + 0.4845, 'OMC B'));



% Mode match calculation to real / nominal OMC
w_OMC = 0.490e-3; % in meters, as-built OMC parameters from Koji

asX.targetWaist(w_OMC, asX.component('OM3').z + dOM3_OMC);
asY.targetWaist(w_OMC, asY.component('OM3').z + dOM3_OMC);

overlapXOMC = asX.targetOverlap;
overlapYOMC = asY.targetOverlap;

disp(['overlap of X at OMC = ',num2str(overlapXOMC)]);
disp(['overlap of Y at OMC = ',num2str(overlapYOMC)]);

X_OMCA_pos = asX.qPropagate(asX.component('OMC A').z);
X_OMCB_pos = asX.qPropagate(asX.component('OMC B').z);

disp(['nominal ITMX beam waist at OMC A = ',num2str(X_OMCA_pos.beamWidth)]);
disp(['nominal ITMX beam waist at OMC B = ',num2str(X_OMCB_pos.beamWidth)]);

Y_OMCA_pos = asY.qPropagate(asY.component('OMC A').z);
Y_OMCB_pos = asY.qPropagate(asY.component('OMC B').z);

disp(['nominal ITMY beam waist at OMC A = ',num2str(Y_OMCA_pos.beamWidth)]);
disp(['nominal ITMY beam waist at OMC B = ',num2str(Y_OMCB_pos.beamWidth)]);


disp(['nominal ITMX Guoy Phase at OM1 = ',num2str(asX.gouyPhase(asX.component('OM1').z))]);
disp(['nominal ITMX Guoy Phase at OM2 = ',num2str(asX.gouyPhase(asX.component('OM2').z))]);
disp(['nominal ITMX Guoy Phase at OM3 = ',num2str(asX.gouyPhase(asX.component('OM3').z))]);

disp(['nominal ITMY Guoy Phase at OM1 = ',num2str(asY.gouyPhase(asY.component('OM1').z))]);
disp(['nominal ITMY Guoy Phase at OM2 = ',num2str(asY.gouyPhase(asY.component('OM2').z))]);
disp(['nominal ITMY Guoy Phase at OM3 = ',num2str(asY.gouyPhase(asY.component('OM3').z))]);




% Construct the nominal AS WFS path - parameters from T1200410-v2 table 2
% Beam is transmitted through the flat mirror OM3
% L101 is 1-inch lens on WFS sled with f=0.334m
% Distances from OM2 to AS WFS sled - these numbers were calculated using as-built measurements from photographs
dOM3_L101 = 0.6376;
dL101_ASWFSA = 0.188;
dL101_ASWFSB = 0.387;

nominalWFSX.addComponent(component.dielectric(inf, inf, tOM1, nOM1, nominalWFSX.component('OM2').z + dOM2_OM3, 'OM3'));
nominalWFSX.addComponent(component.lens(0.334, nominalWFSX.component('OM3').z + dOM3_L101, 'L101'));
nominalWFSX.addComponent(component.flatMirror(nominalWFSX.component('L101').z + dL101_ASWFSA, 'WFSA'));
nominalWFSX.addComponent(component.flatMirror(nominalWFSX.component('L101').z + dL101_ASWFSB, 'WFSB'));

WFSA_posX = nominalWFSX.qPropagate(nominalWFSX.component('WFSA').z);
WFSB_posX = nominalWFSX.qPropagate(nominalWFSX.component('WFSB').z);

disp(['nominal ITMX beam waist at WFSA = ',num2str(WFSA_posX.beamWidth)]);
disp(['nominal ITMX beam waist at WFSB = ',num2str(WFSB_posX.beamWidth)]);




nominalWFSY.addComponent(component.dielectric(inf, inf, tOM1, nOM1, nominalWFSY.component('OM2').z + dOM2_OM3, 'OM3'));
nominalWFSY.addComponent(component.lens(0.334, nominalWFSY.component('OM3').z + dOM3_L101, 'L101'));
nominalWFSY.addComponent(component.flatMirror(nominalWFSY.component('L101').z + dL101_ASWFSA, 'WFSA'));
nominalWFSY.addComponent(component.flatMirror(nominalWFSY.component('L101').z + dL101_ASWFSB, 'WFSB'));

WFSA_posY = nominalWFSY.qPropagate(nominalWFSY.component('WFSA').z);
WFSB_posY = nominalWFSY.qPropagate(nominalWFSY.component('WFSB').z);

disp(['nominal ITMY beam waist at WFSA = ',num2str(WFSA_posY.beamWidth)]);
disp(['nominal ITMY beam waist at WFSB = ',num2str(WFSB_posY.beamWidth)]);



zWFS = nominalWFSX.component('OM2').z-0.1:0.002:nominalWFSX.component('WFSB').z+0.05;
zOMC = asX.component('OM1').z-0.1:0.002:asX.component('OMC waist').z+0.6;

figure(101)
asX.plotSummary(zOMC);

figure(102)
asY.plotSummary(zOMC);

figure(103)
nominalWFSX.plotSummary(zWFS)

figure(104)
nominalWFSY.plotSummary(zWFS)


