ABCD analysis

by C. Durfee, CSM for PHGN480 Laser Physics

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: 
for free space, using only the definition that includes the refractive index

mFreeSpace[z_,n_]:=(({{1, z/n}, {0, 1}})); (* this is the same result as going through a window *)
mSimpleLens[f_]:=(({{1, 0}, {((-1)/f), 1}}));
mSphericalMirror[R_]:=(({{1, 0}, {((-2)/R), 1}}));
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]}}));

See file GaussianBeams ABCD.nb for more calculations relating to linear or nonlinear gradient index. 

additional matrices from Siegman
Tangential plane is where the beam width is measured along the plane of incidence (shortens focal length of tilted mirror).
Sagittal is the other one.

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

may have to fix to account for index inside medium

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];

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

define functions for Gaussian Beam ABCD

Summary of Gaussian beam ABCD functions

Summary of functions:
qinv[R, Z]	 	calculate inverse q from radius and Z=n π (w^2)/λ, where w is local beam size. beam is not collimated unless R=∞
z0[w, n] 		z0=n π (w^2)/λ this is Rayleigh range is beam is at a waist
beamWQ[n, qinv] 	calculate the local beam w from given inverse q. Have to specify local refractive index. 
beamRQ[qinv]		calculate the local beam R from given inverse q. Index not needed.

beamWz[z, w0, n]	beam radial width as function of distance z from waist of radius w0. This is standard definition for Gaussian beam. 
beamRz[z, w0, n]	beam radius of curvature as function of distance z from waist of radius w0, using standard definition for Gaussian

zWaist[w, r, n]		calculate position of waist given local w, R, and n
w0Waist[w, r, n]	calculate size of waist given local w, R, and n

qinvZ[z, n, qinvIn]	do ABCD propagation as function of z. gives new qinv
qinvOptic[z, qinvIn]	gives new qinv after going through an optical system specified by input ABCD 

zWaistQ[n, qinv]	calculate position of waist for given qinv
w0WaistQ[n, qinv]	calculate min waist radius for given qinv

stability[cavityABCD]	stability parameter for a round-trip cavity
StableModeSize[cavityABCD, n]
			stable mode radius at starting position in cavity. Permute matrices to get mode size at other positions.

plotBeamW[n, qinvIn, zStart, dz, zUnit, wUnit, plotColor]
plotBeamR[n, qinvIn, zStart, dz, zUnit, wUnit, plotColor]
			returns plot of the beam size (or beam radius of curvature) vs z, given initial qinv and index

calculations are done with inverse of complex beam parameter q
R is radius of wavefront curvature
Z is a Rayleigh range calculated from the local beam size (even though beam may not be collimated there)

definition of the q parameter is consistent with Svelto and Siegman

qinv[R_,Z_]:=(1/R)-I (1/Z);
z0[w_,n_]:=((n π (w^2))/lam0);

functions to extract beam info from qinv

beamWQ[n_,qinvIn_]:=(Sqrt[-(lam0/(n π Im[qinvIn]))]);
beamRQ[qinvIn_]:=If[Re[qinvIn]==0,∞,1/Re[qinvIn]];

beam prop equations based on distance from waist (z) and waist size (w0).

beamWz[z_,w0_,n_]:=w0(Sqrt[1+(z^2)/((z0[w0,n])^2)]);
beamRz[z_,w0_,n_]:=z((((((z0[w0,n])/z)))^2)+1);

these expressions give the size and location of the waist, given the size and radius of curvature at a given position

Clear[lam0]
Solve[{w1 == beamWz[z, w0, n], r0 == beamRz[z, w0, n]}, {z, w0}]

zWaist[w1_,r1_,n_]:=(((π^2) r1 (w1^4))/(((((lam0/n)))^2) (r1^2)+(π^2) (w1^4)));
w0Waist[w1_,r1_,n_]:=(((lam0/n) Abs[r1] w1)/(Sqrt[((((lam0/n)))^2) (r1^2)+(π^2) (w1^4)]));

propagation

for an optic: 
q=((A q + B)/(C q + D)), which is the same as: 
qinv=((C + D qinv)/(A + B qinv))
These functions are used to calculate the new qinv parameter either by propagating in z or by the application of an optical system ABCD.

qinvZ[z_,n_,qinvIn_]:=(((Subscript[mFreeSpace[z,n], [[2,1]]])+(Subscript[mFreeSpace[z,n], [[2,2]]]) qinvIn)/((Subscript[mFreeSpace[z,n], [[1,1]]])+(Subscript[mFreeSpace[z,n], [[1,2]]]) qinvIn)) ;
qinvOptic[qinvIn_,MatOptic_]:=(((Subscript[MatOptic, [[2,1]]])+(Subscript[MatOptic, [[2,2]]]) qinvIn)/((Subscript[MatOptic, [[1,1]]])+(Subscript[MatOptic, [[1,2]]]) qinvIn)) ;

zWaistQ[qinvIn], w0WaistQ[qinvIn] give the size and location of the waist, given the value of q at a given position
(q can be calculated from the known size and radius of curvature)

qinvZ[z, n, qReIn + I qImIn]

waist is at z where real part of this expression = 0

Expand[(i qImIn+qReIn)(1+(((-i qImIn+qReIn) z)/n))]

Solve[qReIn+(((qImIn^2) z)/n)+(((qReIn^2) z)/n)==0,z]

zWaistQ[n_,qinvIn_]:=-((n Re[qinvIn])/((Abs[qinvIn])^2));
w0WaistQ[n_,qinvIn_]:=beamWQ[n,qinvZ[zWaistQ[n,qinvIn],n,qinvIn]];

stability of an oscillator cavity

In Svelto: -1≤((A + D)/2)≤ 1, with edges of this zone being marginally stable. 
Using the parameter s, defined through
s=1-(((((A + D)/2)))^2), 
the cavity is stable if 0 ≤ s ≤ 1, and maximum stability is at s = 1, marginally stable at s = 0.

stability[MatCav_]:=1-((((((Subscript[MatCav, [[1,1]]])+(Subscript[MatCav, [[2,2]]]))/2)))^2);
StableModeSize[MatCavity_,n0_]:=If[stability[MatCavity]>0,(Sqrt[(lam0/(π n0))])(((Abs[(Subscript[MatCavity, [[1,2]]])])^(1/2))/((stability[MatCavity])^(1/4))),0];

plotting function: plot a propagation segment

plotBeamW[n_,qinvIn_,zStart_,dz_,zUnit_,wUnit_,pltCol_]:=Plot[beamWQ[n,qinvZ[(z zUnit-zStart),n,qinvIn]]/wUnit,{z,zStart/zUnit,(zStart+dz)/zUnit},PlotRange->All,PlotStyle->pltCol];
plotBeamR[n_,qinvIn_,zStart_,dz_,zUnit_,wUnit_,pltCol_]:=Plot[beamRQ[qinvZ[(z zUnit-zStart),n,qinvIn],1]/wUnit,{z,zStart/zUnit,(zStart+dz)/zUnit},PlotRange->Automatic,PlotStyle->pltCol];
plotBeamWGradient[γ_,qinvIn_,zStart_,dz_,zUnit_,wUnit_,pltCol_]:=Plot[beamWQ[1,qinvOptic[qinvIn,mGradientLens[γ,(z zUnit-zStart)]]]/wUnit,{z,zStart/zUnit,(zStart+dz)/zUnit},PlotRange->All,PlotStyle->pltCol];

examples

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]

ABCD propagation examples

propagate from a beam waist

lam0=0.589 um;w0=0.16 mm;R0=∞;col=Blue;
nTst=1;
plotBeamW[nTst,qinv[R0,z0[w0,nTst]],0,500 mm,mm,um,col]

focusing a beam

lam0=0.8 um;w0=0.5 mm;R0=∞;col=Blue;
n1=1;
Print["input beam Rayleigh range: ", z0[w0,n1]/m," m"];

fTst=50mm;
z2=150mm;
qInv0=qinv[R0,z0[w0,n1]];
qInvLens=qinvOptic[qInv0,mSimpleLens[fTst]];

plotBeamW[n1,qInvLens,0,z2,mm,um,col]

relay imaging: single lens

this example illustrates how to join propagation segments together into a single plot

lam0=0.8 um;w0=.041 mm;R0=∞;col=Blue;
n1=1;
Print["input beam Rayleigh range: ", z0[w0,n1]/m," m"];

fTst=50mm;
z1=70mm;
z2=((((1/fTst)-(1/z1)))^(-1));
Print["image plane: ",z2/mm," mm from lens"];
Print["magnification: ",mag=z2/z1];
qInv0=qinv[R0,z0[w0,n1]];
qInv1=qinvZ[z1,n1,qInv0];
qInvLens=qinvOptic[qInv1,mSimpleLens[fTst]];

p1=plotBeamW[n1,qInv0,0,z1,mm,um,col];
p2=plotBeamW[n1,qInvLens,z1,z2,mm,um,col];
Show[p1,p2,AxesOrigin->{0,0},GridLines->{{z1/mm},{w0/um,mag w0/um}}]

Put into a function to use Manipulate:

plotRelay[z1_]:=Block[{},
z2=((((1/fTst)-(1/z1)))^(-1));
mag=z2/z1;
qInv0=qinv[R0,z0[w0,n1]];
qInv1=qinvZ[z1,n1,qInv0];
qInvLens=qinvOptic[qInv1,mSimpleLens[fTst]];p1=plotBeamW[n1,qInv0,0,z1,mm,um,col];
p2=plotBeamW[n1,qInvLens,z1,z2,mm,um,col];
Show[p1,p2,AxesOrigin->{0,0},GridLines->{{z1/mm},{w0/um,mag w0/um}}]
];

Manipulate[plotRelay[z1 mm],{z1,51,200}]

Pick a given z1+z2 and look for focus at that plane. 
Evaluation of this block will change the Manipulate section above.

lentot=250mm;
plotRelay[z1_]:=Block[{},
z2=lentot-z1;
mag=z2/z1;
qInv0=qinv[R0,z0[w0,n1]];
qInv1=qinvZ[z1,n1,qInv0];
qInvLens=qinvOptic[qInv1,mSimpleLens[fTst]];p1=plotBeamW[n1,qInv0,0,z1,mm,um,col];
p2=plotBeamW[n1,qInvLens,z1,z2,mm,um,col];
Show[p1,p2,AxesOrigin->{0,0},GridLines->{{z1/mm},{w0/um,mag w0/um}}]
];

relay imaging: 2 lens 4f system

lam0=0.8 um;w0=0.05 mm;R0=10mm;col=Blue;
n1=1;
Print["input beam Rayleigh range: ", z0[w0,n1]/m," m"];

f1Tst=50mm;
f2Tst=100mm;
z1=f1Tst;
z2=f2Tst;
Print["magnification: ",mag=z2/z1];
qInv0=qinv[R0,z0[w0,n1]];
qInv1=qinvZ[z1,n1,qInv0];
qInvLens1=qinvOptic[qInv1,mSimpleLens[f1Tst]];
qInv3=qinvZ[f1Tst+f2Tst,n1,qInvLens1];
qInvLens2=qinvOptic[qInv3,mSimpleLens[f2Tst]];

p1=plotBeamW[n1,qInv0,0,z1,mm,um,col];
p2=plotBeamW[n1,qInvLens1,z1,f1Tst+f2Tst,mm,um,col];
p3=plotBeamW[n1,qInvLens2,z1+f1Tst+f2Tst,z2,mm,um,col];

Show[p1,p2,p3,AxesOrigin->{0,0},GridLines->{{z1/mm,(z1+f1Tst+f2Tst)/mm},{w0/um,mag w0/um}}]

ABCD analysis of cavities

2 mirror cavity

stability calculation for a 2-mirror cavity:

mSys=mSphericalMirror[R1].mFreeSpace[L,1].mSphericalMirror[R2].mFreeSpace[L,1];
MatrixForm[mSys]

calculate ((A + D)/2):

eq1=FullSimplify[(mSys[[2,2]]+mSys[[1,1]])/2]

define (Subscript[g, 1])=1-L/(Subscript[R, 1]) and (Subscript[g, 2])=1-L/(Subscript[R, 2])

eq2=eq1/.{R1->L/(1-g1),R2->L/(1-g2)};
FullSimplify[eq2]

since the stability condition is -1≤((A + D)/2)≤1, -1≤2 (Subscript[g, 1])(Subscript[g, 2])-1≤1, or  0≤(Subscript[g, 1])(Subscript[g, 2])≤1
In terms of the s parameter defined above, s=1-(((((A + D)/2)))^2)=1-(((2(Subscript[g, 1])(Subscript[g, 2]) -1))^2)=1-(((2(Subscript[g, 1])(Subscript[g, 2]) -1))^2)

FullSimplify[1-(((2 g1 g2-1))^2)]

stabilityg1g2[g1_,g2_]:=-4 g1 g2 (g1 g2-1);

Plot[-4 g1g2 (-1+g1g2),{g1g2,-.1,1.1}]

The following plot is a variation on the standard stability plot. The dark edges are at marginal stability, and the most stable regions are the light bands in the middle.

gMax=3;
pStability=ContourPlot[stabilityg1g2[g1,g2],{g1,-gMax,gMax},{g2,-gMax,gMax},PlotRange->{0,1}]

or, in terms of the cavity length for a given pair of mirrors
The sequence starts at the end mirror 1, so the stable mode size would be at that mirror

Clear[mSys];
mSys[R1_,R2_,L_]=mSphericalMirror[R1].mFreeSpace[L,1].mSphericalMirror[R2].mFreeSpace[L,1];

Example: cavity with unequal curved mirrors

lam0=632.8nm;
R1Tst=80 cm;R2Tst=100 cm;
Lmax=Min[R1Tst+R2Tst,1];
Lmax=2R2Tst+1cm;
Plot[{
StableModeSize[mSys[R1Tst,R2Tst,L],1]/mm,StableModeSize[mSys[R2Tst,R1Tst,L],1]/mm,
stability[mSys[R1Tst,R2Tst,L]]},{L,0,Lmax},
PlotRange->{0,1},PlotStyle->{Red,Blue,GrayLevel[0]}]

StableModeSize[mSys[R1Tst,R2Tst,10cm],1]

pLslice=ParametricPlot[{1-L/R1Tst,1-L/R2Tst},{L,0,Lmax},PlotStyle->Red];
Show[pStability,pLslice]