In the last post I presented the small module for gaussian optics that I’ve created for sympy. Here I’ll try to show a real problem solved with it and with sympy.
[update]: The api is still changing. Also here I’m using ABCD matrices for something that is more easily done with conjugation relations. It’s that way so I can show how the matrices are used in the module.
The problem is finding how a laser beam is modified by passing trough a lens. What is the new waist? Where is it obtained?
Making the boilerplate
Let us first introduce the variables for the input beam and the lens. Notice that as input for the creation of the BeamParameter the Rayleigh length \( z_r\) is used and not the waist \( w_0 = \sqrt{\frac{z_r \lambda}{\pi}}\). That’s the way for now just because it gives simpler expressions:
In [1]: from sympy.physics.gaussopt import *
In [2]:
In [3]: z_r = Symbol("z_r", positive=True)
In [4]: l = Symbol("lambda", positive=True)
In [5]: z = Symbol("z", real=True)
In [6]: f = Symbol("f", real=True)
In [7]: z_r
Out[7]: zᵣ
In [8]: l
Out[8]: λ
Defining the parameter and the lens and then calculating the parameter for the output beam:
In [9]: q_in = BeamParameter(z, z_r, l)
In [10]: tl = ThinLens(f)
In [11]:
In [12]: q_out= tl*q_in
Let’s see what does this give for the waist of the output parameter:
In [13]: w_0_out = q_out.w_0
In [14]: z_out = q_out.z
In [15]: w_0_out
Out[15]:
⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽
⎽⎽⎽ ⎽⎽⎽⎽ ╱ 1
╲╱ λ ⋅╲╱ zᵣ ⋅ ╱ ──────────────────
╱ 2 2
╱ 2⋅z z zᵣ
╱ 1 - ─── + ── + ───
╱ f 2 2
╲╱ f f
─────────────────────────────────────────
⎽⎽⎽
╲╱ π
In [16]: z_out
Out[16]:
2 2
z z zᵣ
────────────────── - ────────────────────── - ──────────────────────
2 2 ⎛ 2 2⎞ ⎛ 2 2⎞
2⋅z z zᵣ ⎜ 2⋅z z zᵣ ⎟ ⎜ 2⋅z z zᵣ ⎟
1 - ─── + ── + ─── f⋅⎜1 - ─── + ── + ───⎟ f⋅⎜1 - ─── + ── + ───⎟
f 2 2 ⎜ f 2 2⎟ ⎜ f 2 2⎟
f f ⎝ f f ⎠ ⎝ f f ⎠
In [17]: together(z_out) # to simplify it a bit
Out[17]:
⎛ 2 2⎞
f⋅⎝f⋅z - z - zᵣ ⎠
─────────────────────
2 2 2
f - 2⋅f⋅z + z + zᵣ
Those are still mostly unreadable. To simplify the expressions further I’ll introduce the dimensionless parameters \( Z_i = \frac{z}{f}\) and \( Z_r = \frac{z_r}{f}\). Also there is the obvious substitution for the input beam waist.
In [19]: Z_in = Symbol("Z_i", real=True)
In [20]: Z_r = Symbol("Z_r", real=True)
In [21]:
In [22]: Z_out = z_out/f
In [23]: Z_out = z_out.subs({z : Z_in*f, z_r : Z_r*f})
In [24]: Z_out = together(Z_out)
In [25]: Z_out
Out[25]:
⎛ 2 2⎞
f⋅⎝- Zᵢ + Zᵢ - Zᵣ ⎠
────────────────────
2 2
Zᵢ - 2⋅Zᵢ + Zᵣ + 1
And same for the waist:
In [26]: w_0_in = Symbol("w_0^i", positive=True)
In [27]: w_0_in
Out[27]: wⁱ₀
In [28]: w_0_in_expression = q_in.w_0
In [29]: w_0_in_expression
Out[29]:
⎽⎽⎽ ⎽⎽⎽⎽
╲╱ λ ⋅╲╱ zᵣ
────────────
⎽⎽⎽
╲╱ π
In [30]:
In [31]: w_0_out = w_0_out.subs({w_0_in_expression : w_0_in, z
: Z_in*f, z_r : Z_r*f})
In [32]: w_0_out
Out[32]:
wⁱ₀
─────────────────────────
⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽
╱ 2 2
╲╱ Zᵢ - 2⋅Zᵢ + Zᵣ + 1
Let’s study how the waist changes. First of all, there better be no solutions to the quadratic equations in the denominator. We want it always positive in order to have real waist. Let’s extract the expression.
In [41]: quad_equ = w_0_out.args[1].args[0]
In [42]: quad_equ
Out[42]:
2 2
Zᵢ - 2⋅Zᵢ + Zᵣ + 1
Ok, obviously the discriminant is $ -Z_r^2 \< 0 $ and we have no real solutions. But just as a crosscheck let’s pass it trough sympy’s solve.
In [44]: solve(quad_equ, Z_in)
Out[44]: [-ⅈ⋅│Zᵣ│ + 1, ⅈ⋅│Zᵣ│ + 1]
Great. The waist is real and positive.
Some useful results
Until now we haven’t done anything really interesting from physical point of view. We have just derived some expressions. Here are some more interesting questions.
What is the maximal possible waist obtainable with a single lens?
To find it we minimize the denominator. I lost a few minutes searching for a maximize function in sympy but didn’t find one. But obviously the maximum is at \( Z_i = 1\) and so \( w_{out}^{max} = \frac{w_{in}}{Z_r}=\frac{f \lambda}{\pi w_{in}}\).
It’s quite interesting that the maximum is for \( Z_i = 1\) (i.e. \( z_i = f\)). Let’s see what does that mean for a convergent lens (it’s mostly the same for a divergent one). The beam is converging to the focal point of the lens so after refraction the beam is parallel. And we know that in order to have parallel beam (i.e. minimize divergence) we must have big waist:
In [39]: q_in.divergence.subs({z_r : waist2rayleigh(w_0_in,l)})
Out[39]:
λ
─────
π⋅wⁱ₀
What about the minimal waist?
The predicted waist goes to zero for very big \( Z_i\). But the gauss beam is a solution to the paraxial equation which is invalid for great curvatures of the beam front. We have the following limit:
In [41]: q_out.waist_aproximation_limit
Out[41]:
2⋅λ
───
π
It would be interesting to know for what distance is this obtainable. The sad thing is that currently sympy can solve only very explicit equations. To rewrite it I do the following:
In [48]: waist_limit = q_out.waist_aproximation_limit
In [49]: equation = (w_0_out.args[0] / waist_limit)**2 -
(1/w_0_out.args[1])**2
In [50]: equation
Out[50]:
2 2
2 2 π ⋅wⁱ₀
- Zᵢ + 2⋅Zᵢ - Zᵣ - 1 + ───────
2
4⋅λ
In [51]:
In [52]: solve(equation, Z_in)
Out[52]:
⎡ ⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽ ⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎤
⎢ ╱ 2 2 2 2 ╱ 2 2 2 2 ⎥
⎢ ╲╱ - 4⋅Zᵣ ⋅λ + π ⋅wⁱ₀ ╲╱ - 4⋅Zᵣ ⋅λ + π ⋅wⁱ₀ ⎥
⎢1 - ─────────────────────────, 1 + ─────────────────────────⎥
⎣ 2⋅λ 2⋅λ ⎦
There are many more interesting questions but that’s all for now.