ABCD analysis

initializations

constants

$TextStyle={FontFamily->"Times",FontSize->14,FontWeight->"Bold"}

m=100;J=1;
cm=(10^(-2))m;mm=(10^(-3))m;um=(10^(-6))m;nm=(10^(-9))m;
uJ=(10^(-6))J;mJ=(10^(-3))J;
m=1;W=1;
cm=(10^(-2))m;mm=(10^(-3))m;um=(10^(-6))m;nm=(10^(-9))m;
J=1;uJ=(10^(-6))J;nJ=(10^(-9))J;mJ=(10^(-3))J;
s=1;fs=(10^(-15))s;
c=3 (10^8);
ϵ0=8.8 (10^(-12));

define abcd elements:

mFreeSpace[z_]:=(({{1, z}, {0, 1}}));
mSimpleLens[f_]:=(({{1, 0}, {((-1)/f), 1}}));
mSphericalMirror[R_]:=(({{1, 0}, {((-2)/R), 1}}));
mSimpleLensX[f_,θ_]:=(({{1, 0}, {((-1)/(f Cos[1.8 θ])), 1}}));(* tilt around y-axis, fx changes *)
mSimpleLensY[f_,θ_]:=(({{1, 0}, {((-1)/(f Cos[0.8 θ])), 1}}));(* tilt around x-axis, fy changes *)
mSphericalMirrorX[R_,θ_]:=(({{1, 0}, {((-2)/(R  Cos[θ])), 1}}));
mSphericalMirrorY[R_,θ_]:=(({{1, 0}, {((-2Cos[θ])/R), 1}}));
mGradientLens[γ_,dz_]:=(({{Cos[γ dz], Sin[γ dz]/γ}, {-γ * Sin[γ dz], Cos[γ dz]}}));

mFreeSpaceIndex[z_,n_]:=(({{1, z/n}, {0, 1}})); (* this is the same result as going through a window *)
mDielFlatInterface[n1_,n2_]:=(({{1, 0}, {0, (n1/n2)}}));

mFlatInterfaceTang[thInc_,n1_,n2_]:=(({{((Cos[ArcSin[n1 Sin[thInc]/n2]])/(Cos[thInc])), 0}, {0, ((Cos[thInc])/(Cos[ArcSin[n1 Sin[thInc]/n2]]))}}));
mFlatInterfaceSagittal[thInc_,n1_,n2_]:=(({{1, 0}, {0, 1}}));

functions for ABCD raytrace

plotting function: plot a propagation segment
y0 = rayIn[[1]] initial height
y' = rayIn[[2]] slope

line to plot: 
y=y0+(z - z0) y' = rayIn[[1]] + (z - z0)  rayIn[[2]]

plotRay[rayIn_,zStart_,dz_,zUnit_,wUnit_,pltCol_]:=Plot[(rayIn[[1]]+(z zUnit-zStart)rayIn[[2]])/wUnit,{z,zStart/zUnit,(zStart+dz)/zUnit},PlotRange->All,PlotStyle->pltCol];

test plotRay

y0=0.6 mm;yPrime=((0.1mm)/(10 mm));
col=Blue;
ray={y0,yPrime};
plotRay[ray,0,300 mm,mm,um,col]

test raytrace

single lens - collimated input

nSegments=2;
rayStart={0.5,0};
z1=10mm;z2=30mm;f1=20mm;
ray1=mSimpleLens[f1].mFreeSpace[z1].rayStart;
rayEnd=mFreeSpace[z2].mSimpleLens[f1].mFreeSpace[z1].rayStart;
Print[rayEnd]

pltLst=Range[nSegments];
pltLst[[1]]=plotRay[rayStart,0,z1,mm,mm,col];
pltLst[[2]]=plotRay[ray1,z1,z2,mm,mm,col];
Show[pltLst]

The following defines a collimated set of rays for the input to the system.

z1=10mm;z2=100mm;f1=20mm;

beamWid=2.;
nRays=10;
dy=beamWid/(nRays-1);
rayHts=Table[(i-1) dy-beamWid/2,{i,1,nRays}];
rayAngles=0 Range[nRays];
raySet0=Transpose[{rayHts,rayAngles}];

This defines a point source in a similar way - starting position at a certain height, then define a set of rays.

beamHt=0;
nRays=10;
angleRange=0.1;
dA=angleRange/(nRays-1);
rayHts=ConstantArray[beamHt,nRays];
rayAngles=Table[(i-1) dA-angleRange/2,{i,1,nRays}];
raySet0=Transpose[{rayHts,rayAngles}];

The following lines apply the matrix dot multiplication to each ray in the list in the raySets
the # is a placeholder for the function argument, the & identifies what precedes it as a "pure function", and the /@ performs the Map[ ] function that does the . product on each of the vectors in the list. 

To do a trace, you make a ray set for the beginning of each propagation segment. 
Here raySet0 is the beginning, propagating forward from z = 0 to z = z1
raySet1 results from forward propagation to the lens, then propagation through the lens. 
raySet2 is not needed for the display of the rays, but could be used to numerically inspect what the rays are like at a certain plane after the lens.

z1=10mm;

raySet1=(mSimpleLens[f1].mFreeSpace[z1].#)&/@raySet0;
raySet2=(mFreeSpace[z2].#)&/@raySet1;

pltLst=ConstantArray[0,{nSegments,nRays}];
pltLst[[1]]=plotRay[#,0,z1,mm,mm,col]&/@raySet0;
pltLst[[2]]=plotRay[#,z1,z2,mm,mm,col]&/@raySet1;
Show[pltLst,GridLines->{{z1/mm},None},GridLinesStyle->Red]