Function Repository Resource:

GaussianBeam

Source Notebook

Represent a Gaussian beam and its parameters

Contributed by: Jaroslav Kysela

ResourceFunction["GaussianBeam"]["prop"]

returns a given property of the Gaussian beam.

ResourceFunction["GaussianBeam"][{"prop",{w0,zR}}]

specifies a beam waist w0 and Rayleigh range zR.

ResourceFunction["GaussianBeam"][{"prop","Normalized"}]

returns property prop for a unit beam waist and Rayleigh range.

ResourceFunction["GaussianBeam"][{"prop",<|k1v1,|>}]

specifies a given set of independent input parameters.

ResourceFunction["GaussianBeam"][{"prop",parspec,z}]

also specifies the longitudinal position.

ResourceFunction["GaussianBeam"][propspec,coords]

specifies the coordinate system.

ResourceFunction["GaussianBeam"][propspec,coords,False]

returns dependent parameters unevaluated.

Details

A Gaussian beam is an important type of optical beam that closely approximates different types of beams encountered in real life as well as various fields of physics.
ResourceFunction["GaussianBeam"] represents a Gaussian field that propagates along the vertical z-axis and that is parametrized by its beam waist w0 and Rayleigh range zR.
Properties are grouped into "Profiles", "Parameters", and "Summary", which groups one can retrieve by calling ResourceFunction["GaussianBeam"]["Properties"]. All available profiles and computed parameters are listed in ResourceFunction["GaussianBeam"]["Profiles"] and ResourceFunction["GaussianBeam"]["Parameters"], respectively.
Available profiles include:
"Amplitude"(real) amplitude of electromagnetic beam
"ComplexAmplitude"(complex) amplitude that includes both the real amplitude and phase
"Intensity"magnitude squared of amplitude
"Phase"phase of electromagnetic beam
Available computed parameters include:
"BeamRadius"beam radius at given longitudinal position
"ComplexBeamParameter"complex beam parameter
"Curvature"beam curvature
"Divergence"beam divergence
"FWHM"full width at half maximum
"GouyPhase"Gouy phase
"RadiusOfCurvature"radius of beam curvature
"RayleighRange"Rayleigh range
The formula for a given computed parameter par is obtained by calling ResourceFunction["GaussianBeam"][par]. For instance, the explicit formula for the beam divergence is retrieved by calling ResourceFunction["GaussianBeam"]["Divergence"].
Summary of all parameters and profiles is returned by ResourceFunction["GaussianBeam"]["Summary"] in the form of Dataset.
Not all beam parameters are independent of each other. The canonical set of independent parameters is the beam waist w0 and Rayleigh range zR, out of which all other properties can be calculated. The relation between the two independent parameters reads , where k=2πn/λ is the beam's wavenumber and where further n is the refractive index of the propagation medium and λ is the beam's wavelength. Instead of w0 and zR one could thus choose any complete set of these characteristics to specify the beam. For that reason, there are several ways of specifying a Gaussian beam:
{w0,zR}specific beam waist and Rayleigh range
"Normalized"equivalent to w0=1 and zR=1
<|k1v1,k2v2,|>Association with a set of independent parameters
When the parameters are specified as an Association, the following named keys are supported:
"BeamWaist"beam waist
"RayleighRange"Rayleigh range
"Wavenumber"wavenumber k=2πn/λ
"Wavelength"wavelength λ
"RefractiveIndex"refractive index n
There are also short keys supported for convenience:
"w0"beam waist
"zR"Rayleigh range
"k"wavenumber k=2πn/λ
"λ"wavelength λ
"n"refractive index n
There are two possibilities for an independent set of parameters: either 1) exactly two of w0, zR, and k are given; or 2) exactly one of w0 and zR together with both λ and n are given.
The input parameters can be specified in the following formats:
vsymbolic variable
0.325numeric value
Quantity[v,unit]symbolic or numeric value with unit
QuantityVariable[v,unit]named physical quantity with unit
Two coordinate systems are supported, which can be specified as the second argument:
"Cartesian"Cartesian coordinates {x,y,z} (default)
"Cylindrical"cylindrical coordinates {r,ϕ,z}
The Gaussian beam formula depends explicitly on several computed parameters such as beam radius and radius of curvature. These are by default automatically expanded into their functional representation. One can suppress this expansion by setting the third argument to False. This turns the parameters into QuantityVariable objects and the beam formula is thus returned in its higher-level form. Such a form might be beneficial when rendering the formula in some text or for reference purposes.

Examples

Basic Examples (4) 

Gaussian beam formula:

In[1]:=
ResourceFunction["GaussianBeam"][{"Amplitude", {w0, zR}}][x, y, z]
Out[1]=

Sample beam intensity along x-axis:

In[2]:=
xint = Table[
  ResourceFunction[
    "GaussianBeam"][{"Intensity", <|"BeamWaist" -> 0.7, "RayleighRange" -> 0.5|>}][x, 0, 0], {x, -1, 1, 0.1}]
Out[2]=
In[3]:=
ListPlot[xint]
Out[3]=

Parameter table for given inputs:

In[4]:=
ResourceFunction[
  "GaussianBeam"][{"Summary", <|"w0" -> Quantity[50, "Micrometers"], "n" -> 1., "\[Lambda]" -> Quantity[808, "Nanometers"]|>, Quantity[1, "Meters"]}]["Parameters"]
Out[4]=

Plot beam intersection along the optical axis:

In[5]:=
Plot[Evaluate[{-1, 1} ResourceFunction["GaussianBeam"][{"BeamRadius", {0.2, 0.3}}][
    z]], {z, -1, 1}, Sequence[
 PlotRange -> 1, PlotStyle -> Blue, Filling -> {1 -> {2}}, Epilog -> {
Arrow[{{0.2, 0.6}, {0., 0.2}}], 
Text["Beam waist", {0.2, 0.6}, {0, -1}]}]]
Out[5]=

Scope (15) 

Profiles (4) 

Available beam profiles:

In[6]:=
ResourceFunction["GaussianBeam"]["Profiles"]
Out[6]=

Beam intensity as a function of coordinates:

In[7]:=
ResourceFunction["GaussianBeam"][{"Intensity", {w0, zR}}]
Out[7]=

Plot it along z-axis:

In[8]:=
Plot[ResourceFunction["GaussianBeam"][{"Intensity", {1.2, 3.5}}][0, 0,
   z], {z, -5, 5}]
Out[8]=

Two-dimensional cross-section of beam amplitude for given z:

In[9]:=
zpos = 0.5;
beam = ResourceFunction["GaussianBeam"][{"Amplitude", {1, 1}, zpos}];
DensityPlot[beam[x, y], {x, -2, 2}, {y, -2, 2}]
Out[11]=

When no parameters are supplied, QuantityVariables are used instead:

In[12]:=
ResourceFunction["GaussianBeam"]["Intensity"]
Out[12]=
In[13]:=
InputForm[%]
Out[13]=

Coordinate systems (3) 

By default, Cartesian coordinate system is used:

In[14]:=
spec = {"ComplexAmplitude", {w0, zR}};
ResourceFunction["GaussianBeam"][spec][x, y, z] == ResourceFunction["GaussianBeam"][spec, "Cartesian"][x, y, z]
Out[15]=

Cylindrical coordinates are also supported:

In[16]:=
ResourceFunction["GaussianBeam"][spec, "Cylindrical"][r, \[Phi], z]
Out[16]=

The two systems are equivalent:

In[17]:=
ResourceFunction["GaussianBeam"][spec, "Cartesian"][x, y, z] == ResourceFunction["GaussianBeam"][spec, "Cylindrical"][##, z] & @@ ToPolarCoordinates[{x, y}] // Simplify
Out[17]=
In[18]:=
ResourceFunction["GaussianBeam"][spec, "Cylindrical"][r, \[Phi], z] ==
     ResourceFunction["GaussianBeam"][spec, "Cartesian"][##, z] & @@ FromPolarCoordinates[{r, \[Phi]}] // Simplify
Out[18]=

Computed parameters (4) 

Available computed parameters:

In[19]:=
ResourceFunction["GaussianBeam"]["Parameters"]
Out[19]=

Summary of all properties:

In[20]:=
Magnify[ResourceFunction["GaussianBeam"]["Summary"], 0.5]
Out[20]=

Beam radius formula:

In[21]:=
ResourceFunction["GaussianBeam"]["BeamRadius"]
Out[21]=

Beam radius as a function of longitudinal position:

In[22]:=
ResourceFunction["GaussianBeam"][{"BeamRadius", {w0, zR}}]
Out[22]=

Beam radius for a specific longitudinal position:

In[23]:=
ResourceFunction["GaussianBeam"][{"BeamRadius", {w0, zR}, 0.53}]
Out[23]=

Numerical parameter values are supported:

In[24]:=
ResourceFunction["GaussianBeam"][{"BeamRadius", {0.5, 1.2}, 0.53}]
Out[24]=

Syntax for most parameters follows that for beam radius:

In[25]:=
{
 ResourceFunction["GaussianBeam"]["FWHM"],
 ResourceFunction["GaussianBeam"][{"FWHM", {w0, zR}}],
 ResourceFunction["GaussianBeam"][{"FWHM", {0.2, 1}}][1.3],
 ResourceFunction["GaussianBeam"][{"FWHM", {0.2, 1}, 1.3}],
 ResourceFunction[
  "GaussianBeam"][{"FWHM", <|"BeamWaist" -> 0.2, "RayleighRange" -> 1|>,
    1.3}]
 }
Out[25]=

Input values (4) 

Input parameters can also have units:

In[26]:=
{
 ResourceFunction[
  "GaussianBeam"][{"Curvature", QuantityVariable["zR", "Distance"], QuantityVariable["z", "Distance"]}],
 ResourceFunction[
  "GaussianBeam"][{"Curvature", Quantity[zR, "Millimeters"], Quantity[z, "Millimeters"]}],
 ResourceFunction[
  "GaussianBeam"][{"Curvature", Quantity[1, "Millimeters"], Quantity[5.0, "Millimeters"]}]
 }
Out[26]=

Define dependent parameters by specifying beam waist w0 and Rayleigh range zR:

In[27]:=
{
 ResourceFunction["GaussianBeam"][{"BeamRadius", {w0, zR}, z}],
 ResourceFunction["GaussianBeam"][{"GouyPhase", zR, z}],
 ResourceFunction["GaussianBeam"][{"Divergence", {w0, zR}}]
 }
Out[27]=

Input the two values as Association:

In[28]:=
ResourceFunction[
 "GaussianBeam"][{"BeamRadius", <|"BeamWaist" -> w0, "RayleighRange" -> zR|>, z}]
Out[28]=

Short keys are supported:

In[29]:=
ResourceFunction[
 "GaussianBeam"][{"BeamRadius", <|"w0" -> w0, "zR" -> zR|>, z}]
Out[29]=

When input as Association, parameters other than beam waist and Rayleigh range can be supplied:

In[30]:=
ResourceFunction[
 "GaussianBeam"][{"BeamRadius", <|"w0" -> w0, "\[Lambda]" -> \[Lambda], "n" -> n|>, z}]
Out[30]=

Compare with:

In[31]:=
ResourceFunction["GaussianBeam"][{"BeamRadius", {w0, zR}, z}] /. zR -> ResourceFunction[
   "GaussianBeam"][{"RayleighRange", {w0, \[Lambda], n}}]
Out[31]=

Applications (2) 

Interactive calculator of Gaussian beam parameters:

In[32]:=
Manipulate[
 UnitConvert@Quiet@ResourceFunction["GaussianBeam"][{"Summary",
      Association["w0" -> Quantity[w0val, "Micrometers"], "n" -> nval,
        "\[Lambda]" -> Quantity[\[Lambda]val, "Nanometers"]]
      , Quantity[zval, "Meters"]}]["Parameters"],
 {{w0val, 50, "beam waist"}, 1, 100}, {{nval, 1, "refractive index"}, 0.1, 2},
 {{\[Lambda]val, 500, "wavelength"}, 100, 1000}, {{zval, 0.1, "longitudinal position"}, -5, 5}
 ]
Out[32]=

Do not evaluate dependent parameters and use default values of independent parameters for nice rendering:

In[33]:=
ResourceFunction[
 "GaussianBeam"]["ComplexAmplitude", "Cartesian", False]
Out[33]=
In[34]:=
TraditionalForm[%]
Out[34]=

Compare with evaluated formula:

In[35]:=
TraditionalForm@
 ResourceFunction["GaussianBeam"]["ComplexAmplitude", "Cartesian", True]
Out[35]=

Properties and Relations (4) 

Formulas are normalized such that the total intensity for a given longitudinal position is unity:

In[36]:=
intensity = ResourceFunction["GaussianBeam"][{"Intensity", {w0, zR}}][x, y, z]
Out[36]=
In[37]:=
totalIntensity = Integrate[
  intensity, {x, -\[Infinity], \[Infinity]}, {y, -\[Infinity], \[Infinity]}, Assumptions -> {w0 > 0, zR > 0, z \[Element] Reals}]
Out[37]=

Beam cross section is given by Gaussian (normal) distribution:

In[38]:=
intensity = ResourceFunction["GaussianBeam"][{"Intensity", {w0, zR}}][x, y, z];
gaussianDistr2D = PDF[NormalDistribution[0, \[Sigma]]][
    x] PDF[NormalDistribution[0, \[Sigma]]][y] /. \[Sigma] -> w0 /2 Sqrt[1 + z^2/zR^2]
Out[39]=
In[40]:=
Simplify[intensity == gaussianDistr2D]
Out[40]=

Intensity is the square of real amplitude:

In[41]:=
With[{i = ResourceFunction["GaussianBeam"][{"Intensity", {w0, zR}}][x, y, z],
   a = ResourceFunction["GaussianBeam"][{"Amplitude", {w0, zR}}][x, y,
     z]},
 Simplify[i == a^2]
 ]
Out[41]=

Complex amplitude contains both the amplitude and phase:

In[42]:=
With[{ca = ResourceFunction["GaussianBeam"][{"ComplexAmplitude", {w0, zR}}][x,
     y, z], a = ResourceFunction["GaussianBeam"][{"Amplitude", {w0, zR}}][x, y, z],
   p = ResourceFunction["GaussianBeam"][{"Phase", {w0, zR}}][x, y, z]},
 Simplify[ca == a Exp[I p]]
 ]
Out[42]=

Curvature is the reciprocal of the radius of curvature:

In[43]:=
With[{curv = ResourceFunction["GaussianBeam"][{"Curvature", {w0, zR}}][z], rad = ResourceFunction[
     "GaussianBeam"][{"RadiusOfCurvature", {w0, zR}}][z]},
 Simplify[1/curv == rad]
 ]
Out[43]=

For numerically heavy calculations with Gaussian profile use GaussianMatrix:

In[44]:=
{time1, mat1} = AbsoluteTiming[GaussianMatrix[20]];
{time2, mat2} = AbsoluteTiming[
   Table[ResourceFunction["GaussianBeam"][{"Intensity", {1, 1}, 0}][x,
      y], {x, -1, 1, 1/20.}, {y, -1, 1, 1/20.}]];
MatrixPlot /@ {mat1, mat2}
Out[46]=

Compare times:

In[47]:=
{time1, time2}
Out[47]=

Possible Issues (4) 

For convenience, radius of curvature at z=0 is set to infinity:

In[48]:=
ResourceFunction["GaussianBeam"][{"RadiusOfCurvature", zR}][0]
Out[48]=

Note however that limits from above and below zero differ:

In[49]:=
Limit[ResourceFunction["GaussianBeam"][{"RadiusOfCurvature", 1}][z], z -> 0, Direction -> "FromAbove"]
Out[49]=
In[50]:=
Limit[ResourceFunction["GaussianBeam"][{"RadiusOfCurvature", 1}][z], z -> 0, Direction -> "FromBelow"]
Out[50]=

Not all sets of input parameters define the beam uniquely:

In[51]:=
ResourceFunction[
 "GaussianBeam"][{"BeamRadius", <|"w0" -> w0, "\[Lambda]" -> \[Lambda]|>, z}]
Out[51]=

Supply all necessary parameters:

In[52]:=
ResourceFunction[
 "GaussianBeam"][{"BeamRadius", <|"w0" -> w0, "\[Lambda]" -> \[Lambda], "n" -> n|>, z}]
Out[52]=

Some sets of input parameters are overdetermined:

In[53]:=
ResourceFunction[
 "GaussianBeam"][{"BeamRadius", <|"w0" -> w0, "\[Lambda]" -> \[Lambda], "n" -> n, "k" -> kv|>, z}]
Out[53]=

Pick only the independent parameters:

In[54]:=
ResourceFunction[
 "GaussianBeam"][{"BeamRadius", <|"w0" -> w0, "n" -> n, "k" -> kv|>, z}]
Out[54]=

Assignment of two different values to the same parameter is not allowed:

In[55]:=
ResourceFunction[
 "GaussianBeam"][{"BeamRadius", <|"w0" -> w0, "BeamWaist" -> w0a, "n" -> n, "k" -> kv|>, z}]
Out[55]=

Pick only one of them:

In[56]:=
ResourceFunction[
 "GaussianBeam"][{"BeamRadius", <|"BeamWaist" -> w0a, "n" -> n, "k" -> kv|>, z}]
Out[56]=

Neat Examples (2) 

Plot high-intensity region around the waist of z-axis-aligned beam:

In[57]:=
intensity = ResourceFunction["GaussianBeam"][{"Intensity", "Normalized"}];
DensityPlot3D[
 intensity[x, y, z], {x, -15, 15}, {y, -15, 15}, {z, -15, 15}, PlotLabel -> "Absolute intensity", PlotPoints -> 100, Axes -> False]
Out[58]=

Plot intensity relative to the optical axis intensity:

In[59]:=
DensityPlot3D[
 intensity[x, y, z]/intensity[0, 0, z], {x, -15, 15}, {y, -15, 15}, {z, -15, 15}, PlotLabel -> "Relative intensity", Axes -> False]
Out[60]=

Plot wavefronts (might take a couple of minutes to compute):

In[61]:=
{time, wavefronts} = AbsoluteTiming[
   ContourPlot3D[
    ResourceFunction["GaussianBeam"][{"Phase", {1, 2}}][x, y, z], {x, -4, 4}, {y, -4, 4}, {z, -10, 10}, RegionFunction -> Function[{x, y, z, f}, ResourceFunction["GaussianBeam"][{"BeamRadius", {1, 2}}][z] > Sqrt[x^2 + y^2]], Sequence[
    Contours -> Subdivide[-32, 32, 10], ContourStyle -> Directive[
Specularity[20], 
Opacity[0.5, Purple]], Mesh -> None, RegionBoundaryStyle -> None]]];
time
Out[62]=
In[63]:=
volume = ResourceFunction["ExtractGraphicsPrimitives"]@
   RevolutionPlot3D[
    ResourceFunction["GaussianBeam"][{"BeamRadius", {1.1, 2}}][
     z], {z, -8, 8}, Sequence[
    RevolutionAxis -> {1, 0, 0}, Mesh -> None, PlotStyle -> Opacity[0.1, Blue], PlotPoints -> 20]];
Show[Graphics3D[{Rotate[
    volume, \[Pi]/2, {0, 1, 0}]}], wavefronts, Sequence[
 ViewPoint -> {1.75, 0.31, 0.9}, ViewVertical -> {0, 1, 0}, ImageSize -> 800, Boxed -> False, PlotRange -> {5 {-1, 1}, 5 {-1, 1}, 10 {-1, 1}}]]
Out[64]=

Publisher

Jaroslav Kysela

Requirements

Wolfram Language 13.0 (December 2021) or above

Version History

  • 1.0.0 – 19 April 2024

Related Resources

Author Notes

One could also add properties GaussianBeam["BeamWaist"] and GaussianBeam["Wavenumber"] akin to GaussianBeam["RayleighRange"].

License Information