T1000179 - v20

MaxDefocusAtITMAll = 2.1*10^-6;

H1-ITMX Vertex Hartmann Sensor
Optical Layout: [AirLens = 1.6994m, VacLens = -1.6994m]@ 820nm

Requirements
1. Image limiting aperture
2. Magnification = 1/20 x
3. Reference beam with own imaging device
4. Determine temperature sensitivity
5. Tolerancing on imaging optics

Get the Zemax Data - Import Units = mm, Convert to m

HAM Tables and BSCs

ConvertToM = 1/1000;

H1fileNAME = "/Users/aidanbrooks/Documents/LIGO/documents/Advanced_LIGO/TCS/Hartmann Sensor/design/optical layout and modeling/optical_paths";
SetDirectory[H1fileNAME];
H1ZEMAX = Import["D0901920-v6_H1_layout.txt", "TSV"];

HAM4L = Position[H1ZEMAX, "HAM4 ref"];
HAM4location = Take[H1ZEMAX[[HAM4L[[1]][[1]]]]];
HAM4XYZ = Reverse[HAM4location[[Range[6, 8]]]]*ConvertToM;


HAM5L = Position[H1ZEMAX, "HAM5 ref"];
HAM5location = Take[H1ZEMAX[[HAM5L[[1]][[1]]]]];
HAM5XYZ = Reverse[HAM5location[[Range[6, 8]]]]*ConvertToM;

BSC1L = Position[H1ZEMAX, "BSC1 REF"];
BSC1location = Take[H1ZEMAX[[BSC1L[[1]][[1]]]]];
BSC1XYZ = Reverse[BSC1location[[Range[6, 8]]]]*ConvertToM;

BSC2L = Position[H1ZEMAX, "BSC2 REF"];
BSC2location = Take[H1ZEMAX[[BSC2L[[1]][[1]]]]];
BSC2XYZ = Reverse[BSC2location[[Range[6, 8]]]]*ConvertToM;

BSC3L = Position[H1ZEMAX, "BSC3 REF"];
BSC3location = Take[H1ZEMAX[[BSC3L[[1]][[1]]]]];
BSC3XYZ = Reverse[BSC3location[[Range[6, 8]]]]*ConvertToM;


L1 = 1110.0*ConvertToM;
HAMTable = Table[Reverse[{0, L1 Cos[ii*Pi/3], L1 Sin[ii*Pi/3]}], {ii, 0.5, 6.5}];
BSCTable = Table[Reverse[{0, L1 Cos[ii*Pi/30], L1 Sin[ii*Pi/30]}], {ii, 0, 60}];
BSCTable = {{-676.275,630.9613938,0},{-228.6, 958.85,0},{228.6,958.85,0},{676.275,630.9613938,0},{676.275,-630.9613938, 0},{228.6,-958.85, 0},{-228.6,-958.85, 0},{-676.275,-630.9613938,0}}*ConvertToM;

HAM4G = Polygon[HAMTable + Table[HAM4XYZ, {ii, 1, 7}]];
HAM5G = Polygon[HAMTable + Table[HAM5XYZ, {ii, 1, 7}]];
BSC1G =Translate[ Rotate[Polygon[BSCTable], 90 Degree, {0, 0 ,1}], BSC1XYZ];
BSC2G = Translate[ Rotate[Polygon[BSCTable], 90 Degree, {0, 0 ,1}], BSC2XYZ];
BSC3G = Polygon[BSCTable + Table[BSC3XYZ, {ii, 1, Length[BSCTable]}]];

IFO = Graphics3D[{HAM4G, HAM5G, BSC1G, BSC2G, BSC3G} ,ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

SR2

SR2L = Position[H1ZEMAX, "SR2_RADIUS"];
SR2location = Take[H1ZEMAX[[SR2L[[1]][[1]]]]];
SR2XYZ =Reverse[SR2location[[Range[6, 8]]]]*ConvertToM;
SR2RRT  = SR2location[[Range[15, 17]]]*ConvertToM;
SR2ang = SR2location[[Range[9, 11]]]
SR2G = Cylinder[{SR2XYZ, SR2XYZ + SR2RRT[[3]]*Reverse[{0, Sin[2Pi - SR2ang[[1]]*Pi/180], Cos[ 2Pi - SR2ang[[1]]*Pi/180]}]}, SR2RRT[[1]]]
 Graphics3D[{SR2G, HAM4G},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

ITMXPO

ITMXPOL = Position[H1ZEMAX, "ITMXPO STEER M1"];
ITMXPOlocation = Take[H1ZEMAX[[ITMXPOL[[1]][[1]]]]];
ITMXPOXYZ =Reverse[ITMXPOlocation[[Range[6, 8]]]]*ConvertToM;
ITMXPORRT  = ITMXPOlocation[[Range[13, 15]]]*ConvertToM;
ITMXPORRT = {ITMXPORRT[[1]], ITMXPORRT[[3]], ITMXPORRT[[2]]};
ITMXPOang = ITMXPOlocation[[Range[9, 11]]]
ITMXPOG = Cylinder[{ITMXPOXYZ, ITMXPOXYZ + ITMXPORRT[[3]]*Reverse[{0, Sin[2Pi - ITMXPOang[[1]]*Pi/180], Cos[ 2Pi - ITMXPOang[[1]]*Pi/180]}]}, ITMXPORRT[[1]]]
 Graphics3D[{ITMXPOG, SR2G, HAM4G},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

SR3

SR3L = Position[H1ZEMAX, "SR3 REF"];
SR3location = Take[H1ZEMAX[[SR3L[[1]][[1]]]]];
SR3XYZ =Reverse[SR3location[[Range[6, 8]]]]*ConvertToM;

SR3L = Position[H1ZEMAX, "SR3_RADIUS"];
SR3location = Take[H1ZEMAX[[SR3L[[1]][[1]]]]];
SR3RRT  = SR3location[[Range[15, 17]]]*ConvertToM;
SR3ang = SR3location[[Range[9, 11]]];
SR3G = Cylinder[{SR3XYZ, SR3XYZ - SR3RRT[[3]]*Reverse[{0, Cos[2Pi - SR3ang[[1]]*Pi/180], Sin[ 2Pi - SR3ang[[1]]*Pi/180]}]}, SR3RRT[[1]]]
 Graphics3D[{SR3G, HAM5G},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

BS

BSL = Position[H1ZEMAX, "BS"];
BSlocation = Take[H1ZEMAX[[BSL[[1]][[1]]]]];
BSXYZ =Reverse[BSlocation[[Range[6, 8]]]]*ConvertToM;
BSRRT  = BSlocation[[Range[13, 14]]]*ConvertToM;
BSang = BSlocation[[Range[9, 11]]];
BSG = Cylinder[{BSXYZ, BSXYZ - BSRRT[[2]]*Reverse[{0, Cos[-Pi/4], Sin[-Pi/4]}]}, BSRRT[[1]]]
 Graphics3D[{BSG, BSC2G, Blue, Sphere[BSXYZ, 25*ConvertToM]},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

ITMY

ITMYL = Position[H1ZEMAX, "ITMY HR REF"];
ITMYlocation = Take[H1ZEMAX[[ITMYL[[1]][[1]]]]];
ITMYXYZ =Reverse[ITMYlocation[[Range[6, 8]]]]*ConvertToM;

ITMYL = Position[H1ZEMAX, "ITMY"];
ITMYlocation = Take[H1ZEMAX[[ITMYL[[1]][[1]]]]];
ITMYRRT  = {170, 200}*ConvertToM;
ITMYang = ITMYlocation[[Range[9, 11]]];
ITMYG = Cylinder[{ITMYXYZ, ITMYXYZ - ITMYRRT[[2]]*{0, 1, 0}}, ITMYRRT[[1]]];
 Graphics3D[{ITMYG, BSC1G, Blue, Sphere[ITMYXYZ, 25*ConvertToM]},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

ITMX

ITMXL = Position[H1ZEMAX, "ITMX HR REF"];
ITMXlocation = Take[H1ZEMAX[[ITMXL[[1]][[1]]]]];
ITMXXYZ =Reverse[ITMXlocation[[Range[6, 8]]]]*ConvertToM;

ITMXL = Position[H1ZEMAX, "ITMX"];
ITMXlocation = Take[H1ZEMAX[[ITMXL[[1]][[1]]]]];
ITMXRRT  = {170, 200}*ConvertToM;
ITMXang = ITMXlocation[[Range[9, 11]]];
ITMXG = Cylinder[{ITMXXYZ, ITMXXYZ - ITMXRRT[[2]]*{1, 0, 0}}, ITMXRRT[[1]]];
 Graphics3D[{ITMXG, BSC3G, Blue, Sphere[ITMXXYZ, 25*ConvertToM]},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

Viewports

HAM4VP6L = Position[H1ZEMAX, "HAM4 VP6"];
HAM4VP6location = Take[H1ZEMAX[[HAM4VP6L[[1]][[1]]]]];
HAM4VP6XYZ =(Reverse[HAM4VP6location[[Range[6, 8]]]]*ConvertToM + HAM4XYZ);
HAM4VP6G = Cylinder[{{0,0,0}, {0,0,0} - {12.7, 0, 0}}*ConvertToM, 63.5*ConvertToM];
HAM4VP6G = Rotate[HAM4VP6G, 13.92 Degree, {0, 0, 1}];
HAM4VP6G = Translate[HAM4VP6G, HAM4VP6XYZ];
 Graphics3D[{HAM4VP6G, HAM4G},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

Hartmann BS Optics

HSBS1L = Position[H1ZEMAX, "HARTMANN-ITMY STEER M1"];
HSBS1location = Take[H1ZEMAX[[HSBS1L[[1]][[1]]]]];
HSBS1XYZ =Reverse[HSBS1location[[Range[6, 8]]]] *ConvertToM;
HSBS1G = Sphere[HSBS1XYZ, 20*ConvertToM];

HSBS2L = Position[H1ZEMAX, "ITMY PO STEER M1"];
HSBS2location = Take[H1ZEMAX[[HSBS2L[[1]][[1]]]]];
HSBS2XYZ =Reverse[HSBS2location[[Range[6, 8]]]] *ConvertToM;
HSBS2G = Sphere[HSBS2XYZ, 20*ConvertToM];

HSBS2L = Position[H1ZEMAX, "ITMY PO TEL M1"];
HSBS2location = Take[H1ZEMAX[[HSBS2L[[1]][[1]]]]];
HSBS3XYZ =Reverse[HSBS2location[[Range[6, 8]]]]*ConvertToM ;
HSBS3G = Sphere[HSBS3XYZ, 20*ConvertToM];

HSBS2L = Position[H1ZEMAX, "ITMY PO TEL M2"];
HSBS2location = Take[H1ZEMAX[[HSBS2L[[1]][[1]]]]];
HSBS4XYZ =Reverse[HSBS2location[[Range[6, 8]]]] *ConvertToM;
HSBS4G = Sphere[HSBS4XYZ, 20*ConvertToM];

HSBS2L = Position[H1ZEMAX, "HARTMANN-ITMY BS STEER"];
HSBS2location = Take[H1ZEMAX[[HSBS2L[[1]][[1]]]]];
HSBS5XYZ =Reverse[HSBS2location[[Range[6, 8]]]]*ConvertToM ;
HSBS5G = Sphere[HSBS5XYZ, 20*ConvertToM];

 Graphics3D[{HSBS5G, HAM4G},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

IFO

Graphics3D[{HAM4VP6G,ITMXG,ITMYG, BSG, SR3G, SR2G, BSC1G, BSC2G, BSC3G, HAM4G, HAM5G, Blue, HSBS1G, HSBS2G, HSBS3G, HSBS4G, HSBS5G},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

HAM4

Graphics3D[{HAM4VP6G, SR2G, HAM4G,Blue, HSBS1G, HSBS2G, HSBS3G, HSBS4G, HSBS5G}, ViewPoint -> {0, 0, Infinity}, ViewVertical -> {0, 0, 1}]

PastSR2 = Norm[HSBS1XYZ - SR2XYZ]+Norm[HSBS1XYZ - HAM4VP6XYZ]

Optical Layout - X Arm
Requirements
1. HWS = Conjugate Plane of BS
2. Beam Size at ITM = Sqrt[2]*ModeSize
3. Magnification = 17.5x
4. Put the telescope outside the vacuum - Can't do it all
5. Beam is approximately collimated at the HWS

Set probe beam size at ITM = 17.5*Diameter_of_CCD/SQRT[2]

BeamSizeAtITM= 127/1000

nFusedSilica840 = 1.445;
CPthick = 130*ConvertToM;
n1 = 1.0;
TotalDeMagnification = -17.5;

Create Function to Determine Output Radius and Beam Size

OutRW[wIN_, RIN_, λIN_,ABCD_]:= Module[{},
Q1Inv = (1/RIN)- ((I λIN)/(Pi wIN^2));
Q1a = 1/Q1Inv;
Q2temp = ((ABCD[[1]][[1]]Q1a + ABCD[[1]][[2]])/(ABCD[[2]][[1]]Q1a + ABCD[[2]][[2]]));
Q2Inv = 1/Q2temp;
ROut = Re[Q2Inv]^-1;
wOut = Sqrt[((-λIN)/(Pi Im[Q2Inv]))];
Return[{ROut, wOut}];
];
FindNearestWaist[wIN_, RIN_, λIN_]:= Module[{},
Q1Inv = (1/RIN)- ((I λIN)/(Pi wIN^2));
Q1a = 1/Q1Inv;


Q2temp = Q1a;
Q2tempR = Re[Q2temp]+x0;
Q2tempI = Im[Q2temp];
Q2InvR = (Q2tempR/(Q2tempR^2 + Q2tempI^2));


ROut = Q2InvR^(-1);

solnIN = Solve[1/ROut == 0, x0];


Return[x0/.solnIN[[1]]];
];

Optical Path - determine solution: AirLens = 1.6994m, VacLens = -1.6994m

Given Physical Parameters

Distances

ITMXPOtoSR3 = Norm[SR3XYZ-ITMXPOXYZ]
SR3toITMX = Norm[SR3XYZ - BSXYZ] + Norm[ITMXXYZ - BSXYZ]  - 2*BSRRT[[2]]+ (nFusedSilica840-1)*(ITMYRRT[[2]] + CPthick)
SR3toBS = Norm[SR3XYZ - BSXYZ]

Curvatures - E0900095, T080268, E080518

ITMRoC = 1934/nFusedSilica840;
SR3RoC = 36;

ABCD Matrices - not relying on solution

Requirement 2: The mode from the ITM, including the size of Sqrt[2]*CavityModeSize is defined below

Test Gaussian beam in and out

MFromITM = (({{1, ITMXPOtoSR3}, {0, 1}})). (({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR3toITMX}, {0, 1}}));
MToITM =  (({{1, SR3toITMX}, {0, 1}})).(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, ITMXPOtoSR3}, {0, 1}}));

wIN = BeamSizeAtITM;
RIN = ITMRoC;
pIn = {RIN, wIN}
λIN = 840*10^-9;

pOut1 = OutRW[wIN, RIN, λIN,MFromITM];
pOut2 = OutRW[pOut1[[2]], -pOut1[[1]], λIN,MToITM];
Print["Beam size at ITM = ", wIN*1000, " mm"];
Print["Beam size at ITMXPO = ", pOut1[[2]]*1000, " mm"];
Print["Beam size of retro-reflected ITMXPO beam at ITM = ", pOut2[[2]]*1000, " mm"];

ITM to ITMXPO: From HR surface outside ITM to AR surface outside ITMXPO [in m]

ITMToITMXPO =(({{1, ITMXPOtoSR3}, {0, 1}})).(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR3toITMX}, {0, 1}}));
MatrixForm[ITMToITMXPO]

BS to ITMXPO: ABCD Matrix - From BS to ITMXPO

BSToITMXPO =(({{1, ITMXPOtoSR3}, {0, 1}})).(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR3toBS}, {0, 1}}));
MatrixForm[ITMToITMXPO]

Solve for the best solution

MOut = Simplify[MToITM.(({{1, L3}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, L2}, {0, 1}})).(({{1, 0}, {-1/f1, 1}})).(({{1, L1a}, {0, 1}}))];
assumptions = {f2-> 1.0, f1-> 1.0}
MOut1 = MOut/. assumptions;
soln1 = Simplify[Solve[MOut1[[1]][[2]]==0 , dx]];
MOut2 = Simplify[MOut1 /. soln1[[1]]];
mag = 17.5;
NSolve[{MOut1[[1]][[2]]==0, MOut1[[1]][[1]]==mag, MOut1[[2]][[1]]==0}, {L1a, L2, L3}]

CVI Lenses - Available Focal Lengths (mm)

CVIFocalLengths  = {3.5,6,10,10,12,15,15,17,18,18,20,20,22,25,25,25,25,25,25,30,30,31,32,32,35,35,35,35,35,38,40,40,40,50,50,50,50,50,50,50,50,60,60,60,60,63,63,70,70,75,75,75,75,75,75,76,76,80,80,85,90,100,100,100,100,100,100,100,105,110,120,125,125,125,125,127,127,127,130,135,140,150,150,150,150,150,152,152,160,160,170,177,200,200,200,200,200,200,250,250,250,250,252,254,269,300,347,350,350,350,360,400,400,400,460,485,500,500,500,500,600,600,600,700,700,750,750,800,800,900,900,1000,1000,1000,1500,1500,2000,2000,3000,3000,4000,4000,5000,5000,7500,7500,10000,10000};
CVIFL2 = {};

For[i=1, i≤ Length[CVIFocalLengths], i++, {
addFL = CVIFocalLengths[[i]];

p1 = Position[CVIFL2, addFL];
If[Length[p1]== 0, CVIFL2 = Append[CVIFL2, addFL]];
}];

CVIFocalLengths = CVIFL2

CVI2 = Drop[CVIFocalLengths, 39]*10^-3.;

CVI3 = {0.2, 0.3, 0.4, 0.5`,0.75,1.`,1.5`, 2.`, 3.`};
Length[CVI2]^3
Length[CVI3]^3
CVI2 = CVI3;
For[ii = 1, ii≤ Length[CVI3], ii++, CVI2= Append[CVI2,CVI3[[ii]]*-1]];

CVILenses820MM = {339.8,396.6,453.1,566.4,679.9,793.1,1133.0,1699.4,2265.9,3398.3,4531.0,5663.7,8496.7,11329.7,-226.6,-283.3,-339.8,-396.6,-454.4,-566.4,-679.9,-906.4,-1084.4,-1133.0,-1699.4}

CVILenses820 = CVILenses820MM/1000;

Loop through all options

MOut = Simplify[MToITM.(({{1, POtoVacLens}, {0, 1}})).(({{1, 0}, {-1/VacLens, 1}})).(({{1, VacLenstoAirLens}, {0, 1}})).(({{1, 0}, {-1/AirLens, 1}})).(({{1, AirLenstoHWS}, {0, 1}}))];
solnList = {};
ncount = 0;
dMag = 0.25;
For[mag = 17.5, mag≤ 17.5, mag = mag + dMag, {

For[mm = 1, mm≤ 2, mm++, {
mag1 = (-1)^mm*mag;

For[ii = 1, ii≤ Length[CVILenses820], ii++, {
For[jj= 1, jj≤ Length[CVILenses820], jj++, {

assumptions = {VacLens-> CVILenses820[[ii]], AirLens->CVILenses820[[jj]]};
MOut1 = MOut/. assumptions;
soln1  = Solve[{MOut1[[1]][[2]]==0, MOut1[[1]][[1]]==mag1, MOut1[[2]][[1]]==0}, {POtoVacLens, VacLenstoAirLens, AirLenstoHWS}];

For[kk=1, kk≤ Length[soln1], kk++, {
results = { AirLenstoHWS, VacLenstoAirLens,  POtoVacLens,  AirLens, VacLens, mag1, "AirLenstoHWS, VacLenstoAirLens,  POtoVacLens,  AirLens, VacLens, mag" }/. soln1[[kk]];
results = results/.assumptions;
If[(results[[1]]>0) && (results[[2]]> 2.7)&&(results[[3]]>0.5), {
If[(results[[1]]<3) && (results[[2]]< 5)&&(results[[3]]< 3), {
solnList = Append[solnList, results];
}];
}];
}];

ncount = ncount + 1;
If[Mod[ncount, 2500]==0, {
Print[{Length[solnList], ncount,21*2* Length[CVILenses820]^2}];
}];

}];
 
}];


}];
}];

MatrixForm[solnList]
H1XSoln = solnList;

MOut = Simplify[MToITM.(({{1, POtoVacLens}, {0, 1}})).(({{1, 0}, {-1/VacLens, 1}})).(({{1, VacLenstoAirLens}, {0, 1}})).(({{1, 0}, {-1/AirLens, 1}})).(({{1, AirLenstoHWS}, {0, 1}}))];
solnList = {};
ncount = 0;
dMag = 0.25;

For[CC = -0.01, CC≤ 0.01, CC= CC+ 0.001, {

mag1 = -17.5;


assumptions = {VacLens-> -1.6994, AirLens->1.6994};
MOut1 = MOut/. assumptions;
soln1  = Solve[{MOut1[[1]][[2]]==0, MOut1[[1]][[1]]==mag1, MOut1[[2]][[1]]==CC}, {POtoVacLens, VacLenstoAirLens, AirLenstoHWS}];

For[kk=1, kk≤ Length[soln1], kk++, {
results = { AirLenstoHWS, VacLenstoAirLens,  POtoVacLens,  AirLens, VacLens, mag1,CC, "AirLenstoHWS, VacLenstoAirLens,  POtoVacLens,  AirLens, VacLens, mag" }/. soln1[[kk]];
results = results/.assumptions;
If[(results[[1]]>0) && (results[[2]]> 2.7)&&(results[[3]]>0.5), {
If[(results[[1]]<3) && (results[[2]]< 5)&&(results[[3]]< 3), {
solnList = Append[solnList, results];
}];
}];
}];

}];


MatrixForm[solnList]

solnListCVI = solnList[[11]]

L3sub = (664+659)+(664-177) + (650-177)
MOut = Simplify[MToITM.(({{1, L3}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, L2}, {0, 1}})).(({{1, 0}, {-1/f1, 1}})).(({{1, L1a}, {0, 1}}))]/. {L3-> solnListCVI[[3]], f2-> solnListCVI[[5]], L2-> solnListCVI[[2]]};
soln1 = Solve[MOut[[1]][[2]] == 0, L1a];
MOut1 = MOut/. soln1[[1]];
solvedvalues={f1, L1a, MOut1[[2]][[1]]}/. soln1;
mvalues = {15.0, 17.5, 20.0, -15.0, -17.5, -20.0};
solnListA = {};
For[ii=1, ii≤ Length[mvalues], ii++, {
msolve = mvalues[[ii]];
soln2 = Solve[MOut1[[1]][[1]]== msolve, f1];
solnListA = Append[solnListA, {msolve, {solvedvalues/. soln2[[1]]}}];

}];

Extract Optical Information - Conjugate plane, beam size at ITMXPO

Get Output Mode from ITMXPO

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840*10^-9;
OutRW[wIN, RIN, λIN,{{1, 100}, {0, 1}}];
pOut = OutRW[wIN, RIN, λIN,ITMToITMXPO];
OutR = pOut[[1]];
OutW = pOut[[2]];
OutModeITMXPO = pOut;
FindNearestWaist[pOut[[2]], pOut[[1]], λIN];

Must create mode - matching optics to produce this mode
Must create imaging optics to image the limiting aperture - the beam splitter - correctly

Print["Probe beam size at ITMXPO = ", OutW*1000, " mm."]
Print["Curvature at ITMXPO = ", OutR*1000, " mm."]

Demonstrate curvature sign by propagating a short distance further

pOutdz = OutRW[wIN, RIN, λIN,{{1, 0.1}, {0, 1}}.ITMToITMXPO];
OutRdz = pOutdz[[1]];
OutWdz = pOutdz[[2]];
Print["Probe beam size at ITMXPO = ", OutW*1000, " mm."]
Print["Probe beam size 100mm out of IFO from ITMXPO = ", OutWdz*1000, " mm."]

Analytic Solution for imaging a conjugate plane with an arbitrary telescope

M1 = Simplify[(({{1, z3}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z2}, {0, 1}})).(({{1, 0}, {-1/f1, 1}})).(({{1, z1}, {0, 1}})).(({{A11, B11}, {C11, D11}}))];

The final curvature is, roughly, the sum of the defocus of the telescope (fixed) and the curvature of the input beam (if, curvature ~= 0, then depends on position of telescope)
Therefore, easiest to vary z1 to set final C = 0

soln1 = Solve[M1[[2]][[1]] == 0, z1][[1]];
ASOLN1 = soln1;

The imaging is most directly solved by moving the final plane along the optical axis.
Therefore, easiest to vary z3 to set final B = 0

M1B= Simplify[M1[[1]][[2]]/. soln1];
soln2 =Simplify[ Solve[M1B == 0, z3][[1]]];
ASOLN2 = soln2;

This leaves z2 to solve for the magnification

M1A = Simplify[M1[[1]][[1]]/. soln1]/.soln2;
soln3 = Solve[M1A == m1a, z2][[1]];
ASOLN3 = soln3;
ZList =Simplify[ ((({z1, z2, z3}/. soln1)/. soln2)/. soln3)/. m1a-> 1/M];
ZListOriginal =Simplify[ ((({z1, z2, z3}/. soln1)/. soln2)/. soln3)];

C = 0 requirement

Simplify[soln1 /.soln3]

A = M requirement

Simplify[soln3]

B = 0 requirement

Simplify[soln2/. soln3]

ZList2 = Simplify[ZList /. {A11-> ITMToITMXPO[[1]][[1]], B11-> ITMToITMXPO[[1]][[2]],C11-> ITMToITMXPO[[2]][[1]],D11-> ITMToITMXPO[[2]][[2]]}]

Solve for the solution - ITMX: Magnification  = +ve!!

Just work out the solutions for the ABCD matrices as follows:
1. B = 0
2. C = 0
3. A = 1/17.5

eq1: Collimation requirement - at output in paraxial approximation
eq2: Total magnification requirement
eq3: Imaging requirement

Use the analytic solution determined above.
Notes: 
1. The combination of lenses was chosen by first constraining the parameter space of possible solutions to those where 
  a. {z1, z2, z3} >= {3, 5, 0}m 
  b. dS/dT of z2 <= 20x smaller than the maximum defocus error of 4E-4 (m^(-1))
2. Once this region was identified the lenses were chosen to have convenient values and to yield a solution close to the minimum Abs[dS/dT]

Hence: f1 = 1.0m, f2 = 3.0m

F1A = -2.0;
F2A = 2.5;

zsolnList = ZList2/. {f1-> F1A, f2-> F2A, M-> TotalDeMagnification}

AllSolns2 = { F1A, F2A, Round[zsolnList[[1]]*1000]/1000., Round[zsolnList[[2]]*1000]/1000., Round[zsolnList[[3]]*1000.0]/1000.0}
Z1A =solnListCVI[[3]];
Z2A = solnListCVI[[2]];
Z3A = solnListCVI[[1]];
F1A = solnListCVI[[5]];
F2A = solnListCVI[[4]];
TestSolutionSubstitution = {f1-> F1A, f2-> F2A, z1->Z1A, z2-> Z2A, z3-> Z3A}

ABCD Matrices - dependent on solution

ITMXPO to HWS: ABCD matrix - get the imaging matrix part - after ITMXPO

ImgMatrix = Simplify[(({{1, z3}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z2}, {0, 1}})).(({{1, 0}, {-1/f1, 1}})).(({{1, z1}, {0, 1}}))];
ITMXPOToHWS = ImgMatrix/. TestSolutionSubstitution;
MatrixForm[ITMXPOToHWS]

HWS to ITMXPO: ABCD matrix - get the imaging matrix part - after ITMXPO

ImgMatrix = Simplify[(({{1, z1}, {0, 1}})).(({{1, 0}, {-1/f1, 1}})).(({{1, z2}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z3}, {0, 1}}))];
HWSToITMXPO = ImgMatrix/. TestSolutionSubstitution;
MatrixForm[HWSToITMXPO]

ITM to HWS: ABCD matrix coming out of the vacuum - For Gaussian beam mode checking

ITMToHWS = ITMXPOToHWS.ITMToITMXPO;
MatrixForm[ITMToHWS]
HWSToITM =(({{1, SR3toITMX}, {0, 1}})) .(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, ITMXPOtoSR3}, {0, 1}})).(({{1, z1 }, {0, 1}})).(({{1, 0}, {-1/f1, 1}})). (({{1, z2}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z3}, {0, 1}}))/. TestSolutionSubstitution;
MITM = (({{1, 0}, {((-2)/(-ITMRoC)), 1}}));
MHWSv = ITMToHWS.MITM.HWSToITM;

BS to HWS: ABCD matrix - Beam splitter to HWS - Shows Imaging

BS2HWS = ITMXPOToHWS.BSToITMXPO;
MatrixForm[BS2HWS]
MatrixForm[ITMToHWS]

Requirement 1: The Beam Splitter to HWS ABCD matrix has a value of B approximately equal to 0.
Requirement 5: Beam is approximately collimated at the HWS

Do a Gaussian beam propagation

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOut = OutRW[wIN, RIN, λIN,ITMToHWS];
OutR = pOut[[1]];
OutW = pOut[[2]];
GaussianMag = wIN/OutW;
Print["Gaussian beam size DeMagnification = ", GaussianMag, "x"]
Print["Beam RoC at HWS = ", OutR, "m"]

Requirement 3: The beam size at the HWS is 17.43x smaller than at the beam splitter
Requirement 5: The beam is approximately flat at the HWS

Define the Input Mode

InputMode = {OutW, -OutR}
InputModeFromHWS = InputMode;

Propagate modes - from ITM to outside and then back again

Propagate the ITM mode from the ITM to the HWS

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
Print["Probe beam size at ITM = ", wIN*1000, " mm."]
Print["Curvature at ITM = ", RIN, " m."]
pOut = OutRW[wIN, RIN, λIN,ITMToHWS];
Print["Probe beam size at HWS = ", pOut[[2]]*1000, " mm."]
Print["Curvature at HWS = ", pOut[[1]]*1000, " mm."]
Print["ITM_Beam_Size/HWS_Beam_Size = ", wIN/pOut[[2]], " x."]

Propagate the HWS mode to the ITM

pOut = OutRW[pOut[[2]], -pOut[[1]], λIN,HWSToITM];
Print["Probe beam size at ITM = ", pOut[[2]]*1000, " mm."]
Print["Curvature at ITM = ", pOut[[1]], " m."]

Propagate from HWS mode to HWS

Print["Probe beam size at ITM = ", InputMode[[1]]*1000, " mm."]
Print["Curvature at ITM = ", InputMode[[2]], " m."]
pOut = OutRW[InputMode[[1]], InputMode[[2]], λIN,MHWSv];
Print["Probe beam size at HWS = ", pOut[[2]]*1000, " mm."]
Print["Curvature at HWS = ", pOut[[1]], " m."]

Plot the Solution

Plot solution from ITMXPO to HWS

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;

p1 = Plot[(OutRW[wIN, RIN, λIN,(({{1, x}, {0, 1}})).ITMToITMXPO])[[2]]*1000, {x, 0, Z1A}, PlotStyle-> {Red,Thickness[0.01]}, PlotRange-> {{0, 14}, {0, 40}}];
BS1 = Table[{x,(OutRW[wIN, RIN, λIN,(({{1, x}, {0, 1}})).ITMToITMXPO])[[2]]*1000},{x, 0, Z1A, Z1A/200}];

p2 = Plot[(OutRW[wIN, RIN, λIN,(({{1, x-Z1A}, {0, 1}})).(({{1, 0}, {-1/F1A, 1}})).(({{1, Z1A}, {0, 1}})).ITMToITMXPO])[[2]]*1000, {x, Z1A, Z1A+Z2A}, PlotStyle-> {Blue,Thickness[0.01]}];
BS2 =  Table[{x,(OutRW[wIN, RIN, λIN,(({{1, x-Z1A}, {0, 1}})).(({{1, 0}, {-1/F1A, 1}})).(({{1, Z1A}, {0, 1}})).ITMToITMXPO])[[2]]*1000}, {x, Z1A, Z1A+Z2A, (Z2A/200)}];
p3 = Plot[(OutRW[wIN, RIN, λIN,(({{1, x-Z1A-Z2A}, {0, 1}})).(({{1, 0}, {-1/F2A, 1}})).(({{1, Z2A}, {0, 1}})).(({{1, 0}, {-1/F1A, 1}})).(({{1, Z1A}, {0, 1}})).ITMToITMXPO])[[2]]*1000, {x, Z1A+Z2A, Z1A+Z2A+Z3A}, PlotStyle-> {Red,Thickness[0.01]}];

BS3 =  Table[{x,(OutRW[wIN, RIN, λIN,(({{1, x-Z1A-Z2A}, {0, 1}})).(({{1, 0}, {-1/F2A, 1}})).(({{1, Z2A}, {0, 1}})).(({{1, 0}, {-1/F1A, 1}})).(({{1, Z1A}, {0, 1}})).ITMToITMXPO])[[2]]*1000}, {x, Z1A+Z2A, Z1A+Z2A+Z3A, Z3A/200}];
ITMXPlottedSolutionFromITMXPO = Show[p1, p2, p3, Frame-> True, FrameLabel-> {"Propagation Distance from ITMXPO (m)", "Beam size (mm)", "Beam propagating from ITMXPO to HWS"}]

BS = BS1;
For[ii=1, ii≤ Length[BS2], ii++, {
BS = Append[BS, BS2[[ii]]];
}];
For[ii=1, ii≤ Length[BS3], ii++, {
BS = Append[BS, BS3[[ii]]];
}];
Export["/Advanced_LIGO/TCS/Hartmann sensor/H1_HWSX_beam_size.txt", BS, "Table"];

Plot from ITMX to ITMXPO

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
plotR = {{0, 55}, {0, 200}};
p1 = Plot[(OutRW[wIN, RIN, λIN,(({{1, x}, {0, 1}}))])[[2]]*1000, {x, 0, SR3toITMX}, PlotStyle-> {Red,Thickness[0.01]}, PlotRange-> plotR];
p2 = Plot[(OutRW[wIN, RIN, λIN,(({{1, x-SR3toITMX}, {0, 1}})).(({{1, 0}, {-2/SR3RoC, 1}})).(({{1, SR3toITMX}, {0, 1}}))])[[2]]*1000, {x, SR3toITMX, SR3toITMX + ITMXPOtoSR3}, PlotStyle-> {Blue,Thickness[0.01]}];
p3 = Plot[(OutRW[wIN, RIN, λIN,(({{1, x-SR3toITMX-ITMXPOtoSR3}, {0, 1}})).ITMToITMXPO])[[2]]*1000, {x,  SR3toITMX + ITMXPOtoSR3,SR3toITMX + ITMXPOtoSR3+ Z1A}, PlotStyle-> {Red,Thickness[0.01]}, PlotRange->plotR ];
ITMXPlottedSolutionFromITMX = Show[p1, p2, p3, Frame-> True, FrameLabel-> {"Propagation Distance from ITMX (m)", "Beam size (mm)", "Beam propagating from ITMX to F1"}]

Tolerancing on X Arm solution

Get the Variable matrix

SR3toHWSVar =(({{1, Z3A+dz3}, {0, 1}})).(({{1, 0}, {-1/F2A -(dSf2), 1}})).(({{1, Z2A+dz2}, {0, 1}})).(({{1, 0}, {-1/F1A -(dSf1), 1}})).(({{1, Z1A+dz1}, {0, 1}})).(({{1, ITMXPOtoSR3}, {0, 1}})).(({{1, 0}, {((-2)/SR3RoC), 1}}));
BStoHWSVar = SR3toHWSVar.(({{1, SR3toBS}, {0, 1}}));
ITMXtoHWSVar = SR3toHWSVar.(({{1, SR3toITMX}, {0, 1}}));

Variation in imaging (B = 0) with length changes: (assume that FL variations can be compensated for) - NO ISSUES - Easy to work with!

M1 = BStoHWSVar.(({{1, x}, {0, 1}}));
B1 = M1[[1]][[2]];

Test variation in position of F1 - solve for the position of the conjugate plane - INSENSITIVE

Set all delta values = 0 except for dz2 and dz1. Set these last two to the same value with opposite sign to model a lens out of position

B1A = Simplify[B1/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> X1, dz1-> -X1}];
Solve[B1A==0, x];
ConjPlaneDisplacement = x/. %[[1]];
NominalConjugatePlane = ConjPlaneDisplacement/. X1-> 0;
Print["Conjugate Plane is nominally out of position by ", NominalConjugatePlane, "m"];
RateCPD = (D[ConjPlaneDisplacement, X1])/. X1-> 0;
Print["Conjugate Plane displacement for +/- 5 mm displacement of F1 = ", RateCPD*1000*0.005, "mm"];
BTestF1OutOfPosition = {RateCPD*1000*0.005, "mm", 5, "mm"};

Hence the conjugate plane is mostly invariant to the position of F1

Test variation in position of F2 - solve for the position of the conjugate plane - INSENSITIVE to distance between F1 and F2

Set all delta values = 0 except for dz2 and dz3. Set these last two to the same value with opposite sign to model a lens out of position

B1A = Simplify[B1/. {dz1-> 0, dSf2-> 0, dSf1-> 0, dz2-> -X1, dz3-> X1}];
Solve[B1A==0, x];
ConjPlaneDisplacement = x/. %[[1]];
NominalConjugatePlane = ConjPlaneDisplacement/. X1-> 0;
Print["Conjugate Plane is nominally out of position by ", NominalConjugatePlane, "m"];
RateCPD = (D[ConjPlaneDisplacement, X1])/. X1-> 0;
Print["Conjugate Plane displacement for +/- 5 mm displacement of F1 = ", RateCPD*0.005, "m"];
BTestF2OutOfPosition = {RateCPD*1000*0.005, "mm", 5, "mm"};

Test variation in F1-F2 distance - INSENSITIVE

B1A = Simplify[B1/. {dz1-> 0, dSf2-> 0, dSf1-> 0, dz2-> -X1, dz3-> 0}];
Solve[B1A==0, x];
ConjPlaneDisplacement = x/. %[[1]];
NominalConjugatePlane = ConjPlaneDisplacement/. X1-> 0;
Print["Conjugate Plane is nominally out of position by ", NominalConjugatePlane, "m"];
RateCPD = (D[ConjPlaneDisplacement, X1])/. X1-> 0;
Print["Conjugate Plane displacement for +/- 5 mm variation in distance between F1 and F2 = ", RateCPD*0.005, "m"];
BTestF1F2Distance = {RateCPD*1000*0.005, "mm", 5, "mm"};

Hence the conjugate plane is totally insensitive to the distance between F2 and F1

Test variation in position of HWS - SENSITIVE = GOOD for tuning. Possibly bad for fringes

Set all delta values = 0 except for dz3.

B1A = Simplify[B1/. {dz1-> 0, dSf2-> 0, dSf1-> 0, dz2-> 0, dz3-> X1}];
Solve[B1A==0, x];
ConjPlaneDisplacement = x/. %[[1]];
NominalConjugatePlane = ConjPlaneDisplacement/. X1-> 0;
Print["Conjugate Plane is nominally out of position by ", NominalConjugatePlane, "m"];
RateCPD = (D[ConjPlaneDisplacement, X1])/. X1-> 0;
Print["Conjugate Plane displacement for +/- 5 mm displacement of F1 = ", RateCPD*0.005, "m"];
BTestF2HWSDistance = {RateCPD*1000*0.005, "mm", 5, "mm"};

Hence the conjugate plane is totally sensitive to the position of HWS - THIS IS A GOOD THING - means it is easy to tune.

Variation in magnification (A = 1/17.5) with length changes: (assume that FL variations can be compensated for) - NO ISSUES - Easy to work with!

M1 =ITMXtoHWSVar;
A1 = M1[[1]][[1]];

Test variation in position of F1 - determine the magnification - INSENSITIVE

Set all delta values = 0 except for dz2 and dz1. Set these last two to the same value with opposite sign to model a lens out of position

A1A = Simplify[A1/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> X1, dz1-> -X1}];
NominalMagnification= 1/A1A/. X1-> 0;
Print["Nominal paraxial magnification is ", NominalMagnification, "m"];
RateVariationInMag = (D[A1A, X1])/. X1-> 0;
Print["Variation in magnification for +/- 5 mm displacement of F1 = ", RateVariationInMag*0.005, "x"];

Hence the paraxial magnification is VERY invariant to the position of F1

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMXtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> 0, dz1-> -0})];
BeamWAtHWSNom = pOutNom[[2]];
MOut = (ITMXtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> X1, dz1-> -X1});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamWAtHWSX1= pOutX1[[2]];
MagNom = BeamSizeAtITM/BeamWAtHWSNom;
MagX1 = BeamSizeAtITM/BeamWAtHWSX1;
Print["Absolute Variation in magnification for +5mm displacement = ", MagX1-MagNom, "x"];
Print["Relative Variation in magnification for +5mm displacement = ", ((MagX1-MagNom)/MagNom)*100, "%"];
ATestF1OutOfPosition = {Round[((MagX1-MagNom)/MagNom)*100*10]/10.0,"%",  5, "mm"};

Hence the Gaussian magnification is VERY invariant to the position of F1

Test variation in position of F2 - determine the magnification - INSENSITIVE to distance between F1 and F2

Set all delta values = 0 except for dL1 and dz3. Set these last two to the same value with opposite sign to model a lens out of position

A1A = Simplify[A1/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dz2-> -X1, dz1-> 0}];
NominalMagnification= 1/A1A/. X1-> 0;
Print["Nominal paraxial magnification is ", NominalMagnification, "m"];
RateVariationInMag = (D[A1A, X1])/. X1-> 0;
Print["Variation in magnification for +/- 5 mm displacement of F2 = ", RateVariationInMag*0.005, "x"];

Hence the paraxial magnification is VERY invariant to the position of F2

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMXtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> 0, dz1-> -0})];
BeamWAtHWSNom = pOutNom[[2]];
MOut = (ITMXtoHWSVar/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dz2-> -X1, dz1-> 0});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamWAtHWSX1= pOutX1[[2]];
MagNom = BeamSizeAtITM/BeamWAtHWSNom;
MagX1 = BeamSizeAtITM/BeamWAtHWSX1;
Print["Absolute Variation in magnification for +5mm displacement = ", MagX1-MagNom, "x"];
Print["Relative Variation in magnification for +5mm displacement = ", ((MagX1-MagNom)/MagNom)*100, "%"];
ATestF2OutOfPosition = {Round[((MagX1-MagNom)/MagNom)*100*10]/10.0,"%",  5, "mm"};

Hence the Gaussian magnification is marginally sensitive to the position of F2

Test variation in F1-F2 distance: determine the magnification - INSENSITIVE

A1A = Simplify[A1/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> -X1, dz1-> 0}];
NominalMagnification= 1/A1A/. X1-> 0;
Print["Nominal paraxial magnification is ", NominalMagnification, "m"];
RateVariationInMag = (D[A1A, X1])/. X1-> 0;
Print["Variation in magnification for +/- 5 mm displacement of F2 = ", RateVariationInMag*0.005, "x"];

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMXtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> 0, dz1-> -0})];
BeamWAtHWSNom = pOutNom[[2]];
MOut = (ITMXtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> -X1, dz1-> 0});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamWAtHWSX1= pOutX1[[2]];
MagNom = BeamSizeAtITM/BeamWAtHWSNom;
MagX1 = BeamSizeAtITM/BeamWAtHWSX1;
Print["Absolute Variation in mag. for +5mm distance change between F1 and F2 = ", MagX1-MagNom, "x"];
Print["Relative Variation in mag. for +5mm distance change between F1 and F2 = ", ((MagX1-MagNom)/MagNom)*100, "%"];
ATestF1F2Distance = {Round[((MagX1-MagNom)/MagNom)*100*10]/10.0,"%",  5, "mm"};

Hence the Gaussian magnification is very insensitive to the distance between F1 and F2

Test variation in position of HWS - determine the magnification - INSENSITIVE to position of HWS: 1% variation per 50mm

Set all delta values = 0 except for dz2 and dz1. Set these last two to the same value with opposite sign to model a lens out of position

A1A = Simplify[A1/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dz2-> 0, dz1-> -0}];
NominalMagnification= 1/A1A/. X1-> 0;
Print["Nominal paraxial magnification is ", NominalMagnification, "m"];
RateVariationInMag = (D[A1A, X1])/. X1-> 0;
Print["Variation in magnification for +/- 5 mm displacement of HWS = ", RateVariationInMag*0.005, "x"];

Hence the paraxial magnification is SUPER invariant to the position of F1

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMXtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> 0, dz1-> -0})];
BeamWAtHWSNom = pOutNom[[2]];
MOut = (ITMXtoHWSVar/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dz2-> 0, dz1-> -0});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamWAtHWSX1= pOutX1[[2]];
MagNom = BeamSizeAtITM/BeamWAtHWSNom;
MagX1 = BeamSizeAtITM/BeamWAtHWSX1;
Print["Absolute Variation in magnification for +5mm displacement = ", MagX1-MagNom, "x"];
Print["Relative Variation in magnification for +5mm displacement = ", ((MagX1-MagNom)/MagNom)*100, "%"];
ATestF2HWSDistance = {Round[((MagX1-MagNom)/MagNom)*100*10]/10.0,"%",  5, "mm"};

Hence the Gaussian magnification is MOSTLY invariant to the position of HWS

Variation in collimation (C = 0) with length changes: (assume that FL variations can be compensated for)  - NO ISSUES - Easy to work with!

M1 =ITMXtoHWSVar;
C1 = M1[[2]][[1]];

Test variation in position of F1 - determine the defocus - INSENSITIVE

Set all delta values = 0 except for dz2 and dz1. Set these last two to the same value with opposite sign to model a lens out of position

C1A = Simplify[C1/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> X1, dz1-> -X1}];
NominalDefocus= C1A/. X1-> 0;
Print["Nominal paraxial defocus is ", NominalDefocus, " m^-1"];
RateVariationInDefocus = (D[C1A, X1])/. X1-> 0;
Print["Variation in defocus for + 5 mm displacement of F1 = ", RateVariationInDefocus*0.005, " m^-1, or ", 1/(RateVariationInDefocus*0.005), "m"];

Hence the relative change in the paraxial defocus is quite large but the absolute change (with displacement of F1) is quite small, about 7km lens for 5mm displacement. At HWS this is not an issue.

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMXtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> 0, dz1-> -0})];
BeamRAtHWSNom = pOutNom[[1]];
MOut = (ITMXtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> X1, dz1-> -X1});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamRAtHWSX1= pOutX1[[1]];
Print["Variation in defocus of Gaussian Beam for +5mm displacement = ", 1/BeamRAtHWSX1-1/BeamRAtHWSNom, " m^-1", "or lens of = ", 1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), " m"];
CTestF1OutOfPosition = {1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), "m", 5 ,"mm"};

Hence the Gaussian curvature at HWS does not experience large changes in the radius of curvature that with motion of F1

Test variation in position of F2 - determine the defocus - INSENSITIVE

Set all delta values = 0 except for dz2 and dz3. Set these last two to the same value with opposite sign to model a lens out of position

C1A = Simplify[C1/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dz2-> -X1, dz1-> 0}];
NominalDefocus= C1A/. X1-> 0;
Print["Nominal paraxial defocus is ", NominalDefocus, " m^-1"];
RateVariationInDefocus = (D[C1A, X1])/. X1-> 0;
Print["Variation in defocus for + 5 mm displacement of F2 = ", RateVariationInDefocus*0.005, " m^-1, or ", 1/(RateVariationInDefocus*0.005), "m"];

Hence the relative change in the paraxial defocus is quite large but the absolute change (with displacement of F2) is quite small, about 8km lens for 5mm displacement. At HWS this is not an issue.

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMXtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> 0, dz1-> -0})];
BeamRAtHWSNom = pOutNom[[1]];
MOut = (ITMXtoHWSVar/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dz2-> -X1, dz1-> 0});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamRAtHWSX1= pOutX1[[1]];
Print["Variation in defocus of Gaussian Beam for +5mm displacement = ", 1/BeamRAtHWSX1-1/BeamRAtHWSNom, " m^-1", "or lens of = ", 1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), " m"];
CTestF2OutOfPosition = {1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), "m", 5 ,"mm"};

Hence the Gaussian curvature at HWS does not experience large changes in the radius of curvature that with motion of F2

Test variation in F1-F2 distance: determine the defocus - INSENSITIVE

C1A = Simplify[C1/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> -X1, dz1-> 0}];
NominalDefocus= C1A/. X1-> 0;
Print["Nominal paraxial defocus is ", NominalDefocus, " m^-1"];
RateVariationInDefocus = (D[C1A, X1])/. X1-> 0;
Print["Variation in defocus for + 5 mm displacement of F2 = ", RateVariationInDefocus*0.005, " m^-1, or ", 1/(RateVariationInDefocus*0.005), "m"];

Hence the relative change in the paraxial defocus is quite large but the absolute change (with variation in F1-F2 distance) is quite small, about 8km lens for 5mm displacement. At HWS this is not an issue.

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMXtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> 0, dz1-> -0})];
BeamRAtHWSNom = pOutNom[[1]];
MOut = (ITMXtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> -X1, dz1-> 0});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamRAtHWSX1= pOutX1[[1]];
Print["Variation in defocus of Gaussian Beam for +5mm displacement = ", 1/BeamRAtHWSX1-1/BeamRAtHWSNom, " m^-1", "or lens of = ", 1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), " m"];
CTestF1F2Distance = {1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), "m", 5 ,"mm"};

Hence the Gaussian curvature at HWS does not experience large changes in the radius of curvature that with variation in F1-F2 distance

Test variation in position of HWS - determine the defocus - INSENSITIVE

Set all delta values = 0 except for dz3. Set these last two to the same value with opposite sign to model a lens out of position

C1A = Simplify[C1/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dz2-> 0, dz1-> -0}];
NominalDefocus= C1A/. X1-> 0;
Print["Nominal paraxial defocus is ", NominalDefocus, " m^-1"];
RateVariationInDefocus = (D[C1A, X1])/. X1-> 0;
Print["Variation in defocus for + 5 mm displacement of F1 = ", RateVariationInDefocus*0.005, " m^-1, or ", 1/(RateVariationInDefocus*0.005), "m"];

Hence the relative change in the paraxial defocus is ZERO

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMXtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dz2-> 0, dz1-> -0})];
BeamRAtHWSNom = pOutNom[[1]];
MOut = (ITMXtoHWSVar/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dz2-> 0, dz1-> -0});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamRAtHWSX1= pOutX1[[1]];
Print["Variation in defocus of Gaussian Beam for +5mm displacement = ", 1/BeamRAtHWSX1-1/BeamRAtHWSNom, " m^-1", "or lens of = ", 1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), " m"];
CTestF2HWSDistance = {1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), "m", 5 ,"mm"};

Hence the Gaussian curvature at HWS does not experience large changes in the radius of curvature that with variation in HWS position - expected since beam is approximately collimated here

Difficulty in mode-matching (MMOut = MMIn) -  Take Care with Mode Matching

InputModeFromHWSWR  = InputModeFromHWS;
InputModeFromHWSRW  = {InputModeFromHWS[[2]],InputModeFromHWS[[1]]}
OutputMode = OutRW[InputModeFromHWSRW[[2]], InputModeFromHWSRW[[1]], λIN,MHWSv]

10% variation in mode size, same curvature - no worries

InputModeFromHWSRW  = {InputModeFromHWS[[2]],InputModeFromHWS[[1]]*1.1}
OutputMode = OutRW[InputModeFromHWSRW[[2]], InputModeFromHWSRW[[1]], λIN,MHWSv]

1.0m lens added to input beam, same mode size

InputModeFromHWSRW  = {((1/1.0) +InputModeFromHWS[[2]]^-1)^-1,InputModeFromHWS[[1]]}
OutputMode = OutRW[InputModeFromHWSRW[[2]], InputModeFromHWSRW[[1]], λIN,MHWSv]
Print["Difference in output curvature = ",((1/(OutputMode[[1]]))+(1/(InputModeFromHWSRW[[1]])))^-1] 
Print["Relative variation in output size with a 1m lens added to input = ",((OutputMode[[2]] - InputModeFromHWSRW[[2]])/(InputModeFromHWSRW[[2]]))*100, "%"]

Temperature change required to scatter 0.01% from HWS error

Simple model first of temperature sensitivity - assuming common Length and focal length changes
requirements: 
1. consider thermal changes in length and defocus
2. ABCD and Gaussian analysis

Look at all the telescopes with defocus and infer maximum acceptable temp. change

Input Mode is ...

InputModeFromHWS

Length sensitive matrices - main telescope between ITMXPO and SR3

ITMToITMXPOdX =Simplify[ (({{1, ITMXPOtoSR3+dx}, {0, 1}})).(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR3toITMX}, {0, 1}}))];ITMXPOToITMdX = Simplify[ (({{1, SR3toITMX}, {0, 1}})).(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, ITMXPOtoSR3+dx}, {0, 1}}))];

ITMXPOToHWSdx = Simplify[(({{1, z3+dx3}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z2+dx2}, {0, 1}})).(({{1, 0}, {-1/f1, 1}})).(({{1, z1 +dx1}, {0, 1}}))]/. TestSolutionSubstitution;
HWSToITMXPOdx = Simplify[(({{1, z1+dx1}, {0, 1}})).(({{1, 0}, {-1/f1, 1}})).(({{1, z2+dx2}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z3+dx3}, {0, 1}}))]/. TestSolutionSubstitution;

HWStoHWSdx = ITMXPOToHWSdx.ITMToITMXPOdX.MITM.ITMXPOToITMdX.HWSToITMXPOdx/. {dx1-> 0, dx2-> 0, dx3-> 0};

Check the matrix

HWStoHWS = HWStoHWSdx/. dx-> 0;
InputModeFromHWS
pOut = OutRW[InputModeFromHWS[[1]], InputModeFromHWS[[2]], λIN,HWStoHWS];
Reverse[pOut]

A1: Get the defocus at ITM for a given displacement between HAM4 and HAM5

To scatter out 0.01% the maximum acceptable defocus is 1.3*10^-6 m^-1;

MaxAcceptableDefocusAtITM = MaxDefocusAtITMAll;

dxTest = 0.00025;
dtM = HWStoHWSdx/. dx-> dxTest;
InputModeFromHWS;
Print["Input defocus = ", 1/InputModeFromHWS[[2]], " m^-1"]
pOut = OutRW[InputModeFromHWS[[1]], InputModeFromHWS[[2]], λIN,dtM];
Print["Output defocus = ", 1/pOut[[1]], " m^-1"]
dS = 1/(-InputModeFromHWS[[2]]) - 1/pOut[[1]]; 
M2 = BeamSizeAtITM/InputModeFromHWS[[1]];

dSAtITM = dS/(M2^2);
Print["Defocus at ITM = ", dSAtITM, "m^-1 for ", Round[dxTest*10^6], "um displacement"]
dSdxITMXPOtoSR3 = dSAtITM/dxTest;
Print["Max Acceptable defocus at ITM = ",MaxAcceptableDefocusAtITM," m^-1"]

A2: Limit apparent temperature change of the concrete between IMTXPO & SR3

Determine how much the previous dx would need to be scaled by to yield the maximum apparent defocus allowed at the ITM

dxScale = MaxAcceptableDefocusAtITM/dSAtITM

CoefficientThermalExpansionOfConcrete = 1.2*10^-5;
dT = dx/(L*CoefficientThermalExpansionOfConcrete)/. {L-> ITMXPOtoSR3, dx->dxTest*dxScale};
dTITMXPOtoSR3 = {Round[dT*100]/100.0, "K"};
Print["Maximum acceptable change in temperature of concrete slab between HAM4 and HAM5 = ", Round[dT*100]*10, " mK"]

Now we need a secondary set of optics that introduce the same level of thermal defocus and give a small beam at the reference retro-reflector.

B1: Get the defocus at ITM for a given displacement between ITMXPO and F1

Print["Magnification between HWS and ITM = ", M2]

HWStoHWSdx1 = ITMXPOToHWSdx.ITMToITMXPOdX.MITM.ITMXPOToITMdX.HWSToITMXPOdx/. {dx-> 0, dx2-> 0, dx3-> 0};

dxTest = 0.00025;
dtM = HWStoHWSdx1/. dx1-> dxTest;
InputModeFromHWS;
Print["Input defocus = ", 1/InputModeFromHWS[[2]], " m^-1"]
pOut = OutRW[InputModeFromHWS[[1]], InputModeFromHWS[[2]], λIN,dtM];
Print["Output defocus = ", 1/pOut[[1]], " m^-1"]
dS = 1/(-InputModeFromHWS[[2]]) - 1/pOut[[1]]; 
M2 = BeamSizeAtITM/InputModeFromHWS[[1]];
dSAtITM = dS/(M2^2);
Print["Defocus at ITM = ", dSAtITM, "m^-1 for ", Round[dxTest*10^6], "um displacement"]
dSdx = dSAtITM/dxTest;
Print["Max Acceptable defocus at ITM = ",MaxAcceptableDefocusAtITM," m^-1"]

B2: Limit apparent temperature change of the stainless between ITMXPO & F1

Determine how much the previous dx would need to be scaled by to yield the maximum apparent defocus allowed at the ITM

dxScale = MaxAcceptableDefocusAtITM/dSAtITM

CoefficientThermalExpansionOfStainlessSteel = 1.73*10^-5;
dT = dx/(L*CoefficientThermalExpansionOfStainlessSteel)/. {L-> Z1A, dx->dxTest*dxScale};
dTITMXPOtoF1 = {Round[dT*100]/100.0, "K"};
Print["Maximum acceptable change in temperature of HAM4 for beam between F1 and ITMXPO = ", Round[dT*100]*10, " mK"]

C1: Get the defocus at ITM for a given displacement between F1 and F2

HWStoHWSdx2 = ITMXPOToHWSdx.ITMToITMXPOdX.MITM.ITMXPOToITMdX.HWSToITMXPOdx/. {dx-> 0, dx1-> 0, dx3-> 0};

dxTest = 0.00025;
dtM = HWStoHWSdx2/. dx2-> dxTest;
InputModeFromHWS;
Print["Input defocus = ", 1/InputModeFromHWS[[2]], " m^-1"]
pOut = OutRW[InputModeFromHWS[[1]], InputModeFromHWS[[2]], λIN,dtM];
Print["Output defocus = ", 1/pOut[[1]], " m^-1"]
dS = 1/(-InputModeFromHWS[[2]]) - 1/pOut[[1]]; 
M2 = BeamSizeAtITM/InputModeFromHWS[[1]];
dSAtITM = dS/(M2^2);
Print["Defocus at ITM = ", dSAtITM, "m^-1 for ", Round[dxTest*10^6], "um displacement"]
dSdx = dSAtITM/dxTest;
Print["Max Acceptable defocus at ITM = ",MaxAcceptableDefocusAtITM," m^-1"]

C2: Limit apparent temperature change of the stainless between F1 & F2

Determine how much the previous dx would need to be scaled by to yield the maximum apparent defocus allowed at the ITM

dxScale = MaxAcceptableDefocusAtITM/dSAtITM

CoefficientThermalExpansionOfStainlessSteel = 1.73*10^-5;
dT = dx/(L*CoefficientThermalExpansionOfStainlessSteel)/. {L-> Z2A, dx->dxTest*dxScale};
dTF1toF2 = {Round[dT*100]/100.0, "K"};
Print["Maximum acceptable change in temperature of HAM4 and ISCT to affect F1 and F2 distance = ", Round[dT*100]*10, " mK"]

D1: Get the defocus at ITM for a given displacement between F2 and HWS

HWStoHWSdx3 = ITMXPOToHWSdx.ITMToITMXPOdX.MITM.ITMXPOToITMdX.HWSToITMXPOdx/. {dx-> 0, dx1-> 0, dx2-> 0};

dxTest = 0.00025;
dtM = HWStoHWSdx3/. dx3-> dxTest;
InputModeFromHWS;
Print["Input defocus = ", 1/InputModeFromHWS[[2]], " m^-1"]
pOut = OutRW[InputModeFromHWS[[1]], InputModeFromHWS[[2]], λIN,dtM];
Print["Output defocus = ", 1/pOut[[1]], " m^-1"]
dS = 1/(-InputModeFromHWS[[2]]) - 1/pOut[[1]]; 
M2 = BeamSizeAtITM/InputModeFromHWS[[1]];
dSAtITM = dS/(M2^2);
Print["Defocus at ITM = ", dSAtITM, "m^-1 for ", Round[dxTest*10^6], "um displacement"]
dSdx = dSAtITM/dxTest;
Print["Max Acceptable defocus at ITM = ",MaxAcceptableDefocusAtITM," m^-1"]

D2: Limit apparent temperature change of the stainless between F2 & HWS

Determine how much the previous dx would need to be scaled by to yield the maximum apparent defocus allowed at the ITM

dxScale = MaxAcceptableDefocusAtITM/dSAtITM

CoefficientThermalExpansionOfStainlessSteel = 1.73*10^-5;
dT = dx/(L*CoefficientThermalExpansionOfStainlessSteel)/. {L-> Z3A, dx->dxTest*dxScale};
dTF2toHWS= {Round[dT*100]/100.0, "K"};
Print["Maximum acceptable change in temperature of ISCT table with HWS = ", Round[dT*100]/100., " K"]

Look at all the telescopes with defocus of lenses and infer maximum acceptable temp. change

Input Mode is ...

InputModeFromHWS

Length sensitive matrices - main telescope between ITMXPO and SR3

ITMToITMXPOdX =Simplify[ (({{1, ITMXPOtoSR3}, {0, 1}})).(({{1, 0}, {((-2)/(SR3RoC(1+ alphaFS * dTSR3))), 1}})).(({{1, SR3toITMX}, {0, 1}}))];ITMXPOToITMdX = Simplify[ (({{1, SR3toITMX}, {0, 1}})).(({{1, 0}, {((-2)/(SR3RoC(1+ alphaFS * dTSR3))), 1}})).(({{1, ITMXPOtoSR3}, {0, 1}}))];

ITMXPOToHWSdx = Simplify[(({{1, z3}, {0, 1}})).(({{1, 0}, {-1/(f2(1+ alphaFS * dTF2)), 1}})).(({{1, z2}, {0, 1}})).(({{1, 0}, {-1/(f1(1+ alphaFS * dTF1)), 1}})).(({{1, z1 }, {0, 1}}))]/. TestSolutionSubstitution;
HWSToITMXPOdx = Simplify[(({{1, z1}, {0, 1}})).(({{1, 0}, {-1/(f1(1+ alphaFS * dTF1)), 1}})).(({{1, z2}, {0, 1}})).(({{1, 0}, {-1/(f2(1+ alphaFS * dTF2)), 1}})).(({{1, z3}, {0, 1}}))]/. TestSolutionSubstitution;

A1: Get the defocus at ITM for a uniform temperature change of all the optics

To scatter out 0.01% the maximum acceptable defocus is 1.3*10^-6 m^-1;

Fused Silica

Assume that the optics are all fused silica

HWStoHWSdT1 = ITMXPOToHWSdx.ITMToITMXPOdX.MITM.ITMXPOToITMdX.HWSToITMXPOdx/. {alphaFS-> 0.55*10^-6, dTSR3-> dT1,  dTF1-> dT1, dTF2-> dT1};

dT1Test = 10.0;
dtM = HWStoHWSdT1/. dT1-> dT1Test;
InputModeFromHWS;
Print["Input defocus = ", 1/InputModeFromHWS[[2]], " m^-1"]
pOut = OutRW[InputModeFromHWS[[1]], InputModeFromHWS[[2]], λIN,dtM];
Print["Output defocus = ", 1/pOut[[1]], " m^-1"]
dS = 1/(-InputModeFromHWS[[2]]) - 1/pOut[[1]]; 
M2 = BeamSizeAtITM/InputModeFromHWS[[1]];

dSAtITM = dS/(M2^2);
Print["Defocus at ITM = ", dSAtITM, "m^-1 for ", Round[dT1Test], "K temperature change"]
dSdxITMXPOtoSR3 = dSAtITM/dxTest;
Print["Max Acceptable defocus at ITM = ",MaxAcceptableDefocusAtITM," m^-1"]

BK7

Assume that the optics are all BK7, except SR3 which is assumed to stay the same temperature

HWStoHWSdT1 = ITMXPOToHWSdx.ITMToITMXPOdX.MITM.ITMXPOToITMdX.HWSToITMXPOdx/. {alphaFS-> 7.5*10^-6, dTSR3-> 0,  dTF1-> dT1, dTF2-> dT1};

MaxAcceptableDefocusAtITM = MaxDefocusAtITMAll;

dT1Test = 10.0;
dtM = HWStoHWSdT1/. dT1-> dT1Test;
InputModeFromHWS;
Print["Input defocus = ", 1/InputModeFromHWS[[2]], " m^-1"]
pOut = OutRW[InputModeFromHWS[[1]], InputModeFromHWS[[2]], λIN,dtM];
Print["Output defocus = ", 1/pOut[[1]], " m^-1"]
dS = 1/(-InputModeFromHWS[[2]]) - 1/pOut[[1]]; 
M2 = BeamSizeAtITM/InputModeFromHWS[[1]];

dSAtITM = dS/(M2^2);
Print["Defocus at ITM = ", dSAtITM, "m^-1 for ", Round[dT1Test], "K temperature change"]
dSdxITMXPOtoSR3 = dSAtITM/dxTest;
Print["Max Acceptable defocus at ITM = ",MaxAcceptableDefocusAtITM," m^-1"]

Report on Solution

Print["=== REPORT ON H1:ITMX-HWS SOLUTION ==="]
Print[" "];
Print["ITMX apparent Radius of Curvature = ", ITMRoC, "m"];
Print["ITM to SR3 optical distance = ", SR3toITMX, "m"];
Print["BS to SR3 optical distance = ", SR3toBS, "m"];
Print["SR3 Radius of Curvature = ", SR3RoC, "m"];
Print["SR3 to ITMXPO distance = ", ITMXPOtoSR3, "m"];
Print[" "];
Print["ITMXPO to F1 distance = ", Round[(Z1A)*1000]/1000.0, "m"];
Print["F1 focal length = ", Round[(F1A)*1000]/1000.0, "m"];
Print["F1 to F2 distance = ", Round[(Z2A )*1000]/1000.0, "m"];
Print["F2 focal length = ", Round[(F2A)*1000]/1000.0, "m"];
Print["F2 to HWS distance = ", Round[(Z3A )*1000]/1000.0, "m"];
Print[" "];
Print["----"]
Print["TEMPERATURE SENSITIVITY"]
Print["Maximum Acceptable defocus at the ITM = ", MaxAcceptableDefocusAtITM];
Print["ITMXPO to SR3 - maximum acceptable temperature change = ", dTITMXPOtoSR3[[1]], dTITMXPOtoSR3[[2]]]
Print["F1 to ITMXPO - maximum acceptable temperature change = ", dTITMXPOtoF1[[1]], dTITMXPOtoF1[[2]]]
Print["F1 to F2 - maximum acceptable temperature change = ", dTF1toF2[[1]], dTF1toF2[[2]]]
Print["F2 to HWS - maximum acceptable temperature change = ", dTF2toHWS[[1]], dTF2toHWS[[2]]]
Print[" "];
Print["----"]
Print["TOLERANCING"]
Print["Imaging (B=0): F1 out of Position by ", BTestF1OutOfPosition[[3]], BTestF1OutOfPosition[[4]], " moves conjugate plane by ", BTestF1OutOfPosition[[1]], BTestF1OutOfPosition[[2]]]
Print["Imaging (B=0): F2 out of Position by ", BTestF2OutOfPosition[[3]], BTestF2OutOfPosition[[4]], " moves conjugate plane by ", BTestF2OutOfPosition[[1]], BTestF2OutOfPosition[[2]]]
Print["Imaging (B=0): F1-F2 distance varies by ", BTestF1F2Distance[[3]], BTestF1F2Distance[[4]], " moves conjugate plane by ", BTestF1F2Distance[[1]], BTestF1F2Distance[[2]]]
Print["Imaging (B=0): HWS out of Position by ", BTestF2HWSDistance[[3]], BTestF2HWSDistance[[4]], " moves conjugate plane by ", BTestF2HWSDistance[[1]], BTestF2HWSDistance[[2]]]
Print[" "]
Print["Magnification (A=1/17.5): F1 out of Position by ", ATestF1OutOfPosition[[3]], ATestF1OutOfPosition[[4]], " varies magnification by ", ATestF1OutOfPosition[[1]], ATestF1OutOfPosition[[2]]]
Print["Magnification (A=1/17.5): F2 out of Position by ", ATestF2OutOfPosition[[3]], ATestF2OutOfPosition[[4]], " varies magnification by ", ATestF2OutOfPosition[[1]], ATestF2OutOfPosition[[2]]]
Print["Magnification (A=1/17.5): F1-F2 distance varies by ", ATestF1F2Distance[[3]], ATestF1F2Distance[[4]], " varies magnification by ", ATestF1F2Distance[[1]], ATestF1F2Distance[[2]]]
Print["Magnification (A=1/17.5): F2-HWS distance varies by ", ATestF2HWSDistance[[3]], ATestF2HWSDistance[[4]], " varies magnification by ", ATestF2HWSDistance[[1]], ATestF2HWSDistance[[2]]]
Print[" "]
Print["Collimation (C=0): F1 out of Position by ", CTestF1OutOfPosition[[3]], CTestF1OutOfPosition[[4]], " looks like an additional lens at HWS with f = ", CTestF1OutOfPosition[[1]], CTestF1OutOfPosition[[2]]]
Print["Collimation (C=0): F2 out of Position by ", CTestF2OutOfPosition[[3]], CTestF2OutOfPosition[[4]], " looks like an additional lens at HWS with f = ", CTestF2OutOfPosition[[1]], CTestF2OutOfPosition[[2]]]
Print["Collimation (C=0): F1-F2 distance varied by ", CTestF1F2Distance[[3]], CTestF1F2Distance[[4]], " looks like an additional lens at HWS with f = ", CTestF1F2Distance[[1]], CTestF1F2Distance[[2]]]
Print["Collimation (C=0): F2-HWS distance varied by ", CTestF2HWSDistance[[3]], CTestF2HWSDistance[[4]], " looks like an additional lens at HWS with f = ", CTestF2HWSDistance[[1]], CTestF2HWSDistance[[2]]]

Show[ITMXPlottedSolutionFromITMX]
Show[ITMXPlottedSolutionFromITMXPO]

Show[ITMXPlottedSolutionFromITMXPO, PlotRange-> {{0, 11}, {0, 14}}]

Beam size at the mirrors

Create a Module to determine the beam size

GetTheBeamSize[ABCDList_, RWInputMode_, λIN_, zIn_]:=Module[{},
z0qq = 0;
plotList = {};
count1 = 0;
zOut = 0.0;
For[ii=1, ii≤ Length[ABCDList], ii++, {
M1IN = ABCDList[[ii]];

plotAll = Plot[0*x, {x, 0, 0.00001}];
If[M1IN[[1]][[2]]≠ 0, {
MOut = (({{1, 0}, {0, 1}}));
For[jj=1, jj≤ ii-1, jj++, {
MOut = ABCDList[[jj]].MOut;
}];
MOut = (({{1, zIn-z0qq}, {0, 1}})).MOut;

eqOut = OutRW[RWInputMode[[2]], RWInputMode[[1]], λIN, MOut];
If[zIn≥  z0qq, zOut = eqOut[[2]]];



count1 = count1 + 1;

z0qq = ABCDList[[ii]][[1]][[2]]+z0qq;



}];
}];
Return[zOut];
];
GetTheGouyPhase[ABCDList_, RWInputMode_, λIN_, zIn_]:=Module[{x0, waistRWqq},
z0qq = 0;
plotList = {};
count1 = 0;
zOut = 0.0;
For[ii=1, ii≤ Length[ABCDList], ii++, {
M1IN = ABCDList[[ii]];

plotAll = Plot[0*x, {x, 0, 0.00001}];
If[M1IN[[1]][[2]]≠ 0, {
MOut = (({{1, 0}, {0, 1}}));
For[jj=1, jj≤ ii-1, jj++, {
MOut = ABCDList[[jj]].MOut;
}];
MOut = (({{1, zIn-z0qq}, {0, 1}})).MOut;

eqOut = OutRW[RWInputMode[[2]], RWInputMode[[1]], λIN, MOut];
If[zIn≥  z0qq, zOut = eqOut];



count1 = count1 + 1;

z0qq = ABCDList[[ii]][[1]][[2]]+z0qq;



}];
}];

x0 = FindNearestWaist[zOut[[2]], zOut[[1]], λIN];
waistRWqq = OutRW[zOut[[2]], zOut[[1]], λIN,(({{1, x0}, {0, 1}}))];
GouyPhaseqq = ArcTan[((-x0)/((((Pi waistRWqq[[2]]^2)/(λIN)))))];
Return[GouyPhaseqq];
];

Import the data on the mirror positions

dataIN = Import["/Users/aidanbrooks/Documents/LIGO/documents/Advanced_LIGO/TCS/Hartmann sensor/modeling/H1HWSX_coords.txt", "TSV"];
totaldistance = 0;
distList = {};
nameList = {};
For[ii = 1, ii≤ Length[dataIN], ii++, {
d1 = dataIN[[ii]][[7]]/1000.0;
totaldistance = totaldistance + d1;
dataIN[[ii]] = Append[dataIN[[ii]], totaldistance];
distList = Append[distList , totaldistance];
nameList = Append[nameList, dataIN[[ii]][[1]]];
}];

dataIN

Distance List from T1100463-v5

distances = {0,972.7,1650.0,2663.7,2926.6,4883.2,5563.0,5811.4,5861.4,6009.9,6959.7,7760.0,7850.0,8278.8,8328.8,8757.6,9172.6,9922.6,10084.4,10234.4,10634.4,10721.4}/1000;
distList = distances;

Determine the beam size at each interface

beamSizeList = {};
For[kk=1, kk≤Length[distList], kk++, {
beamsize =GetTheBeamSize[{(({{1, z1 }, {0, 1}})),(({{1, 0}, {-1/f1, 1}})), (({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z3}, {0, 1}}))}/. TestSolutionSubstitution, OutModeITMXPO, 840*10^-9, distList[[kk]]];

beamSizeList = Append[beamSizeList, {distList[[kk]],beamsize}];

If[kk==3, Print[distList[[kk]]]];
If[kk==7,{
Print[distList[[kk]]-0.877];
wVP = GetTheBeamSize[{(({{1, z1 }, {0, 1}})),(({{1, 0}, {-1/f1, 1}})), (({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z3}, {0, 1}}))}/. TestSolutionSubstitution, OutModeITMXPO, 840*10^-9, distList[[kk]]-0.877];
 wUPM = GetTheBeamSize[{(({{1, z1 }, {0, 1}})),(({{1, 0}, {-1/f1, 1}})), (({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z3}, {0, 1}}))}/. TestSolutionSubstitution, OutModeITMXPO, 840*10^-9, distList[[kk]]];
}];
}];

MtoVP =(({{1, 4.0533 - 1.64901}, {0, 1}})). (({{1, 0}, {-1/f1, 1}})).(({{1, z1 }, {0, 1}})).ITMToITMXPO/. TestSolutionSubstitution;
w1 = Sqrt[2]*54.0*10^-3;

div = ((wVP-wUPM)/0.877);
w0 = ((λIN)/(Pi div))*1000;
f0 = wVP/div;

beamSizeList

ListPlot[beamSizeList]

Export["/Users/aidanbrooks/Documents/LIGO/documents/Advanced_LIGO/TCS/Hartmann sensor/design/optical layout and modeling/H1HWSX_coords_beamsize.txt", beamSizeList,"Table"];

MatrixForm[beamSizeList]

For[ii = 1, ii≤ Length[beamSizeList], ii++, Print[1000*beamSizeList[[ii]][[2]]]]

Determine the beam size at each interface

GouyPhase =GetTheGouyPhase[{(({{1, z1 }, {0, 1}})),(({{1, 0}, {-1/f1, 1}})), (({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z3}, {0, 1}}))}/. TestSolutionSubstitution, OutModeITMXPO, 840*10^-9, distList[[1]]]

beamSizeList = {};
GouyPhaseList = {};
GouyPhaseListDeg = {};
For[kk=1, kk≤Length[distList], kk++, {
beamsize =GetTheBeamSize[{(({{1, z1 }, {0, 1}})),(({{1, 0}, {-1/f1, 1}})), (({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z3}, {0, 1}}))}/. TestSolutionSubstitution, OutModeITMXPO, 840*10^-9, distList[[kk]]];

beamSizeList = Append[beamSizeList, {distList[[kk]],beamsize}];

GouyPhase =GetTheGouyPhase[{(({{1, z1 }, {0, 1}})),(({{1, 0}, {-1/f1, 1}})), (({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z3}, {0, 1}}))}/. TestSolutionSubstitution, OutModeITMXPO, 840*10^-9, distList[[kk]]];

GouyPhaseList = Append[GouyPhaseList, {distList[[kk]],GouyPhase}];
GouyPhaseListDeg = Append[GouyPhaseListDeg, {distList[[kk]],GouyPhase*180/Pi}];

}];

plot1 = ListPlot[GouyPhaseListDeg, Frame-> True, GridLines-> {Automatic, Table[ii, {ii, -90, 90, 15}]}];
plot2 = ListPlot[GouyPhaseListDeg, Frame-> True, Joined-> True];
Show[plot1, plot2]

GouyPhaseOptics = {};
For[jj = 1, jj≤ Length[GouyPhaseListDeg], jj++, {
GouyPhaseOptics = Append[GouyPhaseOptics, {nameList[[jj]], GouyPhaseListDeg[[jj]][[2]]}];
}];
MatrixForm[GouyPhaseOptics]

For[kk=1, kk≤ Length[distList], kk++, {
dataIN[[kk]] = Append[dataIN[[kk]], beamSizeList[[kk]][[2]]*1000.0];
}];

ListPlot[beamSizeList]

Export["/Users/aidanbrooks/Documents/LIGO/documents/Advanced_LIGO/TCS/Hartmann sensor/modeling/H1HWSX_coords_beamsize.txt", dataIN,"Table"];

MatrixForm[dataIN]

Astigmatism

Theta = 6.5;

AstigmaticRadiiCurvature = {(f1*2)*Cos[Theta*(Pi/180)]^2,(f2*2)}/. TestSolutionSubstitution
AstigmaticFocalLengthsMM = AstigmaticRadiiCurvature/2.0

Horizontal Direction - work out change in conj. plane position

ITMToITMXPO;
ImgMatrix = Simplify[(({{1, z3}, {0, 1}})).(({{1, 0}, {-1/AstigmaticFocalLengthsMM[[2]], 1}})).(({{1, z2}, {0, 1}})).(({{1, 0}, {-1/(AstigmaticFocalLengthsMM[[1]]), 1}})).(({{1, z1}, {0, 1}}))];
SR2ARToHWS = ImgMatrix/. TestSolutionSubstitution;
M1 =(({{1, z}, {0, 1}})). SR2ARToHWS.ITMToITMXPO;
zsoln = Solve[M1[[1]][[2]]==0, z]
mag = M1[[2]][[2]]/. zsoln[[1]]
imgHsoln = {z, mag}/. zsoln[[1]]
zSolve = z/.zsoln[[1]]

Horizontal Direction - compensate with second lens

ITMToITMXPO;
ImgMatrix = Simplify[(({{1, z3}, {0, 1}})).(({{1, 0}, {-1/(f2*a1), 1}})).(({{1, z2}, {0, 1}})).(({{1, 0}, {-1/(AstigmaticFocalLengthsMM[[1]]), 1}})).(({{1, z1}, {0, 1}}))];
SR2ARToHWS = ImgMatrix/. TestSolutionSubstitution;
M1 = SR2ARToHWS.ITMToITMXPO;
asoln = Solve[M1[[1]][[2]]==0, a1]
AngleV = (180/Pi)ArcCos[Sqrt[a1]]/. asoln[[1]];
Print["Angle of Incidence for Lens 2 = ", ToString[AngleV], " degrees"];
M1Vertical =Chop[ M1/.asoln[[1]]];
MatrixForm[M1Vertical]
(f2/. TestSolutionSubstitution)*2000*a1/. asoln[[1]]

Major Optics Misalignments - ABCDEFGH formalism

TestSolutionSubstitution

solnMis = {VacLens->-1.699,AirLens->1.699,POtoVacLens->2.156,VacLenstoAirLens->2.808,AirLenstoHWS->2.706}

MITMmis =  (({{1, 0, 0}, {((-2)/(-ITMRoC)), 1, 0}, {0, 0, 1}}));
MAirLensToVacLens = (({{1, VacLenstoAirLens, 0}, {0, 1, 0}, {0, 0, 1}}))
MVacLensToAirLens = (({{1, VacLenstoAirLens, 0}, {0, 1, 0}, {0, 0, 1}}))

MMistoITM =(({{1, 0, 0}, {0, 1, 2*ITMXdtheta}, {0, 0, 1}})).(({{1, 0.2*SR3toITMX, 0}, {0, 1, 0}, {0, 0, 1}})).(({{1, 0, 0}, {0, 1, 2*BSdtheta}, {0, 0, 1}})).(({{1, 0.8*SR3toITMX, 0}, {0, 1, 0}, {0, 0, 1}})).(({{1, 0, 0}, {-2/SR3RoC, 1, 0}, {0, 0, 1}})).(({{1, 0, 0}, {0, 1, 2*SR3dtheta}, {0, 0, 1}})).(({{1, ITMXPOtoSR3, 0}, {0, 1, 0}, {0, 0, 1}})).(({{1, POtoVacLens, 0}, {0, 1, 0}, {0, 0, 1}})).(({{1, 0, 0}, {-1/VacLens, 1, 0}, {0, 0, 1}})).MAirLensToVacLens.(({{1, 0, 0}, {-1/AirLens, 1, 0}, {0, 0, 1}})).(({{1, AirLenstoHWS, 0}, {0, 1, 0}, {0, 0, 1}}));
MMisfromITM =(({{1, AirLenstoHWS, 0}, {0, 1, 0}, {0, 0, 1}})).(({{1, 0, 0}, {-1/AirLens, 1, 0}, {0, 0, 1}})).MVacLensToAirLens.(({{1, 0, 0}, {-1/VacLens, 1, 0}, {0, 0, 1}})).(({{1, POtoVacLens, 0}, {0, 1, 0}, {0, 0, 1}})).(({{1, ITMXPOtoSR3, 0}, {0, 1, 0}, {0, 0, 1}})).(({{1, 0, 0}, {-2/SR3RoC, 1, 0}, {0, 0, 1}})).(({{1, 0, 0}, {0, 1, 2*SR3dtheta}, {0, 0, 1}})).(({{1, 0.8*SR3toITMX, 0}, {0, 1, 0}, {0, 0, 1}})).(({{1, 0, 0}, {0, 1, 2*BSdtheta}, {0, 0, 1}})).(({{1, 0.2*SR3toITMX, 0}, {0, 1, 0}, {0, 0, 1}}));
Mrt =MMisfromITM.MITMmis.MMistoITM.{y,alpha,1}/. solnMis;
Mrt/. {alpha-> 0, ITMXdtheta->0, SR3dtheta->0, BSdtheta-> 10^-3}

17.5*2*10^-3

0.2*SR3toITMX*4*10^-3/17.5

Beam Splitter Misalignment

MBS0 = Mrt/. {alpha-> 0, ITMXdtheta->0, SR3dtheta->0, BSdtheta-> 0*10^-3};
MBS1 = Mrt/. {alpha-> 0, ITMXdtheta->0, SR3dtheta->0, BSdtheta-> 1*10^-3};

MBS1-MBS0

SR3 Misalignment

MSR30 = Mrt/. {alpha-> 0, ITMXdtheta->0, SR3dtheta->0, BSdtheta-> 0*10^-3};
MSR31 = Mrt/. {alpha-> 0*10^-3, ITMXdtheta->0*10^-3, SR3dtheta->1*10^-3, BSdtheta-> 0*10^-3};

MSR31 - MSR30

Input change Misalignment

Min0 = Mrt/. {alpha-> 0, ITMXdtheta->0, SR3dtheta->0, BSdtheta-> 0*10^-3};
Min1 = Mrt/. {alpha-> 1*10^-3, ITMXdtheta->0*10^-3, SR3dtheta->0*10^-3, BSdtheta-> 0*10^-3};

Min1 - Min0

Retro-reflected beam on the table

RRbeam = Simplify[(({{1, AirLenstoHWS}, {0, 1}})).(({{1, 0}, {-1/AirLens, 1}})).(({{1, L2}, {0, 1}})).(({{1, 0}, {1/1.7, 1}})).(({{1, 0.05}, {0, 1}})).(({{1, 0}, {1/1.7, 1}})).(({{1, L2}, {0, 1}})).(({{1, 0}, {-1/AirLens, 1}})).(({{1, AirLenstoHWS}, {0, 1}}))]/. {VacLens->-1.699,AirLens->1.699,POtoVacLens->2.156,VacLenstoAirLens->2.808,AirLenstoHWS->2.706};
RRbeam = Simplify[(({{1, L2}, {0, 1}}))]/. {VacLens->-1.699,AirLens->1.699,POtoVacLens->2.156,VacLenstoAirLens->2.808,AirLenstoHWS->2.706}
RRbeam = Simplify[(({{1, L2}, {0, 1}})).(({{1, 0}, {-1/AirLens, 1}})).(({{1, 0.05}, {0, 1}})).(({{1, 0}, {-1/AirLens, 1}})).(({{1, L2}, {0, 1}}))]/. {VacLens->-1.699,AirLens->1.699,POtoVacLens->2.156,VacLenstoAirLens->2.808,AirLenstoHWS->2.706};

RWOUT = OutRW[InputModeFromHWSRW[[2]], InputModeFromHWSRW[[1]], λIN,RRbeam];
Plot[RWOUT[[1]], {L2, 0, 2}]

RWOUT = OutRW[InputModeFromHWSRW[[2]], InputModeFromHWSRW[[1]], λIN,RRbeam];
Plot[RWOUT[[2]]*10^3, {L2, 0, 4}, PlotRange-> {0, 50}]

H1-ITMY Vertex Hartmann Sensor
Optical Layout: [AirLens = 1.1337m, VacLens = 0.79362m] @ 820nm
30-Apr-2013: This solution was determined on paper after a sign error on SR2 was discovered in -v16. The new solution is analyzed, but not derived from scratch, here. The input beam size was reduced to 127mm based on the availability of collimating lenses.

Requirements
1. Image limiting aperture
2. Magnification = +1/17.5 x
3. Reference beam with own imaging device
4. Determine temperature sensitivity
5. Tolerancing on imaging optics

Get the Zemax Data - Import Units = mm, Convert to m

HAM Tables and BSCs

ConvertToM = 1/1000;

H1fileNAME = "/Users/aidanbrooks/Documents/LIGO/documents/Advanced_LIGO/TCS/Hartmann Sensor/design/optical layout and modeling/optical_paths";
SetDirectory[H1fileNAME];
H1ZEMAX = Import["D0901920-v6_H1_layout.txt", "TSV"];

HAM4L = Position[H1ZEMAX, "HAM4 ref"];
HAM4location = Take[H1ZEMAX[[HAM4L[[1]][[1]]]]];
HAM4XYZ = Reverse[HAM4location[[Range[6, 8]]]]*ConvertToM;


HAM5L = Position[H1ZEMAX, "HAM5 ref"];
HAM5location = Take[H1ZEMAX[[HAM5L[[1]][[1]]]]];
HAM5XYZ = Reverse[HAM5location[[Range[6, 8]]]]*ConvertToM;

BSC1L = Position[H1ZEMAX, "BSC1 REF"];
BSC1location = Take[H1ZEMAX[[BSC1L[[1]][[1]]]]];
BSC1XYZ = Reverse[BSC1location[[Range[6, 8]]]]*ConvertToM;

BSC2L = Position[H1ZEMAX, "BSC2 REF"];
BSC2location = Take[H1ZEMAX[[BSC2L[[1]][[1]]]]];
BSC2XYZ = Reverse[BSC2location[[Range[6, 8]]]]*ConvertToM;

BSC3L = Position[H1ZEMAX, "BSC3 REF"];
BSC3location = Take[H1ZEMAX[[BSC3L[[1]][[1]]]]];
BSC3XYZ = Reverse[BSC3location[[Range[6, 8]]]]*ConvertToM;


L1 = 1110.0*ConvertToM;
HAMTable = Table[Reverse[{0, L1 Cos[ii*Pi/3], L1 Sin[ii*Pi/3]}], {ii, 0.5, 6.5}];
BSCTable = Table[Reverse[{0, L1 Cos[ii*Pi/30], L1 Sin[ii*Pi/30]}], {ii, 0, 60}];
BSCTable = {{-676.275,630.9613938,0},{-228.6, 958.85,0},{228.6,958.85,0},{676.275,630.9613938,0},{676.275,-630.9613938, 0},{228.6,-958.85, 0},{-228.6,-958.85, 0},{-676.275,-630.9613938,0}}*ConvertToM;

HAM4G = Polygon[HAMTable + Table[HAM4XYZ, {ii, 1, 7}]];
HAM5G = Polygon[HAMTable + Table[HAM5XYZ, {ii, 1, 7}]];
BSC1G =Translate[ Rotate[Polygon[BSCTable], 90 Degree, {0, 0 ,1}], BSC1XYZ];
BSC2G = Translate[ Rotate[Polygon[BSCTable], 90 Degree, {0, 0 ,1}], BSC2XYZ];
BSC3G = Polygon[BSCTable + Table[BSC3XYZ, {ii, 1, Length[BSCTable]}]];

IFO = Graphics3D[{HAM4G, HAM5G, BSC1G, BSC2G, BSC3G} ,ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

SR2

SR2L = Position[H1ZEMAX, "SR2_RADIUS"];
SR2location = Take[H1ZEMAX[[SR2L[[1]][[1]]]]];
SR2XYZ =Reverse[SR2location[[Range[6, 8]]]]*ConvertToM;
SR2RRT  = SR2location[[Range[15, 17]]]*ConvertToM;
SR2ang = SR2location[[Range[9, 11]]]
SR2G = Cylinder[{SR2XYZ, SR2XYZ + SR2RRT[[3]]*Reverse[{0, Sin[2Pi - SR2ang[[1]]*Pi/180], Cos[ 2Pi - SR2ang[[1]]*Pi/180]}]}, SR2RRT[[1]]]
 Graphics3D[{SR2G, HAM4G},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

ITMXPO

ITMXPOL = Position[H1ZEMAX, "ITMXPO STEER M1"];
ITMXPOlocation = Take[H1ZEMAX[[ITMXPOL[[1]][[1]]]]];
ITMXPOXYZ =Reverse[ITMXPOlocation[[Range[6, 8]]]]*ConvertToM;
ITMXPORRT  = ITMXPOlocation[[Range[13, 15]]]*ConvertToM;
ITMXPORRT = {ITMXPORRT[[1]], ITMXPORRT[[3]], ITMXPORRT[[2]]};
ITMXPOang = ITMXPOlocation[[Range[9, 11]]]
ITMXPOG = Cylinder[{ITMXPOXYZ, ITMXPOXYZ + ITMXPORRT[[3]]*Reverse[{0, Sin[2Pi - ITMXPOang[[1]]*Pi/180], Cos[ 2Pi - ITMXPOang[[1]]*Pi/180]}]}, ITMXPORRT[[1]]]
 Graphics3D[{ITMXPOG, SR2G, HAM4G},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

SR3

SR3L = Position[H1ZEMAX, "SR3 REF"];
SR3location = Take[H1ZEMAX[[SR3L[[1]][[1]]]]];
SR3XYZ =Reverse[SR3location[[Range[6, 8]]]]*ConvertToM;

SR3L = Position[H1ZEMAX, "SR3_RADIUS"];
SR3location = Take[H1ZEMAX[[SR3L[[1]][[1]]]]];
SR3RRT  = SR3location[[Range[15, 17]]]*ConvertToM;
SR3ang = SR3location[[Range[9, 11]]];
SR3G = Cylinder[{SR3XYZ, SR3XYZ - SR3RRT[[3]]*Reverse[{0, Cos[2Pi - SR3ang[[1]]*Pi/180], Sin[ 2Pi - SR3ang[[1]]*Pi/180]}]}, SR3RRT[[1]]]
 Graphics3D[{SR3G, HAM5G},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

BS

BSL = Position[H1ZEMAX, "BS"];
BSlocation = Take[H1ZEMAX[[BSL[[1]][[1]]]]];
BSXYZ =Reverse[BSlocation[[Range[6, 8]]]]*ConvertToM;
BSRRT  = BSlocation[[Range[13, 14]]]*ConvertToM;
BSang = BSlocation[[Range[9, 11]]];
BSG = Cylinder[{BSXYZ, BSXYZ - BSRRT[[2]]*Reverse[{0, Cos[-Pi/4], Sin[-Pi/4]}]}, BSRRT[[1]]]
 Graphics3D[{BSG, BSC2G, Blue, Sphere[BSXYZ, 25*ConvertToM]},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

ITMY

ITMYL = Position[H1ZEMAX, "ITMY HR REF"];
ITMYlocation = Take[H1ZEMAX[[ITMYL[[1]][[1]]]]];
ITMYXYZ =Reverse[ITMYlocation[[Range[6, 8]]]]*ConvertToM;

ITMYL = Position[H1ZEMAX, "ITMY"];
ITMYlocation = Take[H1ZEMAX[[ITMYL[[1]][[1]]]]];
ITMYRRT  = {170, 200}*ConvertToM;
ITMYang = ITMYlocation[[Range[9, 11]]];
ITMYG = Cylinder[{ITMYXYZ, ITMYXYZ - ITMYRRT[[2]]*{0, 1, 0}}, ITMYRRT[[1]]];
 Graphics3D[{ITMYG, BSC1G, Blue, Sphere[ITMYXYZ, 25*ConvertToM]},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

ITMX

ITMXL = Position[H1ZEMAX, "ITMX HR REF"];
ITMXlocation = Take[H1ZEMAX[[ITMXL[[1]][[1]]]]];
ITMXXYZ =Reverse[ITMXlocation[[Range[6, 8]]]]*ConvertToM;

ITMXL = Position[H1ZEMAX, "ITMX"];
ITMXlocation = Take[H1ZEMAX[[ITMXL[[1]][[1]]]]];
ITMXRRT  = {170, 200}*ConvertToM;
ITMXang = ITMXlocation[[Range[9, 11]]];
ITMXG = Cylinder[{ITMXXYZ, ITMXXYZ - ITMXRRT[[2]]*{1, 0, 0}}, ITMXRRT[[1]]];
 Graphics3D[{ITMXG, BSC3G, Blue, Sphere[ITMXXYZ, 25*ConvertToM]},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

Viewports

HAM4VP6L = Position[H1ZEMAX, "HAM4 VP6"];
HAM4VP6location = Take[H1ZEMAX[[HAM4VP6L[[1]][[1]]]]];
HAM4VP6XYZ =(Reverse[HAM4VP6location[[Range[6, 8]]]]*ConvertToM + HAM4XYZ);
HAM4VP6G = Cylinder[{{0,0,0}, {0,0,0} - {12.7, 0, 0}}*ConvertToM, 63.5*ConvertToM];
HAM4VP6G = Rotate[HAM4VP6G, 13.92 Degree, {0, 0, 1}];
HAM4VP6G = Translate[HAM4VP6G, HAM4VP6XYZ];
 Graphics3D[{HAM4VP6G, HAM4G},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

Hartmann BS Optics

HSBS1L = Position[H1ZEMAX, "HARTMANN-ITMY STEER M1"];
HSBS1location = Take[H1ZEMAX[[HSBS1L[[1]][[1]]]]];
HSBS1XYZ =Reverse[HSBS1location[[Range[6, 8]]]] *ConvertToM;
HSBS1G = Sphere[HSBS1XYZ, 20*ConvertToM];

HSBS2L = Position[H1ZEMAX, "ITMY PO STEER M1"];
HSBS2location = Take[H1ZEMAX[[HSBS2L[[1]][[1]]]]];
HSBS2XYZ =Reverse[HSBS2location[[Range[6, 8]]]] *ConvertToM;
HSBS2G = Sphere[HSBS2XYZ, 20*ConvertToM];

HSBS2L = Position[H1ZEMAX, "ITMY PO TEL M1"];
HSBS2location = Take[H1ZEMAX[[HSBS2L[[1]][[1]]]]];
HSBS3XYZ =Reverse[HSBS2location[[Range[6, 8]]]]*ConvertToM ;
HSBS3G = Sphere[HSBS3XYZ, 20*ConvertToM];

HSBS2L = Position[H1ZEMAX, "ITMY PO TEL M2"];
HSBS2location = Take[H1ZEMAX[[HSBS2L[[1]][[1]]]]];
HSBS4XYZ =Reverse[HSBS2location[[Range[6, 8]]]] *ConvertToM;
HSBS4G = Sphere[HSBS4XYZ, 20*ConvertToM];

HSBS2L = Position[H1ZEMAX, "HARTMANN-ITMY BS STEER"];
HSBS2location = Take[H1ZEMAX[[HSBS2L[[1]][[1]]]]];
HSBS5XYZ =Reverse[HSBS2location[[Range[6, 8]]]]*ConvertToM ;
HSBS5G = Sphere[HSBS5XYZ, 20*ConvertToM];

 Graphics3D[{HSBS5G, HAM4G},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

IFO

Graphics3D[{HAM4VP6G,ITMXG,ITMYG, BSG, SR3G, SR2G, BSC1G, BSC2G, BSC3G, HAM4G, HAM5G, Blue, HSBS1G, HSBS2G, HSBS3G, HSBS4G, HSBS5G},ViewPoint-> {0,0, Infinity}, ViewVertical-> {0, 0, 1}]

HAM4

Graphics3D[{HAM4VP6G, SR2G, HAM4G,Blue, HSBS1G, HSBS2G, HSBS3G, HSBS4G, HSBS5G}, ViewPoint -> {0, 0, Infinity}, ViewVertical -> {0, 0, 1}]

PastSR2 = Norm[HSBS1XYZ - SR2XYZ]+Norm[HSBS1XYZ - HAM4VP6XYZ]

Optical Layout - Y Arm
Requirements
1. HWS = Conjugate Plane of BS
2. Beam Size at ITM = Sqrt[2]*ModeSize
3. Magnification = 17.5x
4. Put the telescope outside the vacuum??
5. Beam is approximately collimated at the HWS

H1-ITMY Probe Beam

Intialize

Set probe beam size at ITM = Sqrt[2]*17.5*CCD_diameter

BeamSizeAtITM=127/1000

nFusedSilica840 = 1.445;
CPthick = 130*ConvertToM;
n1 = 1.0;
TotalDeMagnification = +17.5;

Create Function to Determine Output Radius and Beam Size

OutRW[wIN_, RIN_, λIN_,ABCD_]:= Module[{},
Q1Inv = (1/RIN)- ((I λIN)/(Pi wIN^2));
Q1a = 1/Q1Inv;
Q2temp = ((ABCD[[1]][[1]]Q1a + ABCD[[1]][[2]])/(ABCD[[2]][[1]]Q1a + ABCD[[2]][[2]]));
Q2Inv = 1/Q2temp;
ROut = Re[Q2Inv]^-1;
wOut = Sqrt[((-λIN)/(Pi Im[Q2Inv]))];
Return[{ROut, wOut}];
];
FindNearestWaist[wIN_, RIN_, λIN_]:= Module[{},
Q1Inv = (1/RIN)- ((I λIN)/(Pi wIN^2));
Q1a = 1/Q1Inv;


Q2temp = Q1a;
Q2tempR = Re[Q2temp]+x0;
Q2tempI = Im[Q2temp];
Q2InvR = (Q2tempR/(Q2tempR^2 + Q2tempI^2));


ROut = Q2InvR^(-1);

solnIN = Solve[1/ROut == 0, x0];


Return[x0/.solnIN[[1]]];
];

PlotTheBeamPropagation[ABCDList_, RWInputMode_, λIN_]:=Module[{},
z0qq = 0;
plotList = {};
count1 = 0;
BStempList = {};

For[ii=1, ii≤ Length[ABCDList], ii++, {
M1IN = ABCDList[[ii]];

plotAll = Plot[0*x, {x, 0, 0.00001}];
If[M1IN[[1]][[2]]≠ 0, {
MOut = (({{1, 0}, {0, 1}}));
For[jj=1, jj≤ ii-1, jj++, {
MOut = ABCDList[[jj]].MOut;
}];
MOut = (({{1, zIN-z0qq}, {0, 1}})).MOut;

eqOut = OutRW[RWInputMode[[2]], RWInputMode[[1]], λIN, MOut];

plotA = Plot[eqOut[[2]]*1000, {zIN, z0qq, ABCDList[[ii]][[1]][[2]]+z0qq}, PlotRange-> {{0, ABCDList[[ii]][[1]][[2]]+z0qq}, {0, Automatic}}, Frame->True, FrameLabel-> {"Propagation Distance (m)", "Beam Size (mm)"}, PlotStyle-> If[Mod[count1, 2]==0, Blue, Red]];

For[lll = z0qq, lll≤ ABCDList[[ii]][[1]][[2]]+z0qq, lll = lll + ABCDList[[ii]][[1]][[2]]/100.0, {

BStempList = Append[BStempList, {zIN, eqOut[[2]]*1000}/. zIN-> lll];


}];



plotList = Append[plotList, plotA];
count1 = count1 + 1;

z0qq = ABCDList[[ii]][[1]][[2]]+z0qq;

plotAll = Show[plotA,plotAll];

}];
}];
plotAll = plotList[[1]];
For[ii=2, ii≤ Length[plotList], ii++, {
plotAll = Show[plotList[[ii]],plotAll];
}];

Show[plotAll, PlotRange-> {{0,z0qq}, {0, Automatic}}]
];
z0qq=.;
PlotTheBeamPropagation[{(({{1, 1.0}, {0, 1}})), (({{1, 0}, {-1/1.5, 1}})), (({{1, 1.7}, {0, 1}})), (({{1, 0}, {-1/0.5, 1}})), (({{1, 0.7}, {0, 1}}))}, {1000, 0.001}, 840*10^-9]

Optical Path - determine solution: AirLens = 1.1337 m, VacLens = 0.79362 m

Given Physical Parameters

Distances

D0901179 - Thickness of SR2 = 75mm

SR2RRT
SR2RRT[[3]] = 0.075

SR3XYZ

SR2toSR3 = 15.443 (*T1300438*)
SR3toITMY = 24.7836 (*T1300438*)
SR3toBS = 19.4313 (*T1300438*)

Curvatures - E0900095, T080268, E080518,

ITMRoC = 1934/nFusedSilica840;
SR3RoC = 36;
SR2RoC = 6.43;

ABCD Matrices - Not relying on solution

Requirement 2: The mode from the ITM, including the size of Sqrt[2]*CavityModeSize is defined below

Test Gaussian beam in and out

MFromITM = (({{1, 0}, {0, (nFusedSilica840/n1)}})).(({{1, 0}, {((n1 - nFusedSilica840)/(SR2RoC *nFusedSilica840)), (n1/nFusedSilica840)}})).(({{1, SR2toSR3}, {0, 1}})). (({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR3toITMY}, {0, 1}}));
MToITM =  (({{1, SR3toITMY}, {0, 1}})).(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR2toSR3}, {0, 1}})).(({{1, 0}, {((nFusedSilica840-n1)/(-SR2RoC *n1)), (nFusedSilica840/n1)}})).(({{1, 0}, {0, (n1/nFusedSilica840)}}));

wIN = BeamSizeAtITM;
RIN = ITMRoC;
pIn = {RIN, wIN}
λIN = 840*10^-9;

pOut = OutRW[wIN, RIN, λIN,MFromITM];
pOut = OutRW[pOut[[2]], -pOut[[1]], λIN,MToITM]

ITM to SR2_AR: From HR surface outside ITM to AR surface outside SR2 [in m]

ITMToSR2AR = (({{1, 0}, {0, (nFusedSilica840/n1)}})).(({{1, SR2RRT[[3]]}, {0, 1}})).(({{1, 0}, {((n1 - nFusedSilica840)/(SR2RoC *nFusedSilica840)), (n1/nFusedSilica840)}})).(({{1, SR2toSR3}, {0, 1}})).(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR3toITMY}, {0, 1}}));
MatrixForm[ITMToSR2AR]

BS to SR2_AR: ABCD Matrix - From BS to AR

BSToSR2AR =(({{1, 0}, {0, (nFusedSilica840/n1)}})).(({{1, SR2RRT[[3]]}, {0, 1}})).(({{1, 0}, {((n1 - nFusedSilica840)/(SR2RoC *nFusedSilica840)), (n1/nFusedSilica840)}})).(({{1, SR2toSR3}, {0, 1}})).(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR3toBS}, {0, 1}}));
MatrixForm[BSToSR2AR]

Locate Nearest Conjugate Plane of BS

ConjPlaneM = (({{1, x}, {0, 1}})).(BSToSR2AR);
soln1 = Solve[ConjPlaneM[[1]][[2]]==0];
xConj = x/. soln1[[1]];
magnifcationBS1 = ConjPlaneM[[1]][[1]]/. soln1[[1]];
Print["Nearest conjugate plane of BS is ", xConj, "m from AR surface"]
Print["DeMagnification to First BS Conjugate Plane =  ", 1/magnifcationBS1, " x"]
Print["Residual magnification to still be achieved =  ", TotalDeMagnification*magnifcationBS1, " x"]

ResidualDeMag = TotalDeMagnification*magnifcationBS1
NearestConjPlane = xConj

Extract Optical Information - Conjugate plane, beam size At SR2

solnList = {2.2826, 2.9695, 3.6715, 1.13366, 0.79362}
H1YSoln = solnList;

Get Output Mode from SR2

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840*10^-9;
OutRW[wIN, RIN, λIN,{{1, 100}, {0, 1}}];
pOut = OutRW[wIN, RIN, λIN,ITMToSR2AR];
OutR = pOut[[1]];
OutW = pOut[[2]];
OutModeSR2AR = pOut;
FindNearestWaist[pOut[[2]], pOut[[1]], λIN];

Must create mode - matching optics to produce this mode
Must create imaging optics to image the limiting aperture - the beam splitter - correctly

Print["Probe beam size at AR surface of SR2 = ", OutW*1000, " mm."]
Print["Curvature at AR surface of SR2 = ", OutR*1000, " mm."]

SOutput = OutModeSR2AR[[1]]

Demonstrate curvature sign by propagating a short distance further

pOutdz = OutRW[wIN, RIN, λIN,{{1, 0.1}, {0, 1}}.ITMToSR2AR];
OutRdz = pOutdz[[1]];
OutWdz = pOutdz[[2]];
Print["Probe beam size at AR surface of SR2 = ", OutW*1000, " mm."]
Print["Probe beam size 100mm out of IFO from SR2-AR = ", OutWdz*1000, " mm."]

Analytic Solution for imaging a conjugate plane with an arbitrary telescope

M1 = Simplify[(({{1, z3}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z2}, {0, 1}})).(({{1, 0}, {-1/f1, 1}})).(({{1, z1}, {0, 1}})).(({{A11, B11}, {C11, D11}}))];

The final curvature is, roughly, the sum of the defocus of the telescope (fixed) and the curvature of the input beam (if, curvature ~= 0, then depends on position of telescope)
Therefore, easiest to vary z1 to set final C = 0

soln1 = Solve[M1[[2]][[1]] == 0, z1][[1]];
ASOLN1 = soln1;

The imaging is most directly solved by moving the final plane along the optical axis.
Therefore, easiest to vary z3 to set final B = 0

M1B= Simplify[M1[[1]][[2]]/. soln1];
soln2 =Simplify[ Solve[M1B == 0, z3][[1]]];
ASOLN2 = soln2;

This leaves z2 to solve for the magnification

M1A = Simplify[M1[[1]][[1]]/. soln1]/.soln2;
soln3 = Solve[M1A == m1a, z2][[1]];
ASOLN3 = soln3;
ZList =Simplify[ ((({z1, z2, z3}/. soln1)/. soln2)/. soln3)/. m1a-> 1/M];
ZListOriginal =Simplify[ ((({z1, z2, z3}/. soln1)/. soln2)/. soln3)];

C = 0 requirement

Simplify[soln1 /.soln3]

A = M requirement

Simplify[soln3]

B = 0 requirement

Simplify[soln2/. soln3]

ZList2 = Simplify[ZList /. {A11-> ITMToSR2AR[[1]][[1]], B11-> ITMToSR2AR[[1]][[2]],C11-> ITMToSR2AR[[2]][[1]],D11-> ITMToSR2AR[[2]][[2]]}]

Solve for the solution - ITMY: Magnification  = +ve!!

Just work out the solutions for the ABCD matrices as follows:
1. B = 0
2. C = 0
3. A = 1/17.5

eq1: Collimation requirement - at output in paraxial approximation
eq2: Total magnification requirement
eq3: Imaging requirement

Use the analytic solution determined above.
Notes: 
1. The combination of lenses was chosen by first constraining the parameter space of possible solutions to those where 
  a. {z1, z2, z3} >= {3, 5, 0}m 
  b. dS/dT of z2 <= 20x smaller than the maximum defocus error of 4E-4 (m^(-1))
2. Once this region was identified the lenses were chosen to have convenient values and to yield a solution close to the minimum Abs[dS/dT]

Hence: f1 = 1.0m, f2 = 3.0m

Z1A = solnList[[3]];
Z2A = solnList[[2]];
Z3A = solnList[[1]];
F1A = solnList[[5]];
F2A = solnList[[4]];
AllSolns2 = { F1A, F2A, Z1A, Z2A, Z3A}
TestSolutionSubstitution = {f1-> solnList[[5]], f2-> solnList[[4]], z1->solnList[[3]], z2-> solnList[[2]], z3-> solnList[[1]]}

ABCD Matrices - Dependent on solution

SR2_AR to HWS: ABCD matrix - get the imaging matrix part - after SR2

ImgMatrix = Simplify[(({{1, z3}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z2}, {0, 1}})).(({{1, 0}, {-1/f1, 1}})).(({{1, z1}, {0, 1}}))];
SR2ARToHWS = ImgMatrix/. TestSolutionSubstitution;
MatrixForm[SR2ARToHWS]

HWS to SR2_AR: ABCD matrix - get the imaging matrix part - after SR2

ImgMatrix = Simplify[(({{1, z1 }, {0, 1}})).(({{1, 0}, {-1/f1, 1}})).(({{1, z2}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z3}, {0, 1}}))];
HWSToSR2AR = ImgMatrix/. TestSolutionSubstitution;
MatrixForm[HWStoSR2AR]

ITM to HWS: ABCD matrix coming out of the vacuum - For Gaussian beam mode checking

ITMToHWS = SR2ARToHWS.ITMToSR2AR;
MatrixForm[ITMToHWS];
HWSToITM =(({{1, SR3toITMY}, {0, 1}})) .(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR2toSR3}, {0, 1}})).(({{1, 0}, {((nFusedSilica840-n1)/(-SR2RoC *n1)), (nFusedSilica840/n1)}})).(({{1, SR2RRT[[3]]}, {0, 1}})).(({{1, 0}, {0, (n1/nFusedSilica840)}})).(({{1, z1 }, {0, 1}})).(({{1, 0}, {-1/f1, 1}})). (({{1, z2}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z3}, {0, 1}}))/. TestSolutionSubstitution;
MITM = (({{1, 0}, {((-2)/(-ITMRoC)), 1}}));
MHWSv = ITMToHWS.MITM.HWSToITM

Alignment Laser?

wIN = 0.001
RIN = 1000
OutRW[wIN, RIN, λIN,MHWSv]
0.001*17.5

BS to HWS: ABCD matrix - Beam splitter to HWS - Shows Imaging

BS2HWS = SR2ARToHWS.BSToSR2AR;
MatrixForm[BS2HWS]
MatrixForm[ITMToHWS]

Requirement 1: The Beam Splitter to HWS ABCD matrix has a value of B approximately equal to 0.

Do a Gaussian beam propagation

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOut = OutRW[wIN, RIN, λIN,ITMToHWS]
OutR = pOut[[1]];
OutW = pOut[[2]]
GaussianMag = wIN/OutW

Requirement 3: The beam size at the HWS is 17.5x smaller than at the beam splitter

Define the Input Mode

InputMode = {OutW, -OutR}
InputModeFromHWS = InputMode;

Propagate modes - from ITM to outside and then back again

Propagate the ITM mode from the ITM to the HWS

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
Print["Probe beam size at ITM = ", wIN*1000, " mm."]
Print["Curvature at ITM = ", RIN, " m."]
pOut = OutRW[wIN, RIN, λIN,ITMToHWS];
Print["Probe beam size at HWS = ", pOut[[2]]*1000, " mm."]
Print["Curvature at HWS = ", pOut[[1]]*1000, " mm."]
Print["ITM_Beam_Size/HWS_Beam_Size = ", wIN/pOut[[2]], " x."]

Propagate the HWS mode to the ITM

pOut = OutRW[pOut[[2]], -pOut[[1]], λIN,HWSToITM];
Print["Probe beam size at ITM = ", pOut[[2]]*1000, " mm."]
Print["Curvature at ITM = ", pOut[[1]], " m."]

Propagate from HWS mode to HWS

Print["Probe beam size at ITM = ", InputMode[[1]]*1000, " mm."]
Print["Curvature at ITM = ", InputMode[[2]], " m."]
pOut = OutRW[InputMode[[1]], InputMode[[2]], λIN,MHWSv];
Print["Probe beam size at HWS = ", pOut[[2]]*1000, " mm."]
Print["Curvature at HWS = ", pOut[[1]], " m."]

Plot the Solution

Plot solution from SR2_AR to HWS

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;

p1 = Plot[(OutRW[wIN, RIN, λIN,(({{1, x}, {0, 1}})).ITMToSR2AR])[[2]]*1000, {x, 0, Z1A}, PlotStyle-> {Red,Thickness[0.01]}];
BS1 = Table[{x,(OutRW[wIN, RIN, λIN,(({{1, x}, {0, 1}})).ITMToSR2AR])[[2]]*1000}, {x, 0, Z1A, Z1A/200}];
p2 = Plot[(OutRW[wIN, RIN, λIN,(({{1, x-(Z1A)}, {0, 1}})).(({{1, 0}, {-1/F1A, 1}})).(({{1, (Z1A)}, {0, 1}})).ITMToSR2AR])[[2]]*1000, {x, (Z1A), (Z1A)+Z2A}, PlotStyle-> {Blue,Thickness[0.01]}];
BS2 =  Table[{x,(OutRW[wIN, RIN, λIN,(({{1, x-(Z1A)}, {0, 1}})).(({{1, 0}, {-1/F1A, 1}})).(({{1, (Z1A)}, {0, 1}})).ITMToSR2AR])[[2]]*1000}, {x, (Z1A), (Z1A)+Z2A, (Z2A/200)}];
p3 = Plot[(OutRW[wIN, RIN, λIN,(({{1, x-(Z1A)-Z2A}, {0, 1}})).(({{1, 0}, {-1/F2A, 1}})).(({{1, Z2A}, {0, 1}})).(({{1, 0}, {-1/F1A, 1}})).(({{1, (Z1A)}, {0, 1}})).ITMToSR2AR])[[2]]*1000, {x, (Z1A)+Z2A, (Z1A)+Z2A+Z3A}, PlotStyle-> {Red,Thickness[0.01]}];
BS3 =  Table[{x,(OutRW[wIN, RIN, λIN,(({{1, x-(Z1A)-Z2A}, {0, 1}})).(({{1, 0}, {-1/F2A, 1}})).(({{1, Z2A}, {0, 1}})).(({{1, 0}, {-1/F1A, 1}})).(({{1, (Z1A)}, {0, 1}})).ITMToSR2AR])[[2]]*1000}, {x, (Z1A)+Z2A, (Z1A)+Z2A+Z3A, Z3A/200}];
ITMYPlottedSolutionFromSR2AR = Show[p1, p2, p3, Frame-> True, FrameLabel-> {"Propagation Distance from SR2_AR (m)", "Beam size (mm)", "Beam propagating from SR2_AR to HWS"}, PlotRange-> {{0, 12}, Automatic}]

BS = BS1;
For[ii=1, ii≤ Length[BS2], ii++, {
BS = Append[BS, BS2[[ii]]];
}];
For[ii=1, ii≤ Length[BS3], ii++, {
BS = Append[BS, BS3[[ii]]];
}];
Export["/Advanced_LIGO/TCS/Hartmann sensor/H1_HWSY_beam_size.txt", BS, "Table"];

(OutRW[wIN, RIN, λIN,(({{1, x}, {0, 1}})).ITMToSR2AR])[[2]]*1000/. x-> Z1A
(OutRW[wIN, RIN, λIN,(({{1, x-(Z1A)}, {0, 1}})).(({{1, 0}, {-1/F1A, 1}})).(({{1, (Z1A)}, {0, 1}})).ITMToSR2AR])[[2]]*1000/. x-> Z1A

Z1A

Plot from ITMY to SR2

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
plotR = {{0, 55}, {0, 80}};
p1 = Plot[(OutRW[wIN, RIN, λIN,(({{1, x}, {0, 1}}))])[[2]]*1000, {x, 0, SR3toITMY}, PlotStyle-> {Red,Thickness[0.01]}, PlotRange-> plotR];
p2 = Plot[(OutRW[wIN, RIN, λIN,(({{1, x-SR3toITMY}, {0, 1}})).(({{1, 0}, {-2/SR3RoC, 1}})).(({{1, SR3toITMY}, {0, 1}}))])[[2]]*1000, {x, SR3toITMY, SR3toITMY + SR2toSR3}, PlotStyle-> {Blue,Thickness[0.01]}];
p3 = Plot[(OutRW[wIN, RIN, λIN,(({{1, x-SR3toITMY-SR2toSR3}, {0, 1}})).ITMToSR2AR])[[2]]*1000, {x,  SR3toITMY + SR2toSR3,SR3toITMY + SR2toSR3+ Z1A}, PlotStyle-> {Red,Thickness[0.01]}, PlotRange->plotR ];
ITMYPlottedSolutionFromITMY = Show[p1, p2, p3, Frame-> True, FrameLabel-> {"Propagation Distance from ITMY (m)", "Beam size (mm)", "Beam propagating from ITMY to F1 (substrate of SR2 not shown)"}]

Plot the whole propagation

z1=.;

TestSolutionSubstitution

xConj

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
plot1 = PlotTheBeamPropagation[{(({{1, SR3toITMY}, {0, 1}})) ,(({{1, 0}, {((-2)/SR3RoC), 1}})),(({{1, SR2toSR3}, {0, 1}})),(({{1, 0}, {((n1-nFusedSilica840)/(SR2RoC *nFusedSilica840)), (n1/nFusedSilica840)}})),(({{1, SR2RRT[[3]]}, {0, 1}})),(({{1, 0}, {0, (nFusedSilica840/n1)}})),(({{1, z1 }, {0, 1}})),(({{1, 0}, {-1/f1, 1}})), (({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z3}, {0, 1}}))}/. TestSolutionSubstitution, {RIN, wIN}, 840*10^-9];
Show[plot1, FrameLabel-> {"Propagation Distance (m)", "Beam Size (mm)", "Propagation from ITM to SR3, thru SR2, thru F1 and F2 and onto HWS"}]

MHWSv
RIN = -4.0;
wIN = 0.5*10^-3;
nm = 10^-9;
pOut = OutRW[wIN, RIN, 532nm,MHWSv];
Print["Probe beam size at HWS = ", pOut[[2]]*1000, " mm."]
Print["Curvature at HWS = ", pOut[[1]]*1000, " mm."]
Print["ITM_Beam_Size/HWS_Beam_Size = ", wIN/pOut[[2]], " x."]

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
plot1 = PlotTheBeamPropagation[{(({{1, z3}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f1, 1}})),(({{1, z1 }, {0, 1}})),(({{1, SR2RRT[[3]]}, {0, 1}})),(({{1, 0}, {0, (n1/nFusedSilica840)}})),(({{1, 0}, {((nFusedSilica840-n1)/(-SR2RoC *n1)), (nFusedSilica840/n1)}})),(({{1, SR2toSR3}, {0, 1}})),(({{1, 0}, {((-2)/SR3RoC), 1}})),(({{1, SR3toITMY}, {0, 1}})),MITM,(({{1, SR3toITMY}, {0, 1}})) ,(({{1, 0}, {((-2)/SR3RoC), 1}})),(({{1, SR2toSR3}, {0, 1}})),(({{1, 0}, {((n1-nFusedSilica840)/(SR2RoC *nFusedSilica840)), (n1/nFusedSilica840)}})),(({{1, SR2RRT[[3]]}, {0, 1}})),(({{1, 0}, {0, (nFusedSilica840/n1)}})),(({{1, z1 }, {0, 1}})),(({{1, 0}, {-1/f1, 1}})), (({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z3}, {0, 1}}))}/. TestSolutionSubstitution, {1000.0, 5*10^-4}, 0.0532nm];
A1 = Show[plot1, FrameLabel-> {"Propagation Distance (m)", "Beam Size (mm)", "Propagation from ITM to SR3, thru SR2, thru F1 and F2 and onto HWS"}, PlotRange-> {0, 5}]

Pi*(5*10^-4)^2./(532nm)

plot1 = PlotTheBeamPropagation[{(({{1, SR3toBS-2}, {0, 1}})) ,(({{1, 0}, {((-2)/SR3RoC), 1}})),(({{1, SR2toSR3}, {0, 1}})),(({{1, 0}, {((n1-nFusedSilica840)/(SR2RoC *nFusedSilica840)), (n1/nFusedSilica840)}})),(({{1, SR2RRT[[3]]}, {0, 1}})),(({{1, 0}, {0, (nFusedSilica840/n1)}})),(({{1, z1 }, {0, 1}})),(({{1, 0}, {-1/f1, 1}})), (({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z3}, {0, 1}}))}/. TestSolutionSubstitution, {10^4, 0.004}, 532*10^-9];
plot2 = PlotTheBeamPropagation[{(({{1, SR3toBS-2}, {0, 1}})) ,(({{1, 0}, {((-2)/SR3RoC), 1}})),(({{1, SR2toSR3}, {0, 1}})),(({{1, 0}, {((n1-nFusedSilica840)/(SR2RoC *nFusedSilica840)), (n1/nFusedSilica840)}})),(({{1, SR2RRT[[3]]}, {0, 1}})),(({{1, 0}, {0, (nFusedSilica840/n1)}})),(({{1, z1 }, {0, 1}})),(({{1, 0}, {-1/f1, 1}})), (({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z3}, {0, 1}}))}/. TestSolutionSubstitution, {10^4, 0.004}, 532*10^-12];
Show[plot1, plot2, PlotRange->All]

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
plot1 = PlotTheBeamPropagation[{(({{1, z3}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, 15}, {0, 1}}))}/. TestSolutionSubstitution, {1000.0, 0.5*10^-3}, 532nm];
Show[plot1, PlotRange-> All]

Tolerancing on Y arm solution

Get the Variable matrix

SR3toHWSVar =(({{1, Z3A+dz3}, {0, 1}})).(({{1, 0}, {-1/F2A -(dSf2), 1}})).(({{1, Z2A+dZ2A}, {0, 1}})).(({{1, 0}, {-1/F1A -(dSf1), 1}})).(({{1, Z1A+dz1}, {0, 1}})). (({{1, 0}, {0, (nFusedSilica840/n1)}})).(({{1, SR2RRT[[3]]}, {0, 1}})).(({{1, 0}, {((n1 - nFusedSilica840)/(SR2RoC *nFusedSilica840)), (n1/nFusedSilica840)}})).(({{1, SR2toSR3}, {0, 1}})). (({{1, 0}, {((-2)/SR3RoC), 1}}));
BStoHWSVar = SR3toHWSVar.(({{1, SR3toBS}, {0, 1}}));
ITMYtoHWSVar = SR3toHWSVar.(({{1, SR3toITMY}, {0, 1}}));
ITMXtoHWSVar = .;

Variation in imaging (B = 0) with length changes: (assume that FL variations can be compensated for) - NO ISSUES - Easy to work with!

M1 = BStoHWSVar.(({{1, x}, {0, 1}}));
B1 = M1[[1]][[2]];

Test variation in position of F1 - solve for the position of the conjugate plane - VERY SENSITIVE

Set all delta values = 0 except for dZ2A and dz1. Set these last two to the same value with opposite sign to model a lens out of position

B1A = Simplify[B1/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> X1, dz1-> -X1}];
Solve[B1A==0, x];
ConjPlaneDisplacement = x/. %[[1]];
NominalConjugatePlane = ConjPlaneDisplacement/. X1-> 0;
Print["Conjugate Plane is nominally out of position by ", NominalConjugatePlane, "m"];
RateCPD = (D[ConjPlaneDisplacement, X1])/. X1-> 0;
Print["Conjugate Plane displacement for +/- 5 mm displacement of F1 = ", RateCPD*1000*0.005, "mm"];
BTestF1OutOfPosition = {RateCPD*1000*0.005, "mm", 5, "mm"};

Hence the conjugate plane is VERY dependent on the position of F1: basically see a 10x displacement in the position of the conjugate plane for a displacement of F1 out of position.

Test variation in position of F2 - solve for the position of the conjugate plane - INSENSITIVE to displacement of F2 out of position

Set all delta values = 0 except for dZ2A and dz3. Set these last two to the same value with opposite sign to model a lens out of position

B1A = Simplify[B1/. {dz1-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> -X1, dz3-> X1}];
Solve[B1A==0, x];
ConjPlaneDisplacement = x/. %[[1]];
NominalConjugatePlane = ConjPlaneDisplacement/. X1-> 0;
Print["Conjugate Plane is nominally out of position by ", NominalConjugatePlane, "m"];
RateCPD = (D[ConjPlaneDisplacement, X1])/. X1-> 0;
Print["Conjugate Plane displacement for +/- 5 mm displacement of F2 = ", RateCPD*0.005, "m"];
BTestF2OutOfPosition = {RateCPD*1000*0.005, "mm", 5, "mm"};

Test variation in F1-F2 distance - SENSITIVE

B1A = Simplify[B1/. {dz1-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> -X1, dz3-> 0}];
Solve[B1A==0, x];
ConjPlaneDisplacement = x/. %[[1]];
NominalConjugatePlane = ConjPlaneDisplacement/. X1-> 0;
Print["Conjugate Plane is nominally out of position by ", NominalConjugatePlane, "m"];
RateCPD = (D[ConjPlaneDisplacement, X1])/. X1-> 0;
Print["Conjugate Plane displacement for +/- 5 mm variation in distance between F1 and F2 = ", RateCPD*0.005, "m"];
BTestF1F2Distance = {RateCPD*1000*0.005, "mm", 5, "mm"};

Hence the conjugate plane is sensitive to the distance between F2 and F1 - but that's due to F1 displacement

Test variation in position of HWS - SENSITIVE = GOOD for tuning. Possibly bad for fringes

Set all delta values = 0 except for dz3.

B1A = Simplify[B1/. {dz1-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> 0, dz3-> X1}];
Solve[B1A==0, x];
ConjPlaneDisplacement = x/. %[[1]];
NominalConjugatePlane = ConjPlaneDisplacement/. X1-> 0;
Print["Conjugate Plane is nominally out of position by ", NominalConjugatePlane, "m"];
RateCPD = (D[ConjPlaneDisplacement, X1])/. X1-> 0;
Print["Conjugate Plane displacement for +/- 5 mm displacement of F2 = ", RateCPD*0.005, "m"];
BTestF2HWSDistance = {RateCPD*1000*0.005, "mm", 5, "mm"};

Hence the conjugate plane is totally sensitive to the position of HWS - THIS IS A GOOD THING - means it is easy to tune.

Variation in magnification (A = 1/17.5) with length changes: (assume that FL variations can be compensated for) - NO ISSUES - Easy to work with!

M1 =ITMYtoHWSVar;
A1 = M1[[1]][[1]];

Test variation in position of F1 - determine the magnification - INSENSITIVE

Set all delta values = 0 except for dZ2A and dz1. Set these last two to the same value with opposite sign to model a lens out of position

A1A = Simplify[A1/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> X1, dz1-> -X1}];
NominalMagnification= 1/A1A/. X1-> 0;
Print["Nominal paraxial magnification is ", NominalMagnification, "m"];
RateVariationInMag = (D[A1A, X1])/. X1-> 0;
Print["Variation in magnification for +/- 5 mm displacement of F1 = ", RateVariationInMag*0.005, "x"];

Hence the paraxial magnification is VERY invariant to the position of F1

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMYtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> 0, dz1-> -0})];
BeamWAtHWSNom = pOutNom[[2]];
MOut = (ITMYtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> X1, dz1-> -X1});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamWAtHWSX1= pOutX1[[2]];
MagNom = BeamSizeAtITM/BeamWAtHWSNom;
MagX1 = BeamSizeAtITM/BeamWAtHWSX1;
Print["Absolute Variation in magnification for +5mm displacement = ", MagX1-MagNom, "x"];
Print["Relative Variation in magnification for +5mm displacement = ", ((MagX1-MagNom)/MagNom)*100, "%"];
ATestF1OutOfPosition = {Round[((MagX1-MagNom)/MagNom)*100*10]/10.0,"%",  5, "mm"};

Hence the Gaussian magnification is VERY invariant to the position of F1

Test variation in position of F2 - determine the magnification - MARGINALLY SENSITIVE

Set all delta values = 0 except for dL1 and dz3. Set these last two to the same value with opposite sign to model a lens out of position

A1A = Simplify[A1/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dZ2A-> -X1, dz1-> 0}];
NominalMagnification= 1/A1A/. X1-> 0;
Print["Nominal paraxial magnification is ", NominalMagnification, "m"];
RateVariationInMag = (D[A1A, X1])/. X1-> 0;
Print["Variation in magnification for +/- 5 mm displacement of F2 = ", RateVariationInMag*0.005, "x"];

Hence the paraxial magnification is VERY invariant to the position of F2

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMYtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> 0, dz1-> -0})];
BeamWAtHWSNom = pOutNom[[2]];
MOut = (ITMYtoHWSVar/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dZ2A-> -X1, dz1-> 0});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamWAtHWSX1= pOutX1[[2]];
MagNom = BeamSizeAtITM/BeamWAtHWSNom;
MagX1 = BeamSizeAtITM/BeamWAtHWSX1;
Print["Absolute Variation in magnification for +5mm displacement = ", MagX1-MagNom, "x"];
Print["Relative Variation in magnification for +5mm displacement = ", ((MagX1-MagNom)/MagNom)*100, "%"];
ATestF2OutOfPosition = {Round[((MagX1-MagNom)/MagNom)*100*10]/10.0,"%",  5, "mm"};

Hence the Gaussian magnification is insensitive to the position of F2

Test variation in F1-F2 distance: determine the magnification - INSENSITIVE

A1A = Simplify[A1/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> -X1, dz1-> 0}];
NominalMagnification= 1/A1A/. X1-> 0;
Print["Nominal paraxial magnification is ", NominalMagnification, "m"];
RateVariationInMag = (D[A1A, X1])/. X1-> 0;
Print["Variation in magnification for +/- 5 mm displacement of F2 = ", RateVariationInMag*0.005, "x"];

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMYtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> 0, dz1-> -0})];
BeamWAtHWSNom = pOutNom[[2]];
MOut = (ITMYtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> -X1, dz1-> 0});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamWAtHWSX1= pOutX1[[2]];
MagNom = BeamSizeAtITM/BeamWAtHWSNom;
MagX1 = BeamSizeAtITM/BeamWAtHWSX1;
Print["Absolute Variation in mag. for +5mm distance change between F1 and F2 = ", MagX1-MagNom, "x"];
Print["Relative Variation in mag. for +5mm distance change between F1 and F2 = ", ((MagX1-MagNom)/MagNom)*100, "%"];
ATestF1F2Distance = {Round[((MagX1-MagNom)/MagNom)*100*10]/10.0,"%",  5, "mm"};

Hence the Gaussian magnification is very insensitive to the distance between F1 and F2

Test variation in position of HWS - determine the magnification - INSENSITIVE to position of HWS: 1% variation per 15mm

Set all delta values = 0 except for dZ2A and dz1. Set these last two to the same value with opposite sign to model a lens out of position

A1A = Simplify[A1/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dZ2A-> 0, dz1-> -0}];
NominalMagnification= 1/A1A/. X1-> 0;
Print["Nominal paraxial magnification is ", NominalMagnification, "m"];
RateVariationInMag = (D[A1A, X1])/. X1-> 0;
Print["Variation in magnification for +/- 5 mm displacement of HWS = ", RateVariationInMag*0.005, "x"];

Hence the paraxial magnification is SUPER invariant to the position of F1

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMYtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> 0, dz1-> -0})];
BeamWAtHWSNom = pOutNom[[2]];
MOut = (ITMYtoHWSVar/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dZ2A-> 0, dz1-> -0});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamWAtHWSX1= pOutX1[[2]];
MagNom = BeamSizeAtITM/BeamWAtHWSNom;
MagX1 = BeamSizeAtITM/BeamWAtHWSX1;
Print["Absolute Variation in magnification for +5mm displacement = ", MagX1-MagNom, "x"];
Print["Relative Variation in magnification for +5mm displacement = ", ((MagX1-MagNom)/MagNom)*100, "%"];
ATestF2HWSDistance = {Round[((MagX1-MagNom)/MagNom)*100*10]/10.0,"%",  5, "mm"};

Hence the Gaussian magnification is MOSTLY invariant to the position of HWS

Variation in collimation (C = 0) with length changes: (assume that FL variations can be compensated for)  -

M1 =ITMYtoHWSVar;
C1 = M1[[2]][[1]];

Test variation in position of F1 - determine the defocus - INSENSITIVE

Set all delta values = 0 except for dZ2A and dz1. Set these last two to the same value with opposite sign to model a lens out of position

C1A = Simplify[C1/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> X1, dz1-> -X1}];
NominalDefocus= C1A/. X1-> 0;
Print["Nominal paraxial defocus is ", NominalDefocus, " m^-1"];
RateVariationInDefocus = (D[C1A, X1])/. X1-> 0;
Print["Variation in defocus for + 5 mm displacement of F1 = ", RateVariationInDefocus*0.005, " m^-1, or ", 1/(RateVariationInDefocus*0.005), "m"];

Hence the relative change in the paraxial defocus is quite large but the absolute change (with displacement of F1) is quite small, about 86km lens for 5mm displacement. At HWS this is not an issue.

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMYtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> 0, dz1-> -0})];
BeamRAtHWSNom = pOutNom[[1]];
MOut = (ITMYtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> X1, dz1-> -X1});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamRAtHWSX1= pOutX1[[1]];
Print["Variation in defocus of Gaussian Beam for +5mm displacement = ", 1/BeamRAtHWSX1-1/BeamRAtHWSNom, " m^-1", " or lens of = ", 1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), " m"];
CTestF1OutOfPosition = {1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), "m", 5 ,"mm"};

Hence the Gaussian curvature at HWS does not experience large changes in the radius of curvature that with motion of F1

Test variation in position of F2 - determine the defocus - INSENSITIVE

Set all delta values = 0 except for dZ2A and dz3. Set these last two to the same value with opposite sign to model a lens out of position

C1A = Simplify[C1/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dZ2A-> -X1, dz1-> 0}];
NominalDefocus= C1A/. X1-> 0;
Print["Nominal paraxial defocus is ", NominalDefocus, " m^-1"];
RateVariationInDefocus = (D[C1A, X1])/. X1-> 0;
Print["Variation in defocus for + 5 mm displacement of F2 = ", RateVariationInDefocus*0.005, " m^-1, or ", 1/(RateVariationInDefocus*0.005), "m"];

Hence the relative change in the paraxial defocus is quite large but the absolute change (with displacement of F2) is quite small, about 8km lens for 5mm displacement. At HWS this is not an issue.

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMYtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> 0, dz1-> -0})];
BeamRAtHWSNom = pOutNom[[1]];
MOut = (ITMYtoHWSVar/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dZ2A-> -X1, dz1-> 0});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamRAtHWSX1= pOutX1[[1]];
Print["Variation in defocus of Gaussian Beam for +5mm displacement = ", 1/BeamRAtHWSX1-1/BeamRAtHWSNom, " m^-1", "or lens of = ", 1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), " m"];
CTestF2OutOfPosition = {1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), "m", 5 ,"mm"};

Hence the Gaussian curvature at HWS does not experience large changes in the radius of curvature that with motion of F2

Test variation in F1-F2 distance: determine the defocus - INSENSITIVE

C1A = Simplify[C1/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> -X1, dz1-> 0}];
NominalDefocus= C1A/. X1-> 0;
Print["Nominal paraxial defocus is ", NominalDefocus, " m^-1"];
RateVariationInDefocus = (D[C1A, X1])/. X1-> 0;
Print["Variation in defocus for + 5 mm displacement of F2 = ", RateVariationInDefocus*0.005, " m^-1, or ", 1/(RateVariationInDefocus*0.005), "m"];

Hence the relative change in the paraxial defocus is quite large but the absolute change (with variation in F1-F2 distance) is quite small, about 8km lens for 5mm displacement. At HWS this is not an issue.

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMYtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> 0, dz1-> -0})];
BeamRAtHWSNom = pOutNom[[1]];
MOut = (ITMYtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> -X1, dz1-> 0});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamRAtHWSX1= pOutX1[[1]];
Print["Variation in defocus of Gaussian Beam for +5mm displacement = ", 1/BeamRAtHWSX1-1/BeamRAtHWSNom, " m^-1", "or lens of = ", 1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), " m"];
CTestF1F2Distance = {1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), "m", 5 ,"mm"};

Hence the Gaussian curvature at HWS does not experience large changes in the radius of curvature that with variation in F1-F2 distance

Test variation in position of HWS - determine the defocus - INSENSITIVE

Set all delta values = 0 except for dz3. Set these last two to the same value with opposite sign to model a lens out of position

C1A = Simplify[C1/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dZ2A-> 0, dz1-> -0}];
NominalDefocus= C1A/. X1-> 0;
Print["Nominal paraxial defocus is ", NominalDefocus, " m^-1"];
RateVariationInDefocus = (D[C1A, X1])/. X1-> 0;
Print["Variation in defocus for + 5 mm displacement of F1 = ", RateVariationInDefocus*0.005, " m^-1, or ", 1/(RateVariationInDefocus*0.005), "m"];

Hence the relative change in the paraxial defocus is ZERO

wIN = BeamSizeAtITM;
RIN = ITMRoC;
λIN = 840.*10^-9;
pOutNom = OutRW[wIN, RIN, λIN,(ITMYtoHWSVar/. {dz3-> 0, dSf2-> 0, dSf1-> 0, dZ2A-> 0, dz1-> -0})];
BeamRAtHWSNom = pOutNom[[1]];
MOut = (ITMYtoHWSVar/. {dz3-> X1, dSf2-> 0, dSf1-> 0, dZ2A-> 0, dz1-> -0});
pOutX1 = OutRW[wIN, RIN, λIN,(MOut/. {X1-> 0.005})];
BeamRAtHWSX1= pOutX1[[1]];
Print["Variation in defocus of Gaussian Beam for +5mm displacement = ", 1/BeamRAtHWSX1-1/BeamRAtHWSNom, " m^-1", "or lens of = ", 1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), " m"];
CTestF2HWSDistance = {1/(1/BeamRAtHWSX1-1/BeamRAtHWSNom), "m", 5 ,"mm"};

Hence the Gaussian curvature at HWS does not experience large changes in the radius of curvature that with variation in HWS position - expected since beam is approximately collimated here

Difficulty in mode-matching (MMOut = MMIn) -  Take Care with Mode Matching

InputModeFromHWSWR  = InputModeFromHWS;
InputModeFromHWSRW  = {InputModeFromHWS[[2]],InputModeFromHWS[[1]]}
OutputMode = OutRW[InputModeFromHWSRW[[2]], InputModeFromHWSRW[[1]], λIN,MHWSv]

10% variation in mode size, same curvature - no worries, well inside Rayleigh range ...

InputModeFromHWSRW  = {InputModeFromHWS[[2]],InputModeFromHWS[[1]]*1.1}
OutputMode = OutRW[InputModeFromHWSRW[[2]], InputModeFromHWSRW[[1]], λIN,MHWSv]

1.0m lens added to input beam, same mode size

InputModeFromHWSRW  = {((1/1.0) +InputModeFromHWS[[2]]^-1)^-1,InputModeFromHWS[[1]]}
OutputMode = OutRW[InputModeFromHWSRW[[2]], InputModeFromHWSRW[[1]], λIN,MHWSv]
Print["Difference in output curvature = ",((1/(OutputMode[[1]]))+(1/(InputModeFromHWSRW[[1]])))^-1] 
Print["Relative variation in output size with a 1m lens added to input = ",((OutputMode[[2]] - InputModeFromHWSRW[[2]])/(InputModeFromHWSRW[[2]]))*100, "%"]

Temperature change required to scatter 0.01% from HWS error

Simple model first of temperature sensitivity - assuming common Length and focal length changes
requirements: 
1. beam size at the reference mirror fits on a 3" mirror, maybe 4".
2. the input optics are shared between the probe and reference beams.
3. consider thermal changes in length and defocus
4. ABCD and Gaussian analysis

Look at all the telescopes with defocus and infer maximum acceptable temp. change

Input Mode is ...

OutModeSR2AR

Length sensitive matrices - main telescope between SR2 and SR3

tempLens =(({{1, 0}, {0, (nFusedSilica840/n1)}})).(({{1, SR2RRT[[3]]}, {0, 1}})).(({{1, 0}, {((n1 - nFusedSilica840)/(SR2RoC *nFusedSilica840)), (n1/nFusedSilica840)}}));
EFLSR2 = -1/tempLens[[2]][[1]];

ITMToSR2ARdX =Simplify[ (({{1, 0}, {0, (nFusedSilica840/n1)}})).(({{1, SR2RRT[[3]]}, {0, 1}})).(({{1, 0}, {((n1 - nFusedSilica840)/(SR2RoC *nFusedSilica840)), (n1/nFusedSilica840)}})).(({{1, SR2toSR3+dx}, {0, 1}})).(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR3toITMY}, {0, 1}}))];SR2ToITMARdX = Simplify[ (({{1, SR3toITMY}, {0, 1}})).(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR2toSR3+dx}, {0, 1}})).(({{1, 0}, {(( nFusedSilica840-n1)/((-1*SR2RoC) *n1)), (nFusedSilica840/n1)}})).(({{1, SR2RRT[[3]]}, {0, 1}})).(({{1, 0}, {0, (n1/nFusedSilica840)}}))];
SR2ARToHWSdx = Simplify[(({{1, z3+dx3}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z2+dx2}, {0, 1}})).(({{1, 0}, {-1/f1, 1}})).(({{1, z1 +dx1}, {0, 1}}))]/. TestSolutionSubstitution;
HWSToSR2ARdx = Simplify[(({{1, z1 +dx1}, {0, 1}})).(({{1, 0}, {-1/f1, 1}})).(({{1, z2+dx2}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z3+dx3}, {0, 1}}))]/. TestSolutionSubstitution;

MatrixForm[MITM];
SR2toSR2dx =ITMToSR2ARdX.MITM.SR2ToITMARdX;
HWStoHWSdx = SR2ARToHWSdx.ITMToSR2ARdX.MITM.SR2ToITMARdX.HWSToSR2ARdx/. {dx1-> 0, dx2-> 0, dx3-> 0};

Check the matrices

SR2toSR2 = SR2toSR2dx/. dx-> 0;
OutModeSR2AR
InModeSR2AR= {-OutModeSR2AR[[1]], OutModeSR2AR[[2]]};
pOut = OutRW[InModeSR2AR[[2]], InModeSR2AR[[1]], λIN,SR2toSR2]

HWStoHWS = HWStoHWSdx/. dx-> 0;
InputModeFromHWS
pOut = OutRW[InputModeFromHWS[[1]], InputModeFromHWS[[2]], λIN,HWStoHWS];
Reverse[pOut]

A1: Get the defocus at ITM for a given displacement between SR2 and SR3

To scatter out 0.01% the maximum acceptable defocus is MaxDefocusAtITMAll m^-1; - T1000722

MaxAcceptableDefocusAtITM = MaxDefocusAtITMAll;

dxTest = 0.00025;
dtM = HWStoHWSdx/. dx-> dxTest;
InputModeFromHWS;
Print["Input defocus = ", 1/InputModeFromHWS[[2]], " m^-1"]
pOut = OutRW[InputModeFromHWS[[1]], InputModeFromHWS[[2]], λIN,dtM];
Print["Output defocus = ", 1/pOut[[1]], " m^-1"]
dS = 1/(-InputModeFromHWS[[2]]) - 1/pOut[[1]]; 
M2 = BeamSizeAtITM/InputModeFromHWS[[1]];

dSAtITM = dS/(M2^2);
Print["Defocus at ITM = ", dSAtITM, "m^-1 for ", Round[dxTest*10^6], "um displacement"]
dSdxSR2toSR3 = dSAtITM/dxTest;
Print["Max Acceptable defocus at ITM = ",MaxAcceptableDefocusAtITM," m^-1"]

A2: Limit apparent temperature change of the concrete between SR2 & SR3

Determine how much the previous dx would need to be scaled by to yield the maximum apparent defocus allowed at the ITM

dxScale = MaxAcceptableDefocusAtITM/dSAtITM

CoefficientThermalExpansionOfConcrete = 1.2*10^-5;
dT = dx/(L*CoefficientThermalExpansionOfConcrete)/. {L-> SR2toSR3, dx->dxTest*dxScale};
dTSR2toSR3 =  {Round[dT*100]/100., "K"};
Print["Maximum acceptable change in temperature = ", Round[dT*100]*10, " mK"]

Now we need a secondary set of optics that introduce the same level of thermal defocus and give a small beam at the reference retro-reflector.

B1: Get the defocus at ITM for a given displacement between SR2 and F1

HWStoHWSdx1 = SR2ARToHWSdx.ITMToSR2ARdX.MITM.SR2ToITMARdX.HWSToSR2ARdx/. {dx-> 0, dx2-> 0, dx3-> 0};

dxTest = 0.00025;
dtM = HWStoHWSdx1/. dx1-> dxTest;
InputModeFromHWS;
Print["Input defocus = ", 1/InputModeFromHWS[[2]], " m^-1"]
pOut = OutRW[InputModeFromHWS[[1]], InputModeFromHWS[[2]], λIN,dtM];
Print["Output defocus = ", 1/pOut[[1]], " m^-1"]
dS = 1/(-InputModeFromHWS[[2]]) - 1/pOut[[1]]; 
M2 = BeamSizeAtITM/InputModeFromHWS[[1]];
dSAtITM = dS/(M2^2);
Print["Defocus at ITM = ", dSAtITM, "m^-1 for ", Round[dxTest*10^6], "um displacement"]
dSdx = dSAtITM/dxTest;
Print["Max Acceptable defocus at ITM = ",MaxAcceptableDefocusAtITM," m^-1"]

B2: Limit apparent temperature change of the stainless between SR2 & F1

Determine how much the previous dx would need to be scaled by to yield the maximum apparent defocus allowed at the ITM

dxScale = MaxAcceptableDefocusAtITM/dSAtITM

CoefficientThermalExpansionOfStainlessSteel = 1.73*10^-5;
dT = dx/(L*CoefficientThermalExpansionOfStainlessSteel)/. {L-> Z1A , dx->dxTest*dxScale};
dTSR2toF1 =  {Round[dT*100]/100., "K"};
Print["Maximum acceptable change in temperature = ", Round[dT*100]*10, " mK"]

Now we need a secondary set of optics that introduce the same level of thermal defocus and give a small beam at the reference retro-reflector.

C1: Get the defocus at ITM for a given displacement between F1 and F2

HWStoHWSdx2 = SR2ARToHWSdx.ITMToSR2ARdX.MITM.SR2ToITMARdX.HWSToSR2ARdx/. {dx-> 0, dx1-> 0, dx3-> 0};

dxTest = 0.00025;
dtM = HWStoHWSdx2/. dx2-> dxTest;
InputModeFromHWS;
Print["Input defocus = ", 1/InputModeFromHWS[[2]], " m^-1"]
pOut = OutRW[InputModeFromHWS[[1]], InputModeFromHWS[[2]], λIN,dtM];
Print["Output defocus = ", 1/pOut[[1]], " m^-1"]
dS = 1/(-InputModeFromHWS[[2]]) - 1/pOut[[1]]; 
M2 = BeamSizeAtITM/InputModeFromHWS[[1]];
dSAtITM = dS/(M2^2);
Print["Defocus at ITM = ", dSAtITM, "m^-1 for ", Round[dxTest*10^6], "um displacement"]
dSdx = dSAtITM/dxTest;
Print["Max Acceptable defocus at ITM = ",MaxAcceptableDefocusAtITM," m^-1"]

C2: Limit apparent temperature change of the stainless between F1 & F2

Determine how much the previous dx would need to be scaled by to yield the maximum apparent defocus allowed at the ITM

dxScale = MaxAcceptableDefocusAtITM/dSAtITM

CoefficientThermalExpansionOfStainlessSteel = 1.73*10^-5;
dT = dx/(L*CoefficientThermalExpansionOfStainlessSteel)/. {L-> Z2A, dx->dxTest*dxScale};
dTF1toF2 =  {Round[dT*100]/100., "K"};
Print["Maximum acceptable change in temperature = ", Round[dT*100]*10, " mK"]

Now we need a secondary set of optics that introduce the same level of thermal defocus and give a small beam at the reference retro-reflector.

D1: Get the defocus at ITM for a given displacement between F2 and HWS

HWStoHWSdx3 = SR2ARToHWSdx.ITMToSR2ARdX.MITM.SR2ToITMARdX.HWSToSR2ARdx/. {dx-> 0, dx1-> 0, dx2-> 0};

dxTest = 0.00025;
dtM = HWStoHWSdx3/. dx3-> dxTest;
InputModeFromHWS;
Print["Input defocus = ", 1/InputModeFromHWS[[2]], " m^-1"]
pOut = OutRW[InputModeFromHWS[[1]], InputModeFromHWS[[2]], λIN,dtM];
Print["Output defocus = ", 1/pOut[[1]], " m^-1"]
dS = 1/(-InputModeFromHWS[[2]]) - 1/pOut[[1]]; 
M2 = BeamSizeAtITM/InputModeFromHWS[[1]];
dSAtITM = dS/(M2^2);
Print["Defocus at ITM = ", dSAtITM, "m^-1 for ", Round[dxTest*10^6], "um displacement"]
dSdx = dSAtITM/dxTest;
Print["Max Acceptable defocus at ITM = ",MaxAcceptableDefocusAtITM," m^-1"]

D2: Limit apparent temperature change of the stainless between F2 & HWS

Determine how much the previous dx would need to be scaled by to yield the maximum apparent defocus allowed at the ITM

dxScale = MaxAcceptableDefocusAtITM/dSAtITM

CoefficientThermalExpansionOfStainlessSteel = 1.73*10^-5;
dT = dx/(L*CoefficientThermalExpansionOfStainlessSteel)/. {L-> Z3A, dx->dxTest*dxScale};
dTF2toHWS =  {Round[dT*100]/100., "K"};
Print["Maximum acceptable change in temperature = ", Round[dT*100]*10, " mK"]

Now we need a secondary set of optics that introduce the same level of thermal defocus and give a small beam at the reference retro-reflector.

Report on Solution

Print["=== REPORT ON H1:ITMY-HWS SOLUTION ==="]
Print[" "];
Print["ITMY apparent Radius of Curvature = ", ITMRoC, "m"];
Print["ITM to SR3 optical distance = ", SR3toITMY, "m"];
Print["BS to SR3 optical distance = ", SR3toBS, "m"];
Print["SR3 Radius of Curvature = ", SR3RoC, "m"];
Print["SR3 to SR2 distance = ", SR2toSR3, "m"];
Print["SR2 Thickness = ", SR2RRT[[3]], "m"];
Print["SR2 Effective Focal Length = ",EFLSR2, "m"];
Print["SR2 Radius of Curvature = ", SR2RoC, "m"];
Print[" "];
Print["SR2_AR to F1 distance = ", Round[(Z1A )*1000]/1000.0, "m"];
Print["F1 focal length = ", Round[(F1A)*1000]/1000.0, "m"];
Print["F1 to F2 distance = ", Round[(Z2A )*1000]/1000.0, "m"];
Print["F2 focal length = ", Round[(F2A)*1000]/1000.0, "m"];
Print["F2 to HWS distance = ", Round[(Z3A )*1000]/1000.0, "m"];
Print[" "];
Print["Input beam size from HWS = ", InputModeFromHWS[[1]]*1000, "mm"];
Print["Input curvature from HWS = ", InputModeFromHWS[[2]], "m"];
Print[" "];
Print["----"]
Print["TEMPERATURE SENSITIVITY"]
Print["Maximum Acceptable defocus at the ITM = ", MaxAcceptableDefocusAtITM];
Print["SR2 to SR3 - maximum acceptable temperature change = ", dTSR2toSR3[[1]], dTSR2toSR3[[2]]]
Print["F1 to SR2 - maximum acceptable temperature change = ", dTSR2toF1[[1]], dTSR2toF1[[2]]]
Print["F1 to F2 - maximum acceptable temperature change = ", dTF1toF2[[1]], dTF1toF2[[2]]]
Print["F2 to HWS - maximum acceptable temperature change = ", dTF2toHWS[[1]], dTF2toHWS[[2]]]
Print[" "];
Print["----"]
Print["TOLERANCING"]
Print["Imaging (B=0): F1 out of Position by ", BTestF1OutOfPosition[[3]], BTestF1OutOfPosition[[4]], " moves conjugate plane by ", BTestF1OutOfPosition[[1]], BTestF1OutOfPosition[[2]]]
Print["Imaging (B=0): F2 out of Position by ", BTestF2OutOfPosition[[3]], BTestF2OutOfPosition[[4]], " moves conjugate plane by ", BTestF2OutOfPosition[[1]], BTestF2OutOfPosition[[2]]]
Print["Imaging (B=0): F1-F2 distance varies by ", BTestF1F2Distance[[3]], BTestF1F2Distance[[4]], " moves conjugate plane by ", BTestF1F2Distance[[1]], BTestF1F2Distance[[2]]]
Print["Imaging (B=0): HWS out of Position by ", BTestF2HWSDistance[[3]], BTestF2HWSDistance[[4]], " moves conjugate plane by ", BTestF2HWSDistance[[1]], BTestF2HWSDistance[[2]]]
Print[" "]
Print["Magnification (A=1/17.5): F1 out of Position by ", ATestF1OutOfPosition[[3]], ATestF1OutOfPosition[[4]], " varies magnification by ", ATestF1OutOfPosition[[1]], ATestF1OutOfPosition[[2]]]
Print["Magnification (A=1/17.5): F2 out of Position by ", ATestF2OutOfPosition[[3]], ATestF2OutOfPosition[[4]], " varies magnification by ", ATestF2OutOfPosition[[1]], ATestF2OutOfPosition[[2]]]
Print["Magnification (A=1/17.5): F1-F2 distance varies by ", ATestF1F2Distance[[3]], ATestF1F2Distance[[4]], " varies magnification by ", ATestF1F2Distance[[1]], ATestF1F2Distance[[2]]]
Print["Magnification (A=1/17.5): F2-HWS distance varies by ", ATestF2HWSDistance[[3]], ATestF2HWSDistance[[4]], " varies magnification by ", ATestF2HWSDistance[[1]], ATestF2HWSDistance[[2]]]
Print[" "]
Print["Collimation (C=0): F1 out of Position by ", CTestF1OutOfPosition[[3]], CTestF1OutOfPosition[[4]], " looks like an additional lens at HWS with f = ", CTestF1OutOfPosition[[1]], CTestF1OutOfPosition[[2]]]
Print["Collimation (C=0): F2 out of Position by ", CTestF2OutOfPosition[[3]], CTestF2OutOfPosition[[4]], " looks like an additional lens at HWS with f = ", CTestF2OutOfPosition[[1]], CTestF2OutOfPosition[[2]]]
Print["Collimation (C=0): F1-F2 distance varied by ", CTestF1F2Distance[[3]], CTestF1F2Distance[[4]], " looks like an additional lens at HWS with f = ", CTestF1F2Distance[[1]], CTestF1F2Distance[[2]]]
Print["Collimation (C=0): F2-HWS distance varied by ", CTestF2HWSDistance[[3]], CTestF2HWSDistance[[4]], " looks like an additional lens at HWS with f = ", CTestF2HWSDistance[[1]], CTestF2HWSDistance[[2]]]

Show[ITMYPlottedSolutionFromITMY]
Show[ITMYPlottedSolutionFromSR2AR]

Show[ITMYPlottedSolutionFromSR2AR, PlotRange-> {{0, 14}, {0, 20}}, GridLines-> Automatic]

Show[ITMYPlottedSolutionFromSR2AR, PlotRange-> {0, 0.4}]

Beam size at the mirrors

Create a Module to determine the beam size

GetTheBeamSize[ABCDList_, RWInputMode_, λIN_, zIn_]:=Module[{},
z0qq = 0;
plotList = {};
count1 = 0;
zOut = 0.0;
For[ii=1, ii≤ Length[ABCDList], ii++, {
M1IN = ABCDList[[ii]];

plotAll = Plot[0*x, {x, 0, 0.00001}];
If[M1IN[[1]][[2]]≠ 0, {
MOut = (({{1, 0}, {0, 1}}));
For[jj=1, jj≤ ii-1, jj++, {
MOut = ABCDList[[jj]].MOut;
}];
MOut = (({{1, zIn-z0qq}, {0, 1}})).MOut;

eqOut = OutRW[RWInputMode[[2]], RWInputMode[[1]], λIN, MOut];
If[zIn≥  z0qq, zOut = eqOut[[2]]];



count1 = count1 + 1;

z0qq = ABCDList[[ii]][[1]][[2]]+z0qq;



}];
}];
Return[zOut];
];
GetTheGouyPhase[ABCDList_, RWInputMode_, λIN_, zIn_]:=Module[{x0, waistRWqq},
z0qq = 0;
plotList = {};
count1 = 0;
zOut = 0.0;
For[ii=1, ii≤ Length[ABCDList], ii++, {
M1IN = ABCDList[[ii]];

plotAll = Plot[0*x, {x, 0, 0.00001}];
If[M1IN[[1]][[2]]≠ 0, {
MOut = (({{1, 0}, {0, 1}}));
For[jj=1, jj≤ ii-1, jj++, {
MOut = ABCDList[[jj]].MOut;
}];
MOut = (({{1, zIn-z0qq}, {0, 1}})).MOut;

eqOut = OutRW[RWInputMode[[2]], RWInputMode[[1]], λIN, MOut];
If[zIn≥  z0qq, zOut = eqOut];



count1 = count1 + 1;

z0qq = ABCDList[[ii]][[1]][[2]]+z0qq;



}];
}];

x0 = FindNearestWaist[zOut[[2]], zOut[[1]], λIN];
waistRWqq = OutRW[zOut[[2]], zOut[[1]], λIN,(({{1, x0}, {0, 1}}))];
GouyPhaseqq = ArcTan[((-x0)/((((Pi waistRWqq[[2]]^2)/(λIN)))))];
Return[GouyPhaseqq];
];

Import the data on the mirror positions

dataIN = Import["/Users/aidanbrooks/Documents/LIGO/documents/Advanced_LIGO/TCS/Hartmann sensor/modeling/H1HWSY_coords.txt", "TSV"];
totaldistance = 0;
distList = {};
nameList = {};
For[ii = 1, ii≤ Length[dataIN], ii++, {
d1 = dataIN[[ii]][[7]]/1000.0;
totaldistance = totaldistance + d1;
dataIN[[ii]] = Append[dataIN[[ii]], totaldistance];
distList = Append[distList , totaldistance];
nameList = Append[nameList, dataIN[[ii]][[1]]];
}];

dataIN

distList

distances to the mirrors from handoff at SR2_AR

distances = {0.0,893.8,2096.1,2516.4,3970.1,5793.6,6655.1,7055.1,9725.6,10411.7,10836.7,10886.7,11161.7,11997.4,12697.8,13322.8,13372.8,13472.8,14054.9,14254.9,14330.0,14730.0}/1000;
distList = distances;

Determine the beam size at each interface

beamSizeList = {};
For[kk=1, kk≤Length[distList], kk++, {
beamsize =GetTheBeamSize[{(({{1, z1 }, {0, 1}})),(({{1, 0}, {-1/f1, 1}})), (({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z3}, {0, 1}}))}/. TestSolutionSubstitution, OutModeSR2AR, 840*10^-9, distList[[kk]]];

beamSizeList = Append[beamSizeList, {distList[[kk]],beamsize}];
}];

beamSizeList

ListPlot[beamSizeList]

Export["/Users/aidanbrooks/Documents/LIGO/documents/Advanced_LIGO/TCS/Hartmann sensor/design/optical layout and modeling/H1HWSY_coords_beamsize.txt", beamSizeList,"Table"];

For[ii=1, ii≤Length[beamSizeList], ii++, Print[1000*beamSizeList[[ii]][[2]]]]

Determine the beam size at each interface

GouyPhase =GetTheGouyPhase[{(({{1, z1 }, {0, 1}})),(({{1, 0}, {-1/f1, 1}})), (({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z3}, {0, 1}}))}/. TestSolutionSubstitution, OutModeSR2AR, 840*10^-9, distList[[1]]]

beamSizeList = {};
GouyPhaseList = {};
GouyPhaseListDeg = {};
For[kk=1, kk≤Length[distList], kk++, {
beamsize =GetTheBeamSize[{(({{1, z1 }, {0, 1}})),(({{1, 0}, {-1/f1, 1}})), (({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z3}, {0, 1}}))}/. TestSolutionSubstitution, OutModeSR2AR, 840*10^-9, distList[[kk]]];

beamSizeList = Append[beamSizeList, {distList[[kk]],beamsize}];

GouyPhase =GetTheGouyPhase[{(({{1, z1 }, {0, 1}})),(({{1, 0}, {-1/f1, 1}})), (({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z3}, {0, 1}}))}/. TestSolutionSubstitution, OutModeSR2AR, 840*10^-9, distList[[kk]]];

GouyPhaseList = Append[GouyPhaseList, {distList[[kk]],GouyPhase}];
GouyPhaseListDeg = Append[GouyPhaseListDeg, {distList[[kk]],GouyPhase*180/Pi}];

}];

plot1 = ListPlot[GouyPhaseListDeg, Frame-> True, GridLines-> {Automatic, Table[ii, {ii, -90, 90, 15}]}];
plot2 = ListPlot[GouyPhaseListDeg, Frame-> True, Joined-> True];
Show[plot1, plot2]

GouyPhaseOptics = {};
For[jj = 1, jj≤ Length[GouyPhaseListDeg], jj++, {
GouyPhaseOptics = Append[GouyPhaseOptics, {nameList[[jj]], GouyPhaseListDeg[[jj]][[2]]}];
}];
MatrixForm[GouyPhaseOptics]

For[kk=1, kk≤ Length[distList], kk++, {
dataIN[[kk]] = Append[dataIN[[kk]], beamSizeList[[kk]][[2]]*1000.0];
}];

ListPlot[beamSizeList]

Export["/Users/aidanbrooks/Documents/LIGO/documents/Advanced_LIGO/TCS/Hartmann sensor/modeling/H1HWSY_coords_beamsize.txt", dataIN,"Table"];

MatrixForm[dataIN]

Mis - alignment

Determine the displacement of the beam axis at the HWS if a mirror is misaligned

MatrixList = {(({{1, z3}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f1, 1}})),(({{1, z1}, {0, 1}})), (({{1, 0}, {0, (nFusedSilica840/n1)}})),(({{1, SR2RRT[[3]]}, {0, 1}})),(({{1, 0}, {((n1 - nFusedSilica840)/(SR2RoC *nFusedSilica840)), (n1/nFusedSilica840)}})),(({{1, SR2toSR3}, {0, 1}})),(({{1, 0}, {((-2)/SR3RoC), 1}})),(({{1, SR3toITMY}, {0, 1}})),(({{1, 0}, {((-2)/(-ITMRoC)), 1}})),(({{1, SR3toITMY}, {0, 1}})) ,(({{1, 0}, {((-2)/SR3RoC), 1}})),(({{1, SR2toSR3}, {0, 1}})),(({{1, 0}, {((nFusedSilica840-n1)/(-SR2RoC *n1)), (nFusedSilica840/n1)}})),(({{1, SR2RRT[[3]]}, {0, 1}})),(({{1, 0}, {0, (n1/nFusedSilica840)}})),(({{1, z1 }, {0, 1}})),(({{1, 0}, {-1/f1, 1}})),(({{1, z2}, {0, 1}})),(({{1, 0}, {-1/f2, 1}})),(({{1, z3}, {0, 1}}))}/. z3-> z3+dz;
MatrixList = MatrixList/. TestSolutionSubstitution;
MList2 = Reverse[MatrixList];

Length[MatrixList]

AddAlignmentErrors[Matrices_, alignmentList_, AddPosns_, inputBeam_]:= Module[{}, 
count1 = 1;
outputbeam = inputBeam;
For[ii=1, ii≤ Length[Matrices], ii++, {
outputbeam = Matrices[[ii]].outputbeam;
If[count1≤ Length[AddPosns], {
If[ii== AddPosns[[count1]], {
outputbeam = outputbeam + alignmentList[[count1]];
count1 = count1+1;
}];
}];

}];
Return[outputbeam]]

p1 = Table[AddAlignmentErrors[MatrixList/. dz-> 0.063, {{0, -0.01*Sin[2*Pi*5*t]}}, {14}, {0.01, 0.0}] , {t, 0, 1,0.001}];
ListPlot[p1]

AddAlignmentErrors[MatrixList/. dz-> 0.063, {{0, -0.01}}, {14}, {0.01, 0.0}]

Alignment errors on SR2

SR2_AR to HWS: ABCD matrix - get the imaging matrix part - after SR2

ImgMatrix = Simplify[(({{1, z3}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z2}, {0, 1}})).(({{1, 0}, {-1/f1, 1}})).(({{1, z1}, {0, 1}}))];
SR2ARToHWS = ImgMatrix/. TestSolutionSubstitution;
MatrixForm[SR2ARToHWS]

HWS to SR2_AR: ABCD matrix - get the imaging matrix part - after SR2

ImgMatrix = Simplify[(({{1, z1 }, {0, 1}})).(({{1, 0}, {-1/f1, 1}})).(({{1, z2}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z3}, {0, 1}}))];
HWSToSR2AR = ImgMatrix/. TestSolutionSubstitution;
MatrixForm[HWStoSR2AR]

ITM to HWS: ABCD matrix coming out of the vacuum - For Gaussian beam mode checking

ITMToHWS = SR2ARToHWS.ITMToSR2AR;
MatrixForm[ITMToHWS];
MITM = (({{1, 0}, {((-2)/(-ITMRoC)), 1}}));
MHWSv = ITMToHWS.MITM.HWSToITM;

Line up matrices

HWStoSR2AR =(({{1, z1 }, {0, 1}})).(({{1, 0}, {-1/f1, 1}})). (({{1, z2}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z3}, {0, 1}}))/. TestSolutionSubstitution;
SR2toITM = (({{1, SR3toITMY}, {0, 1}})) .(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR2toSR3}, {0, 1}})).(({{1, 0}, {((nFusedSilica840-n1)/(-SR2RoC *n1)), (nFusedSilica840/n1)}})).(({{1, SR2RRT[[3]]}, {0, 1}})).(({{1, 0}, {0, (n1/nFusedSilica840)}}));
MITM;
ITMToSR2AR;
SR2ARToHWS;

Misalign output beam [m, radians] per radian

MisAlignedISI=Simplify[SR2ARToHWS.({dy, da} +(ITMToSR2AR.MITM.SR2toITM).({dy, da} + HWStoSR2AR.{y0, a0}))]

misalignment = D[MisAlignedISI, da]

misalignment[[1]]*1000*0.001

Alignment errors near F1

F1 to HWS: ABCD matrix - get the imaging matrix part

ImgMatrix = Simplify[(({{1, z3}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z2}, {0, 1}})).(({{1, 0}, {-1/f1, 1}}))];
F1ToHWS = ImgMatrix/. TestSolutionSubstitution;
MatrixForm[F1ToHWS]

HWS to F1: ABCD matrix - get the imaging matrix part - after SR2

ImgMatrix = Simplify[(({{1, 0}, {-1/f1, 1}})).(({{1, z2}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z3}, {0, 1}}))];
HWSToF1 = ImgMatrix/. TestSolutionSubstitution;
MatrixForm[HWStoF1]

ITM to HWS: ABCD matrix coming out of the vacuum - For Gaussian beam mode checking

MITM = (({{1, 0}, {((-2)/(-ITMRoC)), 1}}));

Line up matrices

F1toITM = (({{1, SR3toITMY}, {0, 1}})) .(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR2toSR3}, {0, 1}})).(({{1, 0}, {((nFusedSilica840-n1)/(-SR2RoC *n1)), (nFusedSilica840/n1)}})).(({{1, SR2RRT[[3]]}, {0, 1}})).(({{1, 0}, {0, (n1/nFusedSilica840)}})).(({{1, z1}, {0, 1}}))/. TestSolutionSubstitution;
ITMtoF1 =(({{1, z1}, {0, 1}})).ITMToSR2AR /. TestSolutionSubstitution;
MITM;
ITMToSR2AR;
SR2ARToHWS;

Misalign output beam [m, radians] per radian

MisAlignedISI=Simplify[F1ToHWS.({dy, da} +(ITMtoF1.MITM.F1toITM).({dy, da} + HWSToF1.{y0, a0}))]

misalignment = D[MisAlignedISI, dy]

misalignment[[1]]*1000*0.001

Misalignment of delay line mirrors = delta theta,
maximum range of delay line mirrors = 50 mm, 
displacement at HWS = delta_theta*50 mm*0.1 = 5mm*delta_theta
Maximum allowed displacement at HWS = spot pitch = 450 um,
Therefore delta_theta = 450um/5mm = 0.1 radians

MatrixList

Alignment

Create an ABCD plotting function

PlotABCDMatrixRayTrace[inputBeam_, ABCDList_]:= Module[{},
totalDist = 0;
currentABCDMatrix = {{1,0}, {0,1}};
PlotCount = 0;

For[ii=1, ii≤ Length[ABCDList], ii++, {
If[ABCDList[[ii]][[1]][[2]]==0,
{currentABCDMatrix = ABCDList[[ii]].currentABCDMatrix};
,{
oldABCDMatrix = currentABCDMatrix;
currentABCDMatrix = {{1, zQQ-totalDist}, {0, 1}}.currentABCDMatrix;

dzmax = ABCDList[[ii]][[1]][[2]];

zmin = totalDist;
zmax = totalDist+dzmax;
beamOut = currentABCDMatrix.inputBeam;

plotA = Plot[beamOut[[1]], {zQQ, zmin, zmax}, PlotRange-> {{0, zmax}, {-1,1}}, Frame-> True, PlotStyle-> Thick];


If[PlotCount == 0, {
plotAll = Show[plotA];
},{
plotAll = Show[plotAll, plotA,PlotRange-> {{0, zmax}, Automatic}];
}];


PlotCount = PlotCount+1;
currentABCDMatrix = ABCDList[[ii]].oldABCDMatrix;
totalDist = totalDist + dzmax;
}];

}];
Return[plotAll];
];

Basic Misalignment on Reflection from ITM

Plot the ABCD matrix

dtheta = 100*10^-6;
plot1 = PlotABCDMatrixRayTrace[{0.105, dtheta},MatrixListFromITM];
plot2 = PlotABCDMatrixRayTrace[{-0.0, dtheta},MatrixListFromITM];
plot3 = PlotABCDMatrixRayTrace[{-0.105, dtheta},MatrixListFromITM];
axisPlot = Plot[0,{x, 0, zmax}, PlotStyle-> Dashed] ;
oneinchMt = Plot[{12.5/1000+ 0*x, -12.5/1000+0*x},{x, 0, zmax}, PlotStyle-> Dashed] ;
twoinchMt = Plot[2*{12.5/1000+ 0*x, -12.5/1000+0*x},{x, 0, zmax},PlotStyle-> Dashed] ;
Show[plot1, plot2, plot3,axisPlot,oneinchMt,twoinchMt, PlotRange-> All]

Basic Misalignment on Transmission to ITM

Plot the ABCD matrix

dtheta = 0.010;
plot1 = PlotABCDMatrixRayTrace[{0.006, dtheta},MatrixListToITM];
plot2 = PlotABCDMatrixRayTrace[{-0.00, dtheta},MatrixListToITM];
plot3 = PlotABCDMatrixRayTrace[{-0.006, dtheta},MatrixListToITM];
axisPlot = Plot[0,{x, 0, zmax}, PlotStyle-> Dashed] ;
oneinchMt = Plot[{12.5/1000+ 0*x, -12.5/1000+0*x},{x, 0, zmax}, PlotStyle-> Dashed] ;
twoinchMt = Plot[2*{12.5/1000+ 0*x, -12.5/1000+0*x},{x, 0, zmax},PlotStyle-> Dashed] ;
Show[plot1, plot2, plot3,axisPlot,oneinchMt,twoinchMt, PlotRange-> All]

Full Transmission

Plot the ABCD matrix

MatrixListFull = MatrixListToITM;
For[ii=1, ii≤ Length[MatrixListFromITM], ii++, {MatrixListFull = Append[MatrixListFull, MatrixListFromITM[[ii]] ] ;}];

dtheta = 00.003;
dx = 0.000;
plot1 = PlotABCDMatrixRayTrace[{0.006+dx, dtheta},MatrixListFull];
plot2 = PlotABCDMatrixRayTrace[{-0.0+dx, dtheta},MatrixListFull];
plot3 = PlotABCDMatrixRayTrace[{-0.006+dx, dtheta},MatrixListFull];
axisPlot = Plot[0,{x, 0, zmax}, PlotStyle-> Dashed] ;
oneinchMt = Plot[{12.5/1000+ 0*x, -12.5/1000+0*x},{x, 0, zmax}, PlotStyle-> Dashed] ;
twoinchMt = Plot[2*{12.5/1000+ 0*x, -12.5/1000+0*x},{x, 0, zmax},PlotStyle-> Dashed] ;
Show[plot1, plot2, plot3,axisPlot,oneinchMt,twoinchMt, PlotRange-> All]
Show[plot1, plot2, plot3,axisPlot,oneinchMt,twoinchMt, PlotRange-> {{0, 4}, Automatic}]

Gaussian Width Fitting

GaussianE[xdata_]:= Module[{},
Return[130000Exp[-2*(((xdata-512))/(1024/Sqrt[2]))^2]];
]

data1 = Table[{ii, GaussianE[ii]}, {ii, 1, 1024}];
dx = 430/12.0;
data2 =  Table[{ii, GaussianE[ii] +20*Sqrt[GaussianE[ii]]* (Random[]-0.5)}, {ii, 10, 1024, dx}];

p1 = ListPlot[data1, Joined-> True, Frame-> True, PlotRange-> {0, 150000}];
p2 = ListPlot[data2, Joined-> False, Frame-> True, PlotRange-> {0, 150000}];
p3 = Show[p1, p2];
p4 = p3;

For[kk=1, kk≤ 8, kk++,{
amp1 = GaussianE[(kk-0.345)/8*1024]/130000;
data1 = Table[{ii, amp1*GaussianE[ii]}, {ii, 1, 1024}];
dx = 430/12.0;


data2 =  Table[{ii, amp1*GaussianE[ii] +20*Sqrt[GaussianE[ii]]* (Random[]-0.5)}, {ii, 10, 1024, dx}];

p1 = ListPlot[data1, Joined-> True, Frame-> True, PlotRange-> {0, 150000}, PlotStyle-> Hue[(kk-0.345)/8]];
p2 = ListPlot[data2, Joined-> False, Frame-> True, PlotRange-> {0, 150000},PlotStyle-> Hue[(kk-0.345)/8]];
p4 = Show[p4, p1, p2];
};
];
Show[p3]
Show[p4]

plot5 = DensityPlot[Exp[-2 (((x-512)^2 + (y-512)^2)/((1024/Sqrt[2])^2))], {x, 0, 1024}, {y,0, 1024}, PlotRange-> {0, 1}, Frame-> True,ColorFunction->"Rainbow", PlotPoints-> 100];
plot6 = Plot[0, {x, 0, 1024}, PlotRange-> {{0, 1024}, {0, 1024}}, GridLines-> {{}, Table[((kk-0.345)/8)*1024, {kk,8}]}, Frame-> True];
Show[plot5]

Determine Range and Resolution of the Alignment depending on Actuator Location

Position 1 - Near the Viewport

f1 = .;
f2 = .;

MatrixListFromITMAlign = {(({{1, z3}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, dL1}, {0, 1}})),(({{1, dM1M2}, {0, 1}})),(({{1, z2-dL1-dM1M2}, {0, 1}})). (({{1, 0}, {-1/f1, 1}})).(({{1, z1 }, {0, 1}})).(({{1, 0}, {0, (nFusedSilica840/n1)}})).(({{1, SR2RRT[[3]]}, {0, 1}})).(({{1, 0}, {((n1 - nFusedSilica840)/(SR2RoC *nFusedSilica840)), (n1/nFusedSilica840)}})).(({{1, SR2toSR3}, {0, 1}})).(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR3toITMY}, {0, 1}})).(({{1, 0}, {((-2)/(-ITMRoC)), 1}}))}/. TestSolutionSubstitution;
nElemM = Length[MatrixListFromITMAlign]

Actuator 2 is closer to Viewport. Actuator 1 is close to F1

InputBeam = {yIN, alphaIN};

BeamAtActuator2 = MatrixListFromITMAlign[[nElemM]].InputBeam;
deltaActuator2 = {0, dAlpha2};
BeamAfterActuator2 = Simplify[BeamAtActuator2 + deltaActuator2];

BeamAtActuator1= MatrixListFromITMAlign[[nElemM-1]].BeamAfterActuator2;
deltaActuator1 = {0, dAlpha1};
BeamAfterActuator1 =  Simplify[BeamAtActuator1 + deltaActuator1];

BeamAtHWS = Simplify[MatrixListFromITMAlign[[nElemM-2]].BeamAfterActuator1]

BeamAtHWS1 = Chop[Simplify[BeamAtHWS /. alphaIN-> 0]]

Get the Transfer Matrix Between the Actuators and the HWS

DdAlpha2 = Simplify[D[BeamAtHWS1, dAlpha2]]
DdAlpha1 = Simplify[D[BeamAtHWS1, dAlpha1]]

- If the same angle is added to both actuators but the magnitude is reversed then the displacement at the HWS is just equal to the angle times the distance between the actuators (no magnification?)

Sensing Matrix

MOut = {{DdAlpha2[[1]], DdAlpha1[[1]]}, {DdAlpha2[[2]], DdAlpha1[[2]]}}

Nominal Values for the Actuators: CASE 1 (Position 1)

Case 1 : (HWSY STEER M8 and HWSY STEER M10)
dM1dM2 = 0.325
dL1 = 0.83564

Case1 = {dM1M2-> 0.325, dL1-> 0.83564}

MOutCase = MOut/.Case1;
InvMout = Inverse[MOutCase]

Work Out Resolution

Resolution of angles & displacements on HWS
deltaAlphaMin = 0.7 urad
resolution at HWS = Sqrt[C2A^2 + C1A^2]*0.7 urad

C2A = DdAlpha2/. Case1;
C1A =DdAlpha1/. Case1;
resol = 0.7*10^-6;
ResolutionAtHWS = Sqrt[C1A^2 + C2A^2]*(resol)

mrad = 0.001;
RangeAct = 70 mrad;
RangeAtHWS = (Abs[C2A] + Abs[C1A])*RangeAct

Work Out Ranges

Range of angles & displacements on HWS
deltaAlphaMax = 70 umrad
Range at HWS = (Abs[C2A] + Abs[C1A])*70 mrad

mrad = 0.001;
RangeAct = 70 mrad;
RangeAtHWS = (Abs[C2A] + Abs[C1A])*RangeAct

Angular Range at HWS keeping displacement fixed
dHWS = 0;
Solve dx == 0 and substitute in solution for angles

eqSolve1 = deltaxHWS == MOutCase[[1]].{dAlpha2, dAlpha1};
soln1 = Solve[eqSolve1, dAlpha1][[1]];
angleOut =({MOutCase[[2]].{dAlpha2, dAlpha1}, dAlpha1}/. soln1)/. deltaxHWS-> 0;
angleOut1=angleOut /. dAlpha2-> RangeAct;
Print["Maximum range on Alpha 2 requires " <> ToString[100*Abs[angleOut1[[2]]]/RangeAct]<> "% on Alpha 1"];

soln2 = Solve[eqSolve1, dAlpha2][[1]];
angleOut =({MOutCase[[2]].{dAlpha2, dAlpha1}, dAlpha2}/. soln2)/. deltaxHWS-> 0;
angleOut2=angleOut /. dAlpha1-> RangeAct;
Print["Maximum range on Alpha 1 requires " <> ToString[100*Abs[angleOut2[[2]]]/RangeAct]<> "% on Alpha 2"];
Print[" "];

If[Abs[angleOut2[[2]]]>Abs[angleOut1[[2]]], {
Print["Maximum range of Angle 1 saturates Angle 2"];Print["Angular range is limited by Actuator 2"];
Print["Maximum range (keeping displacement on HWS constant) = "<> ToString[angleOut1[[1]]]<> " radians"];
MaxRangeAngle = Abs[angleOut1[[1]]];
} , {
Print["Maximum range of Angle 2 saturates Angle 1"];
Print["Angular range is limited by Actuator 1"];
Print["Maximum range (keeping displacement on HWS constant) = "<> ToString[angleOut2[[1]]]<> " radians"];
MaxRangeAngle = Abs[angleOut2[[1]]];
}
];

Displacement Range at HWS keeping angle fixed
Solve dtheta == 0 and substitute in solution for angles

eqSolve1 = deltaThetaHWS == MOutCase[[2]].{dAlpha2, dAlpha1};
soln1 = Solve[eqSolve1, dAlpha1][[1]];
dxOut =({MOutCase[[1]].{dAlpha2, dAlpha1}, dAlpha1}/. soln1)/. deltaThetaHWS-> 0;
dxOut1=dxOut /. dAlpha2-> RangeAct;
Print["Maximum range on Alpha 2 requires " <> ToString[100*Abs[dxOut1[[2]]]/RangeAct]<> "% on Alpha 1"];

soln2 = Solve[eqSolve1, dAlpha2][[1]];
dxOut =({MOutCase[[1]].{dAlpha2, dAlpha1}, dAlpha2}/. soln2)/. deltaThetaHWS-> 0;
dxOut2=dxOut /. dAlpha1-> RangeAct;
Print["Maximum range on Alpha 1 requires " <> ToString[100*Abs[dxOut2[[2]]]/RangeAct]<> "% on Alpha 2"];
Print[" "];

If[Abs[dxOut2[[2]]]>Abs[dxOut1[[1]]], {
Print["Maximum range of Angle 1 saturates Angle 2"];Print["Angular range is limited by Actuator 2"];
Print["Maximum range (keeping angle on HWS constant) = "<> ToString[dxOut1[[1]]*1000]<> " mm"];
MaxRangeDisplacement = Abs[dxOut1[[1]]];
} , {
Print["Maximum range of Angle 2 saturates Angle 1"];
Print["Angular range is limited by Actuator 1"];
Print["Maximum range (keeping angle on HWS constant) = "<> ToString[dxOut2[[1]]*1000]<> " mm"];
MaxRangeDisplacement = Abs[dxOut2[[1]]];
}
];

MaxRangeCase1 = {MaxRangeDisplacement, MaxRangeAngle}
Case1Result = {ResolutionAtHWS, MaxRangeCase1}

Nominal Values for the Actuators: CASE 2 (Position 1)

Case 2 : (HWSY STEER M9 and HWSY STEER M10)
dM1dM2 = 0.05
dL1 = 0.83564 + 0.275

Case2 = {dM1M2-> 0.325-0.275, dL1-> 0.83564+0.275}

MOutCase = MOut/.Case2
InvMout = Inverse[MOutCase]

Work Out Resolution

Resolution of angles & displacements on HWS
deltaAlphaMin = 0.7 urad
resolution at HWS = Sqrt[C2A^2 + C1A^2]*0.7 urad

C2A = DdAlpha2/. Case2;
C1A =DdAlpha1/. Case2;
resol = 0.7*10^-6;
ResolutionAtHWS = Sqrt[C1A^2 + C2A^2]*(resol)

mrad = 0.001;
RangeAct = 70 mrad;
RangeAtHWS = (Abs[C2A] + Abs[C1A])*RangeAct

Work Out Ranges

Range of angles & displacements on HWS
deltaAlphaMax = 70 umrad
Range at HWS = (Abs[C2A] + Abs[C1A])*70 mrad

mrad = 0.001;
RangeAct = 70 mrad;
RangeAtHWS = (Abs[C2A] + Abs[C1A])*RangeAct

Angular Range at HWS keeping displacement fixed
dHWS = 0;
Solve dx == 0 and substitute in solution for angles

eqSolve1 = deltaxHWS == MOutCase[[1]].{dAlpha2, dAlpha1};
soln1 = Solve[eqSolve1, dAlpha1][[1]];
angleOut =({MOutCase[[2]].{dAlpha2, dAlpha1}, dAlpha1}/. soln1)/. deltaxHWS-> 0;
angleOut1=angleOut /. dAlpha2-> RangeAct;
Print["Maximum range on Alpha 2 requires " <> ToString[100*Abs[angleOut1[[2]]]/RangeAct]<> "% on Alpha 1"];

soln2 = Solve[eqSolve1, dAlpha2][[1]];
angleOut =({MOutCase[[2]].{dAlpha2, dAlpha1}, dAlpha2}/. soln2)/. deltaxHWS-> 0;
angleOut2=angleOut /. dAlpha1-> RangeAct;
Print["Maximum range on Alpha 1 requires " <> ToString[100*Abs[angleOut2[[2]]]/RangeAct]<> "% on Alpha 2"];
Print[" "];

If[Abs[angleOut2[[2]]]>Abs[angleOut1[[2]]], {
Print["Maximum range of Angle 1 saturates Angle 2"];Print["Angular range is limited by Actuator 2"];
Print["Maximum range (keeping displacement on HWS constant) = "<> ToString[angleOut1[[1]]]<> " radians"];
MaxRangeAngle = Abs[angleOut1[[1]]];
} , {
Print["Maximum range of Angle 2 saturates Angle 1"];
Print["Angular range is limited by Actuator 1"];
Print["Maximum range (keeping displacement on HWS constant) = "<> ToString[angleOut2[[1]]]<> " radians"];
MaxRangeAngle = Abs[angleOut2[[1]]];
}
];

Displacement Range at HWS keeping angle fixed
Solve dtheta == 0 and substitute in solution for angles

eqSolve1 = deltaThetaHWS == MOutCase[[2]].{dAlpha2, dAlpha1};
soln1 = Solve[eqSolve1, dAlpha1][[1]];
dxOut =({MOutCase[[1]].{dAlpha2, dAlpha1}, dAlpha1}/. soln1)/. deltaThetaHWS-> 0;
dxOut1=dxOut /. dAlpha2-> RangeAct;
Print["Maximum range on Alpha 2 requires " <> ToString[100*Abs[dxOut1[[2]]]/RangeAct]<> "% on Alpha 1"];

soln2 = Solve[eqSolve1, dAlpha2][[1]];
dxOut =({MOutCase[[1]].{dAlpha2, dAlpha1}, dAlpha2}/. soln2)/. deltaThetaHWS-> 0;
dxOut2=dxOut /. dAlpha1-> RangeAct;
Print["Maximum range on Alpha 1 requires " <> ToString[100*Abs[dxOut2[[2]]]/RangeAct]<> "% on Alpha 2"];
Print[" "];

If[Abs[dxOut2[[2]]]>Abs[dxOut1[[1]]], {
Print["Maximum range of Angle 1 saturates Angle 2"];Print["Angular range is limited by Actuator 2"];
Print["Maximum range (keeping angle on HWS constant) = "<> ToString[dxOut1[[1]]*1000]<> " mm"];
MaxRangeDisplacement = Abs[dxOut1[[1]]];
} , {
Print["Maximum range of Angle 2 saturates Angle 1"];
Print["Angular range is limited by Actuator 1"];
Print["Maximum range (keeping angle on HWS constant) = "<> ToString[dxOut2[[1]]*1000]<> " mm"];
MaxRangeDisplacement = Abs[dxOut2[[1]]];
}
];

MaxRangeCase2 = {MaxRangeDisplacement, MaxRangeAngle}
Case2Result = {ResolutionAtHWS, MaxRangeCase2}

Position 2 - Near the Viewport

f1 = .;
f2 = .;

MatrixListFromITMAlign = {(({{1, z3-dL1-dM1M2}, {0, 1}})),(({{1, dM1M2}, {0, 1}})),(({{1, dL1}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, z2}, {0, 1}})). (({{1, 0}, {-1/f1, 1}})).(({{1, z1 }, {0, 1}})).(({{1, 0}, {0, (nFusedSilica840/n1)}})).(({{1, SR2RRT[[3]]}, {0, 1}})).(({{1, 0}, {((n1 - nFusedSilica840)/(SR2RoC *nFusedSilica840)), (n1/nFusedSilica840)}})).(({{1, SR2toSR3}, {0, 1}})).(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR3toITMY}, {0, 1}})).(({{1, 0}, {((-2)/(-ITMRoC)), 1}}))}/. TestSolutionSubstitution;
nElemM = Length[MatrixListFromITMAlign]

Actuator 2 is closer to Viewport. Actuator 1 is closer to F1

InputBeam = {yIN, alphaIN};

BeamAtActuator2 = MatrixListFromITMAlign[[nElemM]].InputBeam;
deltaActuator2 = {0, dAlpha2};
BeamAfterActuator2 = Simplify[BeamAtActuator2 + deltaActuator2];

BeamAtActuator1= MatrixListFromITMAlign[[nElemM-1]].BeamAfterActuator2;
deltaActuator1 = {0, dAlpha1};
BeamAfterActuator1 =  Simplify[BeamAtActuator1 + deltaActuator1];

BeamAtHWS = Simplify[MatrixListFromITMAlign[[nElemM-2]].BeamAfterActuator1]

BeamAtHWS1 = Chop[Simplify[BeamAtHWS /. alphaIN-> 0]]

Get the Transfer Matrix Between the Actuators and the HWS

DdAlpha2 = Simplify[D[BeamAtHWS1, dAlpha2]]
DdAlpha1 = Simplify[D[BeamAtHWS1, dAlpha1]]

- If the same angle is added to both actuators but the magnitude is reversed then the displacement at the HWS is just equal to the angle times the distance between the actuators (no magnification?)

Sensing Matrix

MOut = {{DdAlpha2[[1]], DdAlpha1[[1]]}, {DdAlpha2[[2]], DdAlpha1[[2]]}}

Nominal Values for the Actuators: CASE 3 (Position 2)

Case 3 : (HWSY STEER M11 and HWSY STEER M12)
dM1dM2 = 0.625
dL1 = 0.7

Case3 = {dM1M2-> 0.625, dL1-> 0.7}

MOutCase = MOut/.Case3;
InvMout = Inverse[MOutCase]

Work Out Resolution

Resolution of angles & displacements on HWS
deltaAlphaMin = 0.7 urad
resolution at HWS = Sqrt[C2A^2 + C1A^2]*0.7 urad

C2A = DdAlpha2/. Case3;
C1A =DdAlpha1/. Case3;
resol = 0.7*10^-6;
ResolutionAtHWS = Sqrt[C1A^2 + C2A^2]*(resol)

mrad = 0.001;
RangeAct = 70 mrad;
RangeAtHWS = (Abs[C2A] + Abs[C1A])*RangeAct

Work Out Ranges

Range of angles & displacements on HWS
deltaAlphaMax = 70 umrad
Range at HWS = (Abs[C2A] + Abs[C1A])*70 mrad

mrad = 0.001;
RangeAct = 70 mrad;
RangeAtHWS = (Abs[C2A] + Abs[C1A])*RangeAct

Angular Range at HWS keeping displacement fixed
dHWS = 0;
Solve dx == 0 and substitute in solution for angles

eqSolve1 = deltaxHWS == MOutCase[[1]].{dAlpha2, dAlpha1};
soln1 = Solve[eqSolve1, dAlpha1][[1]];
angleOut =({MOutCase[[2]].{dAlpha2, dAlpha1}, dAlpha1}/. soln1)/. deltaxHWS-> 0;
angleOut1=angleOut /. dAlpha2-> RangeAct;
Print["Maximum range on Alpha 2 requires " <> ToString[100*Abs[angleOut1[[2]]]/RangeAct]<> "% on Alpha 1"];

soln2 = Solve[eqSolve1, dAlpha2][[1]];
angleOut =({MOutCase[[2]].{dAlpha2, dAlpha1}, dAlpha2}/. soln2)/. deltaxHWS-> 0;
angleOut2=angleOut /. dAlpha1-> RangeAct;
Print["Maximum range on Alpha 1 requires " <> ToString[100*Abs[angleOut2[[2]]]/RangeAct]<> "% on Alpha 2"];
Print[" "];

If[Abs[angleOut2[[2]]]>Abs[angleOut1[[2]]], {
Print["Maximum range of Angle 1 saturates Angle 2"];Print["Angular range is limited by Actuator 2"];
Print["Maximum range (keeping displacement on HWS constant) = "<> ToString[angleOut1[[1]]]<> " radians"];
MaxRangeAngle = Abs[angleOut1[[1]]];
} , {
Print["Maximum range of Angle 2 saturates Angle 1"];
Print["Angular range is limited by Actuator 1"];
Print["Maximum range (keeping displacement on HWS constant) = "<> ToString[angleOut2[[1]]]<> " radians"];
MaxRangeAngle = Abs[angleOut2[[1]]];
}
];

Displacement Range at HWS keeping angle fixed
Solve dtheta == 0 and substitute in solution for angles

eqSolve1 = deltaThetaHWS == MOutCase[[2]].{dAlpha2, dAlpha1};
soln1 = Solve[eqSolve1, dAlpha1][[1]];
dxOut =({MOutCase[[1]].{dAlpha2, dAlpha1}, dAlpha1}/. soln1)/. deltaThetaHWS-> 0;
dxOut1=dxOut /. dAlpha2-> RangeAct;
Print["Maximum range on Alpha 2 requires " <> ToString[100*Abs[dxOut1[[2]]]/RangeAct]<> "% on Alpha 1"];

soln2 = Solve[eqSolve1, dAlpha2][[1]];
dxOut =({MOutCase[[1]].{dAlpha2, dAlpha1}, dAlpha2}/. soln2)/. deltaThetaHWS-> 0;
dxOut2=dxOut /. dAlpha1-> RangeAct;
Print["Maximum range on Alpha 1 requires " <> ToString[100*Abs[dxOut2[[2]]]/RangeAct]<> "% on Alpha 2"];
Print[" "];

If[Abs[dxOut2[[2]]]>Abs[dxOut1[[1]]], {
Print["Maximum range of Angle 1 saturates Angle 2"];Print["Angular range is limited by Actuator 2"];
Print["Maximum range (keeping angle on HWS constant) = "<> ToString[dxOut1[[1]]*1000]<> " mm"];
MaxRangeDisplacement = Abs[dxOut1[[1]]];
} , {
Print["Maximum range of Angle 2 saturates Angle 1"];
Print["Angular range is limited by Actuator 1"];
Print["Maximum range (keeping angle on HWS constant) = "<> ToString[dxOut2[[1]]*1000]<> " mm"];
MaxRangeDisplacement = Abs[dxOut2[[1]]];
}
];

MaxRangeCase3 = {MaxRangeDisplacement, MaxRangeAngle}
Case3Result = {ResolutionAtHWS, MaxRangeCase3}

Nominal Values for the Actuators: CASE 4 (Position 2)

Case 3 : (HWSY STEER M12 and HWSY STEER M13)
dM1dM2 = 0.625
dL1 = 0.7+0625

Case4 = {dM1M2-> 0.05, dL1-> 0.7+0.625}

MOutCase = MOut/.Case4;
InvMout = Inverse[MOutCase]

Work Out Resolution

Resolution of angles & displacements on HWS
deltaAlphaMin = 0.7 urad
resolution at HWS = Sqrt[C2A^2 + C1A^2]*0.7 urad

C2A = DdAlpha2/. Case4;
C1A =DdAlpha1/. Case4;
resol = 0.7*10^-6;
ResolutionAtHWS = Sqrt[C1A^2 + C2A^2]*(resol);

ResolutionCase4 = ResolutionAtHWS;

mrad = 0.001;
RangeAct = 70 mrad;
RangeAtHWS = (Abs[C2A] + Abs[C1A])*RangeAct

Work Out Ranges

Range of angles & displacements on HWS
deltaAlphaMax = 70 umrad
Range at HWS = (Abs[C2A] + Abs[C1A])*70 mrad

mrad = 0.001;
RangeAct = 70 mrad;
RangeAtHWS = (Abs[C2A] + Abs[C1A])*RangeAct

Angular Range at HWS keeping displacement fixed
dHWS = 0;
Solve dx == 0 and substitute in solution for angles

eqSolve1 = deltaxHWS == MOutCase[[1]].{dAlpha2, dAlpha1};
soln1 = Solve[eqSolve1, dAlpha1][[1]];
angleOut =({MOutCase[[2]].{dAlpha2, dAlpha1}, dAlpha1}/. soln1)/. deltaxHWS-> 0;
angleOut1=angleOut /. dAlpha2-> RangeAct;
Print["Maximum range on Alpha 2 requires " <> ToString[100*Abs[angleOut1[[2]]]/RangeAct]<> "% on Alpha 1"];

soln2 = Solve[eqSolve1, dAlpha2][[1]];
angleOut =({MOutCase[[2]].{dAlpha2, dAlpha1}, dAlpha2}/. soln2)/. deltaxHWS-> 0;
angleOut2=angleOut /. dAlpha1-> RangeAct;
Print["Maximum range on Alpha 1 requires " <> ToString[100*Abs[angleOut2[[2]]]/RangeAct]<> "% on Alpha 2"];
Print[" "];

If[Abs[angleOut2[[2]]]>Abs[angleOut1[[2]]], {
Print["Maximum range of Angle 1 saturates Angle 2"];Print["Angular range is limited by Actuator 2"];
Print["Maximum range (keeping displacement on HWS constant) = "<> ToString[angleOut1[[1]]]<> " radians"];
MaxRangeAngle = Abs[angleOut1[[1]]];
} , {
Print["Maximum range of Angle 2 saturates Angle 1"];
Print["Angular range is limited by Actuator 1"];
Print["Maximum range (keeping displacement on HWS constant) = "<> ToString[angleOut2[[1]]]<> " radians"];
MaxRangeAngle = Abs[angleOut2[[1]]];
}
];

Displacement Range at HWS keeping angle fixed
Solve dtheta == 0 and substitute in solution for angles

eqSolve1 = deltaThetaHWS == MOutCase[[2]].{dAlpha2, dAlpha1};
soln1 = Solve[eqSolve1, dAlpha1][[1]];
dxOut =({MOutCase[[1]].{dAlpha2, dAlpha1}, dAlpha1}/. soln1)/. deltaThetaHWS-> 0;
dxOut1=dxOut /. dAlpha2-> RangeAct;
Print["Maximum range on Alpha 2 requires " <> ToString[100*Abs[dxOut1[[2]]]/RangeAct]<> "% on Alpha 1"];

soln2 = Solve[eqSolve1, dAlpha2][[1]];
dxOut =({MOutCase[[1]].{dAlpha2, dAlpha1}, dAlpha2}/. soln2)/. deltaThetaHWS-> 0;
dxOut2=dxOut /. dAlpha1-> RangeAct;
Print["Maximum range on Alpha 1 requires " <> ToString[100*Abs[dxOut2[[2]]]/RangeAct]<> "% on Alpha 2"];
Print[" "];

If[Abs[dxOut2[[2]]]>Abs[dxOut1[[1]]], {
Print["Maximum range of Angle 1 saturates Angle 2"];Print["Angular range is limited by Actuator 2"];
Print["Maximum range (keeping angle on HWS constant) = "<> ToString[dxOut1[[1]]*1000]<> " mm"];
MaxRangeDisplacement = Abs[dxOut1[[1]]];
} , {
Print["Maximum range of Angle 2 saturates Angle 1"];
Print["Angular range is limited by Actuator 1"];
Print["Maximum range (keeping angle on HWS constant) = "<> ToString[dxOut2[[1]]*1000]<> " mm"];
MaxRangeDisplacement = Abs[dxOut2[[1]]];
}
];

MaxRangeCase4 = {MaxRangeDisplacement, MaxRangeAngle}
Case4Result = {ResolutionAtHWS, MaxRangeCase4}

Nominal Values for the Actuators: CASE 5 (Position 2)

Case 5 : (HWSY STEER M11 and HWSY STEER M12)
dM1dM2 = 1.357
dL1 = 0.7

Case5 = {dM1M2-> 0.625+0.05+0.1+0.582, dL1-> 0.7}

MOutCase = MOut/.Case5;
InvMout = Inverse[MOutCase]
Det[MOutCase]

Work Out Resolution

Resolution of angles & displacements on HWS
deltaAlphaMin = 0.7 urad
resolution at HWS = Sqrt[C2A^2 + C1A^2]*0.7 urad

C2A = DdAlpha2/. Case5
C1A =DdAlpha1/. Case5
resol = 0.7*10^-6;
ResolutionAtHWS = Sqrt[C1A^2 + C2A^2]*(resol)

mrad = 0.001;
RangeAct = 70 mrad;
RangeAtHWS = (Abs[C2A] + Abs[C1A])*RangeAct

Work Out Ranges

Range of angles & displacements on HWS
deltaAlphaMax = 70 umrad
Range at HWS = (Abs[C2A] + Abs[C1A])*70 mrad

mrad = 0.001;
RangeAct = 70 mrad;
RangeAtHWS = (Abs[C2A] + Abs[C1A])*RangeAct

Angular Range at HWS keeping displacement fixed
dHWS = 0;
Solve dx == 0 and substitute in solution for angles

eqSolve1 = deltaxHWS == MOutCase[[1]].{dAlpha2, dAlpha1};
soln1 = Solve[eqSolve1, dAlpha1][[1]];
angleOut =({MOutCase[[2]].{dAlpha2, dAlpha1}, dAlpha1}/. soln1)/. deltaxHWS-> 0;
angleOut1=angleOut /. dAlpha2-> RangeAct;
Print["Maximum range on Alpha 2 requires " <> ToString[100*Abs[angleOut1[[2]]]/RangeAct]<> "% on Alpha 1"];

soln2 = Solve[eqSolve1, dAlpha2][[1]];
angleOut =({MOutCase[[2]].{dAlpha2, dAlpha1}, dAlpha2}/. soln2)/. deltaxHWS-> 0;
angleOut2=angleOut /. dAlpha1-> RangeAct;
Print["Maximum range on Alpha 1 requires " <> ToString[100*Abs[angleOut2[[2]]]/RangeAct]<> "% on Alpha 2"];
Print[" "];

If[Abs[angleOut2[[2]]]>Abs[angleOut1[[2]]], {
Print["Maximum range of Angle 1 saturates Angle 2"];Print["Angular range is limited by Actuator 2"];
Print["Maximum range (keeping displacement on HWS constant) = "<> ToString[angleOut1[[1]]]<> " radians"];
MaxRangeAngle = Abs[angleOut1[[1]]];
} , {
Print["Maximum range of Angle 2 saturates Angle 1"];
Print["Angular range is limited by Actuator 1"];
Print["Maximum range (keeping displacement on HWS constant) = "<> ToString[angleOut2[[1]]]<> " radians"];
MaxRangeAngle = Abs[angleOut2[[1]]];
}
];

Displacement Range at HWS keeping angle fixed
Solve dtheta == 0 and substitute in solution for angles

eqSolve1 = deltaThetaHWS == MOutCase[[2]].{dAlpha2, dAlpha1};
soln1 = Solve[eqSolve1, dAlpha1][[1]];
dxOut =({MOutCase[[1]].{dAlpha2, dAlpha1}, dAlpha1}/. soln1)/. deltaThetaHWS-> 0;
dxOut1=dxOut /. dAlpha2-> RangeAct;
Print["Maximum range on Alpha 2 requires " <> ToString[100*Abs[dxOut1[[2]]]/RangeAct]<> "% on Alpha 1"];

soln2 = Solve[eqSolve1, dAlpha2][[1]];
dxOut =({MOutCase[[1]].{dAlpha2, dAlpha1}, dAlpha2}/. soln2)/. deltaThetaHWS-> 0;
dxOut2=dxOut /. dAlpha1-> RangeAct;
Print["Maximum range on Alpha 1 requires " <> ToString[100*Abs[dxOut2[[2]]]/RangeAct]<> "% on Alpha 2"];
Print[" "];

If[Abs[dxOut2[[2]]]>Abs[dxOut1[[1]]], {
Print["Maximum range of Angle 1 saturates Angle 2"];Print["Angular range is limited by Actuator 2"];
Print["Maximum range (keeping angle on HWS constant) = "<> ToString[dxOut1[[1]]*1000]<> " mm"];
MaxRangeDisplacement = Abs[dxOut1[[1]]];
} , {
Print["Maximum range of Angle 2 saturates Angle 1"];
Print["Angular range is limited by Actuator 1"];
Print["Maximum range (keeping angle on HWS constant) = "<> ToString[dxOut2[[1]]*1000]<> " mm"];
MaxRangeDisplacement = Abs[dxOut2[[1]]];
}
];

MaxRangeCase5 = {MaxRangeDisplacement, MaxRangeAngle}
Case5Result = {ResolutionAtHWS, MaxRangeCase5}

Position 3 - Around F1

f1 = .;
f2 = .;

MatrixListFromITMAlign = {(({{1, z3-dL2}, {0, 1}})),(({{1, dL2}, {0, 1}})).(({{1, 0}, {-1/f2, 1}})).(({{1, dL1}, {0, 1}})),(({{1, z2-dL1}, {0, 1}})). (({{1, 0}, {-1/f1, 1}})).(({{1, z1 }, {0, 1}})).(({{1, 0}, {0, (nFusedSilica840/n1)}})).(({{1, SR2RRT[[3]]}, {0, 1}})).(({{1, 0}, {((n1 - nFusedSilica840)/(SR2RoC *nFusedSilica840)), (n1/nFusedSilica840)}})).(({{1, SR2toSR3}, {0, 1}})).(({{1, 0}, {((-2)/SR3RoC), 1}})).(({{1, SR3toITMY}, {0, 1}})).(({{1, 0}, {((-2)/(-ITMRoC)), 1}}))}/. TestSolutionSubstitution;
nElemM = Length[MatrixListFromITMAlign]

Actuator 2 is closer to Viewport. Actuator 1 is closer to F1

InputBeam = {yIN, alphaIN};

BeamAtActuator2 = MatrixListFromITMAlign[[nElemM]].InputBeam;
deltaActuator2 = {0, dAlpha2};
BeamAfterActuator2 = Simplify[BeamAtActuator2 + deltaActuator2];

BeamAtActuator1= MatrixListFromITMAlign[[nElemM-1]].BeamAfterActuator2;
deltaActuator1 = {0, dAlpha1};
BeamAfterActuator1 =  Simplify[BeamAtActuator1 + deltaActuator1];

BeamAtHWS = Simplify[MatrixListFromITMAlign[[nElemM-2]].BeamAfterActuator1]

BeamAtHWS1 = Chop[Simplify[BeamAtHWS /. alphaIN-> 0]]

Get the Transfer Matrix Between the Actuators and the HWS

DdAlpha2 = Simplify[D[BeamAtHWS1, dAlpha2]]
DdAlpha1 = Simplify[D[BeamAtHWS1, dAlpha1]]

- If the same angle is added to both actuators but the magnitude is reversed then the displacement at the HWS is just equal to the angle times the distance between the actuators (no magnification?)

Sensing Matrix

MOut = {{DdAlpha2[[1]], DdAlpha1[[1]]}, {DdAlpha2[[2]], DdAlpha1[[2]]}}

Chop[Det[MOut]]

Nominal Values for the Actuators: CASE 6 (Position 3)

Case 6: (HWSY STEER M11 and HWSY STEER M12)
dL1 = 0.83564
dL2 = 0.700

Case6 = {dL1-> 0.835, dL2-> 0.7}

MOutCase = MOut/.Case6;
InvMout = Inverse[MOutCase]

Work Out Resolution

Resolution of angles & displacements on HWS
deltaAlphaMin = 0.7 urad
resolution at HWS = Sqrt[C2A^2 + C1A^2]*0.7 urad

C2A = DdAlpha2/. Case6;
C1A =DdAlpha1/. Case6;
resol = 0.7*10^-6;
ResolutionAtHWS = Sqrt[C1A^2 + C2A^2]*(resol)

mrad = 0.001;
RangeAct = 70 mrad;
RangeAtHWS = (Abs[C2A] + Abs[C1A])*RangeAct

Work Out Ranges

Range of angles & displacements on HWS
deltaAlphaMax = 70 umrad
Range at HWS = (Abs[C2A] + Abs[C1A])*70 mrad

mrad = 0.001;
RangeAct = 70 mrad;
RangeAtHWS = (Abs[C2A] + Abs[C1A])*RangeAct

Angular Range at HWS keeping displacement fixed
dHWS = 0;
Solve dx == 0 and substitute in solution for angles

eqSolve1 = deltaxHWS == MOutCase[[1]].{dAlpha2, dAlpha1};
soln1 = Solve[eqSolve1, dAlpha1][[1]];
angleOut =({MOutCase[[2]].{dAlpha2, dAlpha1}, dAlpha1}/. soln1)/. deltaxHWS-> 0;
angleOut1=angleOut /. dAlpha2-> RangeAct;
Print["Maximum range on Alpha 2 requires " <> ToString[100*Abs[angleOut1[[2]]]/RangeAct]<> "% on Alpha 1"];

soln2 = Solve[eqSolve1, dAlpha2][[1]];
angleOut =({MOutCase[[2]].{dAlpha2, dAlpha1}, dAlpha2}/. soln2)/. deltaxHWS-> 0;
angleOut2=angleOut /. dAlpha1-> RangeAct;
Print["Maximum range on Alpha 1 requires " <> ToString[100*Abs[angleOut2[[2]]]/RangeAct]<> "% on Alpha 2"];
Print[" "];

If[Abs[angleOut2[[2]]]>Abs[angleOut1[[2]]], {
Print["Maximum range of Angle 1 saturates Angle 2"];Print["Angular range is limited by Actuator 2"];
Print["Maximum range (keeping displacement on HWS constant) = "<> ToString[angleOut1[[1]]]<> " radians"];
MaxRangeAngle = Abs[angleOut1[[1]]];
} , {
Print["Maximum range of Angle 2 saturates Angle 1"];
Print["Angular range is limited by Actuator 1"];
Print["Maximum range (keeping displacement on HWS constant) = "<> ToString[angleOut2[[1]]]<> " radians"];
MaxRangeAngle = Abs[angleOut2[[1]]];
}
];

Displacement Range at HWS keeping angle fixed
Solve dtheta == 0 and substitute in solution for angles

eqSolve1 = deltaThetaHWS == MOutCase[[2]].{dAlpha2, dAlpha1};
soln1 = Solve[eqSolve1, dAlpha1][[1]];
dxOut =({MOutCase[[1]].{dAlpha2, dAlpha1}, dAlpha1}/. soln1)/. deltaThetaHWS-> 0;
dxOut1=dxOut /. dAlpha2-> RangeAct;
Print["Maximum range on Alpha 2 requires " <> ToString[100*Abs[dxOut1[[2]]]/RangeAct]<> "% on Alpha 1"];

soln2 = Solve[eqSolve1, dAlpha2][[1]];
dxOut =({MOutCase[[1]].{dAlpha2, dAlpha1}, dAlpha2}/. soln2)/. deltaThetaHWS-> 0;
dxOut2=dxOut /. dAlpha1-> RangeAct;
Print["Maximum range on Alpha 1 requires " <> ToString[100*Abs[dxOut2[[2]]]/RangeAct]<> "% on Alpha 2"];
Print[" "];

If[Abs[dxOut2[[2]]]>Abs[dxOut1[[1]]], {
Print["Maximum range of Angle 1 saturates Angle 2"];Print["Angular range is limited by Actuator 2"];
Print["Maximum range (keeping angle on HWS constant) = "<> ToString[dxOut1[[1]]*1000]<> " mm"];
MaxRangeDisplacement = Abs[dxOut1[[1]]];
} , {
Print["Maximum range of Angle 2 saturates Angle 1"];
Print["Angular range is limited by Actuator 1"];
Print["Maximum range (keeping angle on HWS constant) = "<> ToString[dxOut2[[1]]*1000]<> " mm"];
MaxRangeDisplacement = Abs[dxOut2[[1]]];
}
];

MaxRangeCase6 = {MaxRangeDisplacement, MaxRangeAngle}
Case6Result = {ResolutionAtHWS, MaxRangeCase3}

Results

Case1Result*1000
Case2Result*1000
Case3Result*1000
Case4Result*1000
Case5Result*1000
Case6Result*1000

4.5/17.5*180/Pi mrad

Final results

H1X

MatrixForm[H1XSoln]

Print["POtoVacLens = ",H1XSoln[[1]][[3]]]
Print["VacLens = ",H1XSoln[[1]][[5]]]
Print["VacLenstoAirLens = ",H1XSoln[[1]][[2]]]
Print["AirLens = ",H1XSoln[[1]][[4]]]
Print["AirLenstoHWS = ",H1XSoln[[1]][[1]]]

(({{1, 2.706}, {0, 1}})).(({{1, 0}, {-1/1.6994, 1}})).(({{1, 2.808}, {0, 1}})).(({{1, 0}, {1/1.6994, 1}})).(({{1, 2.156}, {0, 1}})).(({{1, 15.247}, {0, 1}})).(({{1, 0}, {-1/18, 1}})).(({{1, 24.666}, {0, 1}}))

inputbeam = OutRW[ 152/1000, 1330,820*10^-9,(({{1, 15.247}, {0, 1}})).(({{1, 0}, {-1/18, 1}})).(({{1, 24.666}, {0, 1}}))];

OutRW[wIN_, RIN_, λIN_, ABCD_]

PlotTheBeamPropagation[ABCDList_, RWInputMode_, λIN_]

DistancesFromSR3 = {15241.4,16498.9,16901.3,17067.6,17402.1,18992.4,19680.8,20210.2,20534.3,20834.3,20884.3,21982.5,22782.5,22916.1}/1000;
DistancesFromPO = DistancesFromSR3- DistancesFromSR3[[1]];

ABCDList =Reverse[ {(({{1, 2.706}, {0, 1}})),(({{1, 0}, {-1/1.6994, 1}})),(({{1, 2.808}, {0, 1}})),(({{1, 0}, {1/1.6994, 1}})),(({{1, 2.156}, {0, 1}}))}];

LambdaIn = 820*10^-9;

For[ii=1, ii≤ Length[DistancesFromPO], ii++, {
If[ii==1, plotL = ListPlot[{{DistancesFromPO[[ii]],0}, { DistancesFromPO[[ii]], 30}}, Joined-> True, PlotStyle-> Green], {
plotL1 =  ListPlot[{{DistancesFromPO[[ii]], 0}, { DistancesFromPO[[ii]], 30}}, Joined-> True, PlotStyle-> Green];
plotL = Show[plotL, plotL1];
}];
}];

plot1 = PlotTheBeamPropagation[ABCDList, inputbeam,LambdaIn];
Show[plot1, plotL, GridLines-> {None, Automatic}, PlotRange-> {0, 25}]
Export["/Advanced_LIGO/TCS/Hartmann sensor/H1_HWSX_beam_size.txt", BStempList, "Table"];
Export["/Advanced_LIGO/TCS/Hartmann sensor/H1_HWSX_beam_size_mirror_posns.txt", DistancesFromPO, "Table"];

Picomotors

LPM

MLPM = (({{1, 2.808-509.6/1000}, {0, 1}})).(({{1, 0}, {1/1.6994, 1}})).(({{1, 2.156}, {0, 1}})).(({{1, 15.247}, {0, 1}})).(({{1, 0}, {-1/18, 1}})).(({{1, 24.666}, {0, 1}}));
MLPMRev = Inverse[MLPM]

DLPM = MLPMRev.{0, alpha1}

UPM

MUPM = (({{1, 2.808-509.6/1000 - 630/1000}, {0, 1}})).(({{1, 0}, {1/1.6994, 1}})).(({{1, 2.156}, {0, 1}})).(({{1, 15.247}, {0, 1}})).(({{1, 0}, {-1/18, 1}})).(({{1, 24.666}, {0, 1}}))
MUPMRev = Inverse[MUPM];
DUPM = MUPMRev.{0, alpha2}

Hold alpha constant

Solve[DLPM[[2]] + DUPM[[2]]==0, alpha1]
DLPM[[1]] + DUPM[[1]]/. %[[1]]

H1Y

MatrixForm[H1YSoln]

Print["SR2ARtoVacLens = ",H1YSoln[[3]]]
Print["VacLens = ",H1YSoln[[5]]]
Print["VacLenstoAirLens = ",H1YSoln[[2]]]
Print["AirLens = ",H1YSoln[[4]]]
Print["AirLenstoHWS = ",H1YSoln[[1]]]

(({{1, 4.005}, {0, 1}})).(({{1, 0}, {-1/2.2659, 1}})).(({{1, 3.640}, {0, 1}})).(({{1, 0}, {1/1.6994, 1}})).(({{1, 2.333}, {0, 1}})).(({{1, 0}, {0, 1.453}})).(({{1, 0.075}, {0, 1}})).(({{1, 0}, {-0.04849, 0.6882}})).(({{1, 15.461}, {0, 1}})).(({{1, 0}, {-1/18, 1}})).(({{1, 24.792}, {0, 1}}));
inputbeam = OutRW[ 127/1000, 1330,820*10^-9,(({{1, 0}, {0, 1.453}})).(({{1, 0.075}, {0, 1}})).(({{1, 0}, {-0.04849, 0.6882}})).(({{1, 15.461}, {0, 1}})).(({{1, 0}, {-1/18, 1}})).(({{1, 24.792}, {0, 1}}))]

DistancesFromSR2 = {0,893.5,2095.7,2170.7,2333.0,2855.1,4334.2,5020.3,5973.0,6267.8,6549.4,6674.4,7756.0,8281.0,8478.1,8794.0}/1000.0;

ABCDList =Reverse[ {(({{1, H1YSoln[[1]]}, {0, 1}})),(({{1, 0}, {-1/H1YSoln[[4]], 1}})),(({{1, H1YSoln[[2]]}, {0, 1}})),(({{1, 0}, {-1/H1YSoln[[5]], 1}})),(({{1, H1YSoln[[3]]}, {0, 1}}))}];
LambdaIn = 820*10^-9;

For[ii=1, ii≤ Length[DistancesFromSR2], ii++, {
If[ii==1, plotL = ListPlot[{{DistancesFromSR2[[ii]],0}, { DistancesFromSR2[[ii]], 30}}, Joined-> True, PlotStyle-> Green], {
plotL1 =  ListPlot[{{DistancesFromSR2[[ii]], 0}, { DistancesFromSR2[[ii]], 30}}, Joined-> True, PlotStyle-> Green];
plotL = Show[plotL, plotL1];
}];
}];

plot1 = PlotTheBeamPropagation[ABCDList, inputbeam,LambdaIn];
Export["/Advanced_LIGO/TCS/Hartmann sensor/H1_HWSY_beam_size.txt", BStempList, "Table"];
Export["/Advanced_LIGO/TCS/Hartmann sensor/H1_HWSY_beam_size_mirror_posns.txt", DistancesFromSR2, "Table"];
Show[plot1, plotL, GridLines-> {None, Automatic}, PlotRange-> {0, 25}]

Periscope

UPM

MUPM = (({{1, 3.640 - (637 + 50+180+136.4+381.7+125+60)/1000}, {0, 1}})).(({{1, 0}, {1/1.6994, 1}})).(({{1, 2.333}, {0, 1}})).(({{1, 0}, {0, 1.453}})).(({{1, 0.075}, {0, 1}})).(({{1, 0}, {-0.04849, 0.6882}})).(({{1, 15.461}, {0, 1}})).(({{1, 0}, {-1/18, 1}})).(({{1, 24.792}, {0, 1}}));
MUPMRev = Inverse[MUPM]
MUPMRev.{0, alpha}

LPM

MLPM = (({{1, 3.640 - ( 50+180+136.4+381.7+125+60)/1000}, {0, 1}})).(({{1, 0}, {1/1.6994, 1}})).(({{1, 2.333}, {0, 1}})).(({{1, 0}, {0, 1.453}})).(({{1, 0.075}, {0, 1}})).(({{1, 0}, {-0.04849, 0.6882}})).(({{1, 15.461}, {0, 1}})).(({{1, 0}, {-1/18, 1}})).(({{1, 24.792}, {0, 1}}));
MLPMRev = Inverse[MLPM]
MLPMRev.{0, alpha}

Alignment Matrices

VacLenstoUPM = (637 + 50+180+136.4+381.7+125+60)/1000;
VacLenstoLPM = (50+180+136.4+381.7+125+60)/1000;
HWStoBS = 38.5*25.4/1000;
BStoLens1 = 18*25.4/1000;
BStoLens2 = 21*25.4/1000;
LensQPD = 200/1000;
LensQPDtoQPD1 = 9*25.4/1000;
LensQPDtoQPD2 = 5*25.4/1000;

alignMat1 =Simplify[ (({{1, LensQPDtoQPD}, {0, 1}})).(({{1, 0}, {-1/LensQPD, 1}})).(({{1, H1YSoln[[1]]-HWStoBS + BStoLens}, {0, 1}})).(({{1, 0}, {-1/H1YSoln[[4]], 1}})).(({{1, H1YSoln[[2]] - VacLenstoPM}, {0, 1}}))];

UPM

alignMatUPM1 = alignMat1/. {LensQPDtoQPD-> LensQPDtoQPD1, BStoLens-> BStoLens1, VacLenstoPM-> VacLenstoUPM};
alignMatUPM2 = alignMat1/. {LensQPDtoQPD-> LensQPDtoQPD2, BStoLens-> BStoLens2, VacLenstoPM-> VacLenstoUPM}; 
dxUPM1 =alignMatUPM1[[1]][[2]]
dxUPM2 =alignMatUPM2[[1]][[2]]

LPM

alignMatUPM1 = alignMat1/. {LensQPDtoQPD-> LensQPDtoQPD1, BStoLens-> BStoLens1, VacLenstoPM-> VacLenstoLPM};
alignMatUPM2 = alignMat1/. {LensQPDtoQPD-> LensQPDtoQPD2, BStoLens-> BStoLens2, VacLenstoPM-> VacLenstoLPM}; 
dxLPM1 =alignMatUPM1[[1]][[2]]
dxLPM2 =alignMatUPM2[[1]][[2]]

MatrixPeriscope

periM = {{dxUPM1, dxUPM2}, {dxLPM1, dxLPM2}}
periMInv = Inverse[periM]

MatrixForm[periM]
MatrixForm[periMInv]

Det[periM]

Det[IdentityMatrix[2]]