ORIGINAL ARTICLE: Takatoshi Yamada, (2022), Practical numerical solutions of Wassermann–Wolf equations, Optik. 257: 168725. https://doi.org/10.1016/j.ijleo.2022.168725
Article abstract: In this paper, we present a complete method for determining the two surfaces of an aspheric singlet by numerical integration based on the classical Wassermann–Wolf equations (hereinafter WW). Although aplanatic optics is mainly treated in WW, we also show some other examples that can be determined by the WW. This indicates that surface shape determination by optimization software such as Zemax or CodeV is not necessary for some aspherical singlet designs, including the aplanatic case. A link to an executable Mathematica code for this purpose is also provided.

Aplanatic aspherical singlet design demonstration

The surface shapes of both sides are fixed by numerically solving the Wassermann-Wolf equation. The focal length always normalized to 1.0. You can specify some parameters, numerical aperture, lens thickness, CSF (Coddington Shape Factor that represents lens bending), wavelength and glass medium (extracted from SCHOTT catalog). Some sets of the parameters can not reach solutions with specified N.A. For these cases, this notebook shows maximum N.A. under the specified paraxial conditions.
​
Target N.A.
Thickness
CSF
Number of rays
Wavelength
Glass name
N-SF15

Contents

This notebook summarizes the derivation of expressions in the original paper,
“Practical numerical solutions of Wassermann–Wolf equations”:
​https://doi.org/10.1016/j.ijleo.2022.168725
This notebook contains derivations of :
1. Wassermann-Wolf equations from Snell’s law 2. explicit equations for infinite conjugate singlet 3. explicit equations for infinite conjugate singlet with a cover glass 4. explicit equations for finite conjugate singlet without cover glasses 5. explicit equations for finite conjugate singlet with cover glass in the image space 6. explicit equations for finite conjugate singlet with cover glass in the object space 7. explicit equations for finite conjugate singlet with cover glass in both spaces 8. explicit equations for afocal singlet 9. explicit equations for one-dimensional Gaussian-to-tophat converter 10. explicit equations for pseudo-Bessel beam generator
Please also refer the classical Wassermann–Wolf equations (hereinafter WW) paper :
Wassermann G.D., Wolf E.
“On the theory of aplanatic aspheric systems” Proc. Phys. Soc. B, 62 (1) (1949), p. 2

Notations

in WW, symbols concerning the second surface are specified by adding dash "'" but in Mathematica "'" symbol is used for derivative of Lagrange' s notation.
Therefore, subscript is used for assigning surfaces, e.g.,
In[]:=
{
ω
1
[t],
ω
s
[t],
ω
2
[t]};
mean
ω,
*
ω
,ω'
in WW.
And all Mathematica symbols defined by the author are started by small case characters because we can easily distinguish them from system symbols.
In Mathematica, dependent variables are better to be expressed explicitly with independent variable,e.g.,
In[]:=
ω[t];
not
In[]:=
ω;

Deriving Wassermann-Wolf equation

Preliminarily, Snell’s law in parametric surface expression

Here is well known Snell' s law expression.
Please put your insertion cursor inside the text box bellow and then input shift-return key to evaluate the expression.
After the evaluation, your Mathematica kernel has learned the expression for later use.
In[]:=
snellsLaw=(
n
1
*Sin[
θ
1
]
n
2
*Sin[
θ
2
]);
Please check the prompt header changed to “In[1]:=” from “In[*]:=” after evaluation.
Hereafter please evaluate text boxes with prompt “In[*]:=” in turn.
And please be careful that the order of evaluation may affect the results.
transforming it from the figure above,
In[]:=
replaceAngles={
θ
1

ω
1
-ϕ,
θ
2

ω
2
-ϕ};​​replaceSurfaceNormalAngle=(ϕArcTan[-x'[t]/y'[t]]);
leads the same mathematical meaning as an expression in Eq.(2.5) of WW. assuming y[t] is monotonous and concerning part of surface axis above
In[]:=
Simplify[TrigExpand[snellsLaw/.replaceAngles]/.replaceSurfaceNormalAngle,{Element[{x'[t]},Reals],y'[t]>0}]
Out[]=
n
1
(Cos[
ω
1
]
′
x
[t]+Sin[
ω
1
]
′
y
[t])
n
2
(Cos[
ω
2
]
′
x
[t]+Sin[
ω
2
]
′
y
[t])
The equation above is same meaning of Snell’s law expressed in Eq.(2.5) of WW, also same as Eq.(1) in the original paper.

Start from Snell's law for each surfaces, Eq.(2.5) in WW.

From here, dependent variables of 't' are expressed explicitly, i.e., ‘f[t]’ not ‘f’.
In[]:=
snell
1
=
n
1
(Cos[
ω
1
[t]]
x
1
'[t]+Sin[
ω
1
[t]]
y
1
'[t])n(Cos[
ω
s
[t]]
x
1
'[t]+Sin[
ω
s
[t]]
y
1
'[t]);​​
snell
2
=
n
2
(Cos[
ω
2
[t]]
x
2
'[t]+Sin[
ω
2
[t]]
y
2
'[t])n(Cos[
ω
s
[t]]
x
2
'[t]+Sin[
ω
s
[t]]
y
2
'[t]);
y dependency formulae, Eq.(2.8) and (2.9) in WW.
In[]:=
yDependencies[baseExpression]={
y
1
[t]
h
1
[t]+
x
1
[t]Tan[
ω
1
[t]],
y
2
[t]
h
2
[t]+
x
2
[t]Tan[
ω
2
[t]]};
and their derivatives
In[]:=
yDependencies[differentials]=Map[D[#,t]&,yDependencies[baseExpression]];
R definitions, Eq.(2.6) and Eq.(2.7) in WW.
In[]:=
omegaAstariscs={Sin[
ω
s
[t]]
R
y
/R,Cos[
ω
s
[t]]
R
x
/R};​​vectorRDefinition={
R
x
->
x
2
[t]-
x
1
[t]+d,
R
y

y
2
[t]-
y
1
[t],RSqrt[
R
x
^2+
R
y
^2]};
Rewrite expressions for derivative of x1 and x2 to get the equivalents to Eq.(2.10) and Eq.(2.11) in WW
In[]:=
wweq1011[original]=Thread[{
x
1
'[t],
x
2
'[t]}==Flatten[{
x
1
'[t],
x
2
'[t]}/.Simplify[Solve[{
snell
1
,
snell
2
}/.yDependencies[differentials]/.omegaAstariscs,{
x
1
'[t],
x
2
'[t]}]]]]
Out[]=

′
x
1
[t]-
(RSin[
ω
1
[t]]
n
1
-n
R
y
)(
′
h
1
[t]+
2
Sec[
ω
1
[t]]
x
1
[t]
′
ω
1
[t])
RSec[
ω
1
[t]]
n
1
-n(
R
x
+
R
y
Tan[
ω
1
[t]])
,
′
x
2
[t]-
(RSin[
ω
2
[t]]
n
2
-n
R
y
)(
′
h
2
[t]+
2
Sec[
ω
2
[t]]
x
2
[t]
′
ω
2
[t])
RSec[
ω
2
[t]]
n
2
-n(
R
x
+
R
y
Tan[
ω
2
[t]])


For confirmation of the equivalency, the first factor of Eq.(2.10) in WW.

before that, a preparatory helper function : fractionalRationalize[_,_]
In[]:=
Clear[fractionalRationalize];fractionalRationalize[fraction_,factor_]:=Apply[Divide,Map[Expand[#*factor]&,NumeratorDenominator[fraction]]]
here is a denominator factor of Eq.(2.10) and Eq.(2.11) in WW
In[]:=
originalFirstFactor=1/((n*
R
x
-R*Cos[ω[t]])/(n*
R
y
-R*Sin[ω[t]])+Tan[ω[t]])
Out[]=
1
-RCos[ω[t]]+n
R
x
-RSin[ω[t]]+n
R
y
+Tan[ω[t]]
In[]:=
Simplify[Map[Together,originalFirstFactor]]
and the second factor of Eq.(2.10)
lead the expression above.

Further calculation beyond Schwarzschild example in WW.

Initial conditions

and explicit expression of θ (t)
Hereafter, all singlets is in air.

Deriving equations for infinite conjugate singlet

The first example is an infinite conjugate lens.
It's an easy geometric exercise.
Eq.(2.10) and Eq.(2.11) in WW (equivalent expressions Eq(1) and Eq.(2) in the original) are derived only from Snell’s law, therefore the equations does not include conditions of stigmatism or aplanatism apparently.
Sine conditions or stigmatism need to formularize separately from WW equations as conditions for ray congruence.
The desired expressions for the infinite conjugate singlet can be obtained by substituting the conditions of the ray congruence into the equations.
transform to an explicit expression about variable t, eliminating dependencies with θ,ω etc...
The expressions are rather large and It is difficult to solve them analytically because it is nonlinear.

Trying to solve with numerical values by evaluating the equation above.

The standard function “NDSolve” can be used to solve this simultaneous nonlinear differential equations
And “ParametricPlot” function can be used to show the solved lens shape.

Infinite conjugate singlet with cover glass

Same as above except existence of a cover glass in the object space, right-side of a lens.

The ray congruence rules bellow include the refraction at the both surfaces of the cover glass.
transforming to an explicit expressions
evaluating the expressions numerically
and plot.

Finite Conjugate

Next example is a finite conjugate case.
Ray congruence conditions :
Then the numerical solution can be obtained :

Finite conjugate with cover glass in image space

Next example is a finite conjugate case with a cover glass.
evaluation of an expression bellow may be slow.

Finite conjugate with cover glass in object space

Next example is a finite conjugate case with a cover glass but it is in object space.
evaluation of an expression bellow may also be slow.

Finite conjugate with cover glass in both image and object space

Next example is a finite conjugate case with a cover glass in both space.
This configuration is commonly seen on optical disc pickups, i.e., an objective lens with a beam splitter and disc.
It may take very log time to evaluate the expression bellow because it will be very large.
But it can be numerically solved.

Afocal singlet

The next example is an afocal optics by aspherical singlet lens. It may only be used for beam expander, etc...?

All previous examples were completely aplastic singlets. However, the WW equations do not contain aplanatism in the equation itself, as mentioned earlier. Therefore, the equations and their methods above can be applied other optical design cases.
Combination of the equations and the ray congruence rules can also be solved numerically.

Tophat converter

An example of afocal system case.

Laser beam with gaussian profile is converted to tophat (flat - top) profile by the converter . This example is two-dimensional, in z-y plane. This case is not rotational symmetric but the equations and the methods can accept the case.
From two-dimensional geometrical intensity conservation,
leads a differential equation,

Bessel beam generator

The other non-aplanatic example.

Pseudo - Bessel beam can be produced by an axicon lens. This pseudo-Bessel beam generator example is better to get longer focal depth and working distance. The solution is rotational symmetric, not two-dimensional.
From three- dimensional geometrical intensity conservation, for small area ds,
where Ii(h) and Io(h2(h)) are intensities at height h respectively.
It leads a differential equation,
From total intensity conservation,the constant Υ0 can be determined,
And fix the constant of integration from the inside most ray,
Take positive solution,
Check the solution to lead marginal relation,
Angle in the image space ω1 be set to a constant κ.

Conclusion

WW equations and “NDSolve” function of Mathematica should have applications other than aplanatic lenses as shown in the afocal optics or the Bessel beam generator.

CITE THIS NOTEBOOK

Practical numerical solution of Wassermann - Wolf equations​
by Takatoshi Yamada@AYASE Co.,Ltd.​
Wolfram Community, STAFF PICKS, March 21, 2022
​https://community.wolfram.com/groups/-/m/t/2494155