clear all
close all
% The script first calculates the deflection angle caused by the ETM and ERM
% wedge for both green and red beam.
% It also calculates the ABCD matrix from the ETM to the two red QPDs.
% Using these two, it tells the beam position of the green beam reflected
% off of the dichroic (which is very small in power) on QPD1 and QPD2.

%% Telescope path parameters.
d_ETM_Pri = 47*25.4e-3; % physical distance from ETM HR side to primary is about 48 inches.
                        % The exact number doesn't matter that much, but
                        % anyway this was obtained from Figure 1 and 22 of TMS
                        % preliminary design document E100040 using ISC table 
                        % dimension (30x30").
zETM = 0;

ROC_Pri=4; %primary focal length 2m
zPri=zETM + d_ETM_Pri;

ROC_Sec = -0.2; %secondary focal length -10cm;
detuning = 2.58636e-3; %nominal detuning parameter;

zSec = zPri + (ROC_Pri+ROC_Sec)/2+detuning;

dSec_L101 = 70*25.4e-3; % From the secondary to the first lens on the red QPD sled is about 70 inches, see e.g. G1100873 for the on-table distances and D1000484 for the side view.
zL101 = zSec + dSec_L101;
fL101 = 333.6e-3; % This is the standard in-vac 2" ISC lens.

dL101_L102 = 240e-3; % from T1000247
zL102 = zL101 + dL101_L102;
fL102 = -111e-3;

dL102_QPD1 = 410e-3; % from T1000247
zQPD1 = zL102 + dL102_QPD1;

dL102_QPD2 = 710e-3;  % from T1000247
zQPD2 = zL102 + dL102_QPD2;

%% define the beam path
mypath = beamPath();
mypath.addComponent(component.flatMirror(zETM, 'ETM')); % curvature of the ETM doesn't matter because the wedge effect is after the curved surface.
mypath.addComponent(component.curvedMirror(ROC_Pri, zPri, 'Primary'));
mypath.addComponent(component.curvedMirror(ROC_Sec, zSec, 'Secondary'));
mypath.addComponent(component.lens(fL101, zL101, 'L101'));
mypath.addComponent(component.lens(fL101, zL102, 'L102'));
mypath.addComponent(component.flatMirror(zQPD1, 'QPD1'));
mypath.addComponent(component.flatMirror(zQPD2, 'QPD2'));

%% deflection angle
nR = 1.44963; % refractive index of fused silica for 1.064um, Handbook of optics
nG = 1.46071; % for green, Handbook of Optics


wedgeETM = 0.08 *pi/180; % Vertical wedge of ETM. 
                         % Nominal wedge is 0.07 +0.03 -0.0, so it tends to
                         % be bigger.
                         % In reality this was measured to be between 0.07
                         % and 0.08 degrees depending on the optic.

wedgeERM = 0.04 * pi/180; % Horizontal wedge of ERM.
                          % No measurement, the spec says 0.04 + 0.04 -
                          % 0.03. In the worst case this should be 0.08,
                          % however it's more likely that the vendor
                          % decided to shoot for 0.04.
thetaRY = asin(nR*sin(wedgeETM)) - wedgeETM; % Red deflection angle due to ERM and ETM
thetaGY = asin(nG*sin(wedgeETM)) -wedgeETM;

dthetaY = thetaRY - thetaGY; 
% This should be a negative number as nG>nR.
% This means that the green beam is deflected down more. 

thetaRX = asin(nR*sin(wedgeERM)) - wedgeERM; % Red deflection angle due to ERM and ETM
thetaGX = asin(nG*sin(wedgeERM)) -wedgeERM;

dthetaX = thetaRX - thetaGX;

dtheta = [dthetaY, dthetaX];

%% Now we calculate the beam position on each of the optical components;
z=zeros(size(mypath.components));
d=zeros(length(z), 2);
for ii=1:length(z)
    z(ii) = mypath.component(ii).z;
    ABCD = mypath.getTransferMatrix(-1e-6, z(ii));
    for jj=1:2
        dummy = ABCD * [0; dtheta(jj)];
        d(ii, jj) = dummy(1);
    end
end

figure
plot(z, d(:, 1), 'r.-')
hold on
plot(z, d(:,2), 'b.-')
for ii = 1: length(z)
    text(z(ii), d(ii,2)+0.1e-3, mypath.component(ii).label);
end
grid on
title ('Position of the green beam in RED path')
ylabel('lateral position [m]')
xlabel('distance from ETM [m]');
legend('Vertical', 'Horizontal', 'location', 'southwest');
print('-depsc', 'redpathAlignment.eps')