Tutorial 3: performing physical optics propagations.

In this third tutorial, we take the same optical system as in tutorial 2. However, instead of performing ray-traces through the system, we perform physical optics (PO).

%matplotlib widget 
# Uncomment for interactive plots when running the notebook!
# Requires installation of the ipympl package

import numpy as np

from PyPO.System import System
from PyPO.Enums import FieldComponents, Projections, Units, Scales

s = System()

2026-03-01 17:15:34 - WARNING - System override set to True. 

plane_focus = {
            "name"      : "plane_focus",
            "gmode"     : "uv",
            "lims_u"    : np.array([0, 0.1])*Units.MM,
            "lims_v"    : np.array([0, 360])*Units.DEG,
            "gridsize"  : np.array([101, 100])
            }

s.addPlane(plane_focus)

GPODict = {                                                                                                                                                                     
    "name"      : "focus",                                                                                                                                  
    "lam"       : 0.01*Units.MM,                                                                                                      
    "w0x"       : 0.05*Units.MM,                                                                                             
    "w0y"       : 0.05*Units.MM,                                                                                             
    "n"         : 1,                                                                                                                             
    "E0"        : 1,                                                                                                                                  
    "dxyz"      : 0,                                                                                                 
    "pol"       : np.array([1, 0, 0])                                                                                                          
}

s.createGaussian(GPODict, "plane_focus")
s.plotBeam2D("focus", FieldComponents.Ex, vmin=-30, vmax=0)

s.translateGrids("plane_focus", np.array([0, 0, 100])*Units.MM)

2026-03-01 17:15:34 - INFO - Added plane plane_focus to system. 

Figure

2026-03-01 17:15:35 - INFO - Translated element plane_focus by ('0.000e+00', '0.000e+00', '1.000e+02') millimeters. 

We start by defining a plane in the upper focus of the ellipsoid. We give it a radius of 0.1 mm. However, we do not place it in the focus just yet. First, we define a Gaussian PO beam on the plane. Because the Gaussian PO beams are always defined with their focus at x=y=z=0, we need to be careful to define this before we translate the focal plane.

The Gaussian beam itself is created by passing a filled dictionary, called a GPODict here, to the s.createGaussian() method. The dictionary consists of several fields. The "name" field contains the name by which the field will be stored in the internal fields dictionary of System. "lam" sets the wavelength, 0.01 mm in this case. After that, "w0x" and "w0y" set the beamwaist sizes along the x- and y-axis in mm, respectively. The refractive index of the medium in which the beam is defined is set by "n". The peak electric field value at the focus is set by "E0". The astigmatic z-axis distance between the focus of the beam in the x-plane and the y-plane in mm is set by "dxyz". Note that, if this option is used, the focus in the x-plane is placed at z=0 and the focus in the y-plane at z=-dxyz. The final field, "pol", sets the polarisation vector of the Gaussian beam.

ellipse = {
            "name"      : "ellipsoid",
            "pmode"     : "focus",
            "gmode"     : "uv",
            "flip"      : True,
            "focus_1"   : np.array([0, 0, 100])*Units.MM,
            "focus_2"   : np.array([0, 0, -100])*Units.MM,
            "orient"    : "z",
            "ecc"       : 0.5,
            "lims_u"    : np.array([0, 10])*Units.MM,
            "lims_v"    : np.array([0, 360]),
            "gridsize"  : np.array([11, 10])
            }

s.addEllipse(ellipse)
s.autoConverge("focus", "ellipsoid")

plane_t = {
            "name"      : "plane_t",
            "gmode"     : "uv",
            "lims_u"    : np.array([0, 5])*Units.MM,
            "lims_v"    : np.array([0, 360]),
            "gridsize"  : np.array([301, 300])
            }

s.addPlane(plane_t)
s.rotateGrids("plane_t", np.array([45, 0, 0]))

parabola = {
            "name"      : "paraboloid",
            "pmode"     : "focus",
            "gmode"     : "uv",
            "focus_1"   : np.array([0, -100, 0])*Units.MM,
            "vertex"    : np.array([-10, -100, 0])*Units.MM,
            "lims_u"    : np.array([0, .5])*Units.MM,
            "lims_v"    : np.array([0, 360]),
            "gcenter"   : np.array([0,-20])*Units.MM,
            "gridsize"  : np.array([301, 300])
            }

s.addParabola(parabola)

plane_out = {
            "name"      : "plane_out",
            "gmode"     : "xy",
            "lims_x"    : np.array([-2, 2])*Units.MM,
            "lims_y"    : np.array([-2, 2])*Units.MM,
            "gridsize"  : np.array([101, 101])
            }

s.addPlane(plane_out)
s.rotateGrids("plane_out", np.array([0, -90, 0]))
s.translateGrids("plane_out", np.array([50, -120, 0])*Units.MM)

2026-03-01 17:15:35 - INFO - Added ellipsoid ellipsoid to system. 
2026-03-01 17:15:35 - WORK - *** Starting auto-convergence ***  
2026-03-01 17:15:35 - WORK - Difference : 1e+99 at gridsize ('1.100e+01', '1.100e+01') 
2026-03-01 17:15:35 - WORK - Difference : 1.161640683809516 at gridsize ('2.100e+01', '2.100e+01') 
2026-03-01 17:15:35 - WORK - Difference : 0.3879551634916183 at gridsize ('3.100e+01', '3.100e+01') 
2026-03-01 17:15:35 - WORK - Difference : 0.19411849824491156 at gridsize ('4.100e+01', '4.100e+01') 
2026-03-01 17:15:35 - WORK - Difference : 0.11650743654744389 at gridsize ('5.100e+01', '5.100e+01') 
2026-03-01 17:15:35 - WORK - Difference : 0.07769657833499011 at gridsize ('6.100e+01', '6.100e+01') 
2026-03-01 17:15:35 - WORK - Difference : 0.05551009978387711 at gridsize ('7.100e+01', '7.100e+01') 
2026-03-01 17:15:36 - WORK - Difference : 0.041644036698322395 at gridsize ('8.100e+01', '8.100e+01') 
2026-03-01 17:15:36 - WORK - Difference : 0.032385592546248176 at gridsize ('9.100e+01', '9.100e+01') 
2026-03-01 17:15:36 - WORK - Difference : 0.025893077277380883 at gridsize ('1.010e+02', '1.010e+02') 
2026-03-01 17:15:36 - WORK - Difference : 0.02118800292914713 at gridsize ('1.110e+02', '1.110e+02') 
2026-03-01 17:15:36 - WORK - Difference : 0.017693320593757278 at gridsize ('1.210e+02', '1.210e+02') 
2026-03-01 17:15:36 - WORK - Difference : 0.014941571861857739 at gridsize ('1.310e+02', '1.310e+02') 
2026-03-01 17:15:36 - WORK - Difference : 0.012817017391348884 at gridsize ('1.410e+02', '1.410e+02') 
2026-03-01 17:15:36 - WORK - Difference : 0.011101457678936555 at gridsize ('1.510e+02', '1.510e+02') 
2026-03-01 17:15:37 - WORK - Difference : 0.009727354566013702 at gridsize ('1.610e+02', '1.610e+02') 
2026-03-01 17:15:37 - RESULT - Found converged solution, gridsize: ('1.449e+03', '1.449e+03') 
2026-03-01 17:15:37 - INFO - Added plane plane_t to system. 
2026-03-01 17:15:37 - INFO - Rotated element plane_t by ('4.500e+01', '0.000e+00', '0.000e+00') degrees around ('0.000e+00', '0.000e+00', '0.000e+00'). 
2026-03-01 17:15:37 - INFO - Added paraboloid paraboloid to system. 
2026-03-01 17:15:37 - INFO - Added plane plane_out to system. 
2026-03-01 17:15:37 - INFO - Rotated element plane_out by ('0.000e+00', '-9.000e+01', '0.000e+00') degrees around ('0.000e+00', '0.000e+00', '0.000e+00'). 
2026-03-01 17:15:37 - INFO - Translated element plane_out by ('5.000e+01', '-1.200e+02', '0.000e+00') millimeters. 

Here, we basically re-create the system from the previous tutorial. However, as we are doing PO now, the size of the reflectors become important. In the previous tutorial we purposefully oversized the off-axis paraboloid reflector for illustrative purposes. This time, we size the paraboloid in such a way that the illuminating beam has an edge taper between -10 and -15 dB.

Also, note that we put the gridsize of the ellipsoid at numpy.array([1,1]). This is, of course, way too low for doing a sensible propagation. However, also note we call the s.autoConverge() method with our Gaussian beam and ellipsoid as arguments. This method calculates the necessary gridsize on a surface in order for a converged solution to occur. In this case, it calculates the difference in power on the target as the gridsizes are increased. When this difference falls below a specified threshold ($10^{-2}$ in this case), the method finishes.

The method returns the gridsizes for reference, but also automatically adjusts the gridsize of the target surface, so that we do not have to manually adjust after the call.

device = 'GPU'

focus_to_ell_PO = {
    "t_name"      : "ellipsoid",
    "s_current"   : "focus",
    "mode"        : "JM",
    "name_JM"     : "JM_ell",
    "device"      : device
}

s.runPO(focus_to_ell_PO)

2026-03-01 17:15:37 - WORK - *** Starting PO propagation *** 
2026-03-01 17:15:37 - WORK - Propagating focus on plane_focus to ellipsoid, propagation mode: JM. 
2026-03-01 17:15:37 - WORK - Hardware: running 256 CUDA threads per block. 
2026-03-01 17:15:37 - WORK - ... Calculating ... 
2026-03-01 17:15:42 - WORK - *** Finished: 5.053 seconds *** 

<PyPO.PyPOTypes.currents at 0x7c3597063250>

ell_to_plane_t_PO = {
    "t_name"      : "plane_t",
    "s_current"   : "JM_ell",
    "mode"        : "JM",
    "name_JM"     : "JM_pt",
    "device"      : device
}

s.runPO(ell_to_plane_t_PO)

plane_t_to_par_PO = {
    "t_name"      : "paraboloid",
    "s_current"   : "JM_pt",
    "mode"        : "JMEH",
    "name_JM"     : "JM_par",
    "name_EH"     : "EH_par",
    "device"      : device
}

s.runPO(plane_t_to_par_PO)

2026-03-01 17:15:42 - WORK - *** Starting PO propagation *** 
2026-03-01 17:15:42 - WORK - Propagating JM_ell on ellipsoid to plane_t, propagation mode: JM. 
2026-03-01 17:15:42 - WORK - Hardware: running 256 CUDA threads per block. 
2026-03-01 17:15:42 - WORK - ... Calculating ... 
2026-03-01 17:16:54 - WORK - *** Finished: 72.352 seconds *** 
2026-03-01 17:16:54 - WORK - *** Starting PO propagation *** 
2026-03-01 17:16:54 - WORK - Propagating JM_pt on plane_t to paraboloid, propagation mode: JMEH. 
2026-03-01 17:16:54 - WORK - Hardware: running 256 CUDA threads per block. 
2026-03-01 17:16:54 - WORK - ... Calculating ... 
2026-03-01 17:16:56 - WORK - *** Finished: 1.658 seconds *** 

[<PyPO.PyPOTypes.currents at 0x7c367076cc30>,
 <PyPO.PyPOTypes.fields at 0x7c358ea06850>]

par_to_plane_out_PO = {
    "t_name"      : "plane_out",
    "s_current"   : "JM_par",
    "mode"        : "JMEH",
    "name_JM"     : "JM_out",
    "name_EH"     : "EH_out",
    "device"      : device
}

s.runPO(par_to_plane_out_PO)

2026-03-01 17:16:56 - WORK - *** Starting PO propagation *** 
2026-03-01 17:16:56 - WORK - Propagating JM_par on paraboloid to plane_out, propagation mode: JMEH. 
2026-03-01 17:16:56 - WORK - Hardware: running 256 CUDA threads per block. 
2026-03-01 17:16:56 - WORK - ... Calculating ... 
2026-03-01 17:16:56 - WORK - *** Finished: 0.384 seconds *** 

[<PyPO.PyPOTypes.currents at 0x7c367076cb00>,
 <PyPO.PyPOTypes.fields at 0x7c358ea06990>]

s.plotBeam2D("EH_out", FieldComponents.Ey, project=Projections.yz, vmin=-30, vmax=0)

Figure

We propagate the beam in the focus through the entire system onto the paraboloid reflector. We do not specify the "device" field in the PO dictionaries. This field defaults to "CPU" if the CUDA libraries are not compiled. If they are, the field defaults to "GPU". Similarly, the number of threads "nThreads" is not specified. If "device" = "CPU", "nThreads" defaults to the total number of CPU threads in your computer. If "device" = "GPU", "nThreads" defaults to 256 threads per CUDA block. In practice, this number is not critical and therefore does not have to be set often, but if set it should be a multiple of 32 (see this link) for a good overview about CUDA blocks and threads).

We set, for the first two propagations, the "mode" parameter to "JM". This means we only store the calculated electric (J) and magnetic (M) currents on the target surface. If we specify "EH", such as for the last propagation, we only save the illuminating electric (E) and magnetic (H) field on the target surface. If we want both, we specify "mode" as "JMEH". Another option, "FF" for far-field, will be explained in more detail below. The last option, "EHP", stores the reflected field and corresponding Poynting vectors. With this option it is possible to do a combined ray-trace and PO approach. This will be introduced in a later tutorial.

plane_ff = {
            "name"      : "plane_ff",
            "gmode"     : "AoE",
            "lims_Az"    : np.array([-0.7, 0.7]) * 6,
            "lims_El"    : np.array([-0.7, 0.7]) * 6 + 90,
            "gridsize"  : np.array([301, 301])
            }
s.addPlane(plane_ff)

par_to_plane_ff_PO = {
    "t_name"      : "plane_ff",
    "s_current"   : "JM_par",
    "mode"        : "FF",
    "name_EH"     : "EH_FF",
    "device"      : device
}

s.runPO(par_to_plane_ff_PO)

2026-03-01 17:16:57 - INFO - Added plane plane_ff to system. 
2026-03-01 17:16:57 - WORK - *** Starting PO propagation *** 
2026-03-01 17:16:57 - WORK - Propagating JM_par on paraboloid to plane_ff, propagation mode: FF. 
2026-03-01 17:16:57 - WORK - Hardware: running 256 CUDA threads per block. 
2026-03-01 17:16:57 - WORK - ... Calculating ... 
2026-03-01 17:16:58 - WORK - *** Finished: 1.312 seconds *** 

<PyPO.PyPOTypes.fields at 0x7c358ea41350>

s.plotBeam2D("EH_FF", FieldComponents.Ey, vmin=-30, scale=Scales.dB)

Units.MM
mm

Figure

In this final section, we propagated the field from the paraboloid to the far-field. We did this by specifying the "gmode" parameter to "AoE", which stands for Azimuth-over-Elevation (the only far-field co-ordinate system currently present in PyPO). The limits in "AoE" mode were given in degrees. Note that we subtracted 90 degrees from the Elevation limits. This is because the paraboloid illuminates the far-field along the x-axis. Also, note that the xyz components of the resulting far-field object do not align with the axes shown in the system plot. This is because, in the resulting fields object, x now aligns with the field along the Azimuthal direction (the y-axis in the original co-ordinate system) and y points along the Elevation direction (the z-axis in the original co-ordinate system).