Demo 6: building and validating a Dragonian dual reflector¶

In this demo, we build a Dragonian reflector and validate its design by simulating the propagation of a Gaussian beam through the setup. Specifically, we propagate from one focus of the Dragonian to the other. We compare the beamwidths of the Gaussian at both focii and see if they match.

The Dragonian in this demo is an existing model, currently used for aligning the DESHIMA 2.0 spectrometer. It has a magnification factor of 1.89, and it is this factor that we will check now.

WARNING: this demo is quite heavy on the hardware. If you do not have a CUDA capable graphics card, this demo can take a long time to finish. If you do have a CUDA graphics card, this demo should be no problem to run.

In [1]:
#%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

# Ellipsoidal reflector: note that these parameters are for the 2D case.
_A_ELLIPSE  = 3689.3421*Units.MM / 2             # Semi-major axis in mm
_C_ELLIPSE  = 1836.4965*Units.MM / 2             # Focii distance in mm
_B_ELLIPSE  = 3199.769638*Units.MM / 2           # Semi-minor axis in mm

_X_LIM_ELL  = np.array([1435, 1545])*Units.MM    # In frame where lower ellipse vertex is at origin
_Y_LIM_ELL  = np.array([-200, 200])*Units.MM

# Hyperboloid reflector: note that these parameters are for the 2D case.
_A_HYPERBO  = 1226.5776*Units.MM / 2             # Vertex distance in mm
_C_HYPERBO  = 2535.878*Units.MM / 2              # Focii distance in mm
_B_HYPERBO  = 2219.500985*Units.MM / 2           # Semi-minor axis in mm

_X_LIM_HYP  = np.array([160, 480])*Units.MM
_Y_LIM_HYP  = np.array([-160, 160])*Units.MM

# First, define focii of hyperbola and ellipsoid
h_f0 = np.zeros(3) # Horizontal focus of Dragonian = lower focus of hyperboloid
h_f1 = np.array([2298.285978374926, 0.0, 1071.7083523464814])*Units.MM # Upper hyperboloid focus - 
                                                              # co-incides with far ellipsoid focus

e_f0 = np.array([486.3974883317985, 0.0, 1371.340617233771])*Units.MM # Upper focus of Dragonian = close ellipsoid focus
e_f1 = h_f1 # Far ellipsoid focus

e_center = (e_f1 + e_f0) / 2 
diff = (e_f1 - e_f0) / np.sqrt(np.dot(e_f1 - e_f0, e_f1 - e_f0))
theta = np.degrees(np.arccos(np.dot(np.array([1,0,0]), diff)))
#print(theta)

# Initialize system
s = System()

# Add parabolic reflector and hyperbolic reflector by focus, vertex and two foci and eccentricity
h_coeff = np.array([_B_HYPERBO, _B_HYPERBO, _A_HYPERBO])
h_gridsize = np.array([201, 201])

e_coeff = np.array([_B_ELLIPSE, _B_ELLIPSE, _A_ELLIPSE])
e_gridsize = np.array([201, 201])

h_wo = {
        "name"      : "h_wo",
        "pmode"     : "manual",
        "gmode"     : "xy",
        "coeffs"    : h_coeff,
        "flip"      : False,
        "lims_x"    : _X_LIM_HYP,
        "lims_y"    : _Y_LIM_HYP,
        "gridsize"  : h_gridsize
        }

e_wo = {
        "name"      : "e_wo",
        "pmode"     : "manual",
        "gmode"     : "xy",
        "coeffs"    : e_coeff,
        "flip"      : True,
        "lims_x"    : _X_LIM_ELL,
        "lims_y"    : _Y_LIM_ELL,
        "gridsize"  : e_gridsize
        }

s.addHyperbola(h_wo)
s.addEllipse(e_wo)

s.rotateGrids("h_wo", np.array([180, 0, 0]), pivot=np.zeros(3))
s.translateGrids("h_wo", np.array([0, 0, _C_HYPERBO]))
s.rotateGrids("h_wo", np.array([0, 65, 0]), pivot=np.zeros(3))

s.rotateGrids("e_wo", np.array([180, 90, 0]), pivot=np.zeros(3))
s.translateGrids("e_wo", e_center)
s.rotateGrids("e_wo", np.array([0, 1*theta, 0]), pivot=e_center)

# Build the plane in the horizontal focus for the source Gaussian beam.
focus1 = {
        "name"     : "focus1",
        "gmode"    : "xy",
        "lims_x"   : np.array([-20, 20])*Units.MM,
        "lims_y"   : np.array([-20, 20])*Units.MM,
        "gridsize" : np.array([101, 101])
        }

s.addPlane(focus1)

# Define Gaussian on focus1 before rotating
GPOfield =  {
        "name"      : "Gauss1",
        "lam"       : 1*Units.MM,
        "w0x"       : 10*Units.MM,
        "w0y"       : 10*Units.MM,
        "n"         : 1,
        "E0"        : 1, 
        "dxyz"      : 0,
        "pol"       : np.array([1,0,0])
}

s.createGaussian(GPOfield, "focus1")

# Rotate the plane with 90 degrees around the y axis so that it faces the hyperboloid mirror.
s.rotateGrids("focus1", np.array([0, 90, 0]))

s.autoConverge("Gauss1", "h_wo", tol=1e-3)

# Build the plane in the vertical focus for calculation.
focus2 = {
        "name"     : "focus2",
        "gmode"    : "xy",
        "lims_x"   : np.array([-40, 40])*Units.MM,
        "lims_y"   : np.array([-40, 40])*Units.MM,
        "gridsize" : np.array([101, 101])
        }

s.addPlane(focus2)
s.translateGrids("focus2", e_f0)

s.plotSystem()

2026-03-20 15:09:59 - WARNING - System override set to True. 
2026-03-20 15:09:59 - INFO - Added hyperboloid h_wo to system. 
2026-03-20 15:09:59 - INFO - Added ellipsoid e_wo to system. 
2026-03-20 15:09:59 - INFO - Rotated element h_wo by ('1.800e+02', '0.000e+00', '0.000e+00') degrees around ('0.000e+00', '0.000e+00', '0.000e+00'). 
2026-03-20 15:09:59 - INFO - Translated element h_wo by ('0.000e+00', '0.000e+00', '1.268e+03') millimeters. 
2026-03-20 15:09:59 - INFO - Rotated element h_wo by ('0.000e+00', '6.500e+01', '0.000e+00') degrees around ('0.000e+00', '0.000e+00', '0.000e+00'). 
2026-03-20 15:09:59 - INFO - Rotated element e_wo by ('1.800e+02', '9.000e+01', '0.000e+00') degrees around ('0.000e+00', '0.000e+00', '0.000e+00'). 
2026-03-20 15:09:59 - INFO - Translated element e_wo by ('1.392e+03', '0.000e+00', '1.222e+03') millimeters. 
2026-03-20 15:09:59 - INFO - Rotated element e_wo by ('0.000e+00', '9.390e+00', '0.000e+00') degrees around ('1.392e+03', '0.000e+00', '1.222e+03'). 
2026-03-20 15:09:59 - INFO - Added plane focus1 to system. 
2026-03-20 15:09:59 - INFO - Rotated element focus1 by ('0.000e+00', '9.000e+01', '0.000e+00') degrees around ('0.000e+00', '0.000e+00', '0.000e+00'). 
2026-03-20 15:09:59 - WORK - *** Starting auto-convergence ***  
2026-03-20 15:09:59 - WORK - Difference : 1e+99 at gridsize ('1.100e+01', '1.100e+01') 
2026-03-20 15:09:59 - WORK - Difference : 0.0690600616330952 at gridsize ('2.100e+01', '2.100e+01') 
2026-03-20 15:09:59 - WORK - Difference : 0.02416343647335717 at gridsize ('3.100e+01', '3.100e+01') 
2026-03-20 15:09:59 - WORK - Difference : 0.009342896868539396 at gridsize ('4.100e+01', '4.100e+01') 
2026-03-20 15:09:59 - WORK - Difference : 0.007032639637988725 at gridsize ('5.100e+01', '5.100e+01') 
2026-03-20 15:09:59 - WORK - Difference : 0.005043108895898429 at gridsize ('6.100e+01', '6.100e+01') 
2026-03-20 15:10:00 - WORK - Difference : 0.0028070152566510576 at gridsize ('7.100e+01', '7.100e+01') 
2026-03-20 15:10:00 - WORK - Difference : 0.0020824505995415965 at gridsize ('8.100e+01', '8.100e+01') 
2026-03-20 15:10:00 - WORK - Difference : 0.0016662599343750806 at gridsize ('9.100e+01', '9.100e+01') 
2026-03-20 15:10:00 - WORK - Difference : 0.0014319183183884032 at gridsize ('1.010e+02', '1.010e+02') 
2026-03-20 15:10:00 - WORK - Difference : 0.0012609546510740088 at gridsize ('1.110e+02', '1.110e+02') 
2026-03-20 15:10:00 - WORK - Difference : 0.001238709780188918 at gridsize ('1.210e+02', '1.210e+02') 
2026-03-20 15:10:00 - WORK - Difference : 0.0005955903267107487 at gridsize ('1.310e+02', '1.310e+02') 
2026-03-20 15:10:00 - RESULT - Found converged solution, gridsize: ('1.179e+03', '1.179e+03') 
2026-03-20 15:10:00 - INFO - Added plane focus2 to system. 
2026-03-20 15:10:00 - INFO - Translated element focus2 by ('4.864e+02', '0.000e+00', '1.371e+03') millimeters.
No description has been provided for this image

We plot the Dragonian reflector. We also plot the planes we defined in the two focii. Now we will start the calculations. WARNING! If not using CUDA, this can take a long time!

In [2]:
device="GPU"

runPODict1 = {
            "t_name"    : "h_wo",
            "s_current" : "Gauss1",
            "exp"       : "fwd",
            "mode"      : "JMEH",
            "name_JM"   : "JM_h",
            "name_EH"   : "EH_h",
            "device"    : device
            }

s.runPO(runPODict1)

s.autoConverge("EH_h", "e_wo", tol=1e-3)

runPODict2 = {
            "t_name"    : "e_wo",
            "s_current" : "JM_h",
            "exp"       : "fwd",
            "mode"      : "JMEH",
            "name_JM"   : "JM_e",
            "name_EH"   : "EH_e",
            "device"    : device
            }

s.runPO(runPODict2)

runPODict3 = {
            "t_name"    : "focus2",
            "s_current" : "JM_e",
            "exp"       : "fwd",
            "mode"      : "EH",
            "name_EH"   : "EH_f2",
            "device"    : device
            }

s.runPO(runPODict3)

2026-03-20 15:10:03 - WORK - *** Starting PO propagation *** 
2026-03-20 15:10:03 - WORK - Propagating Gauss1 on focus1 to h_wo, propagation mode: JMEH. 
2026-03-20 15:10:03 - WORK - Hardware: running 256 CUDA threads per block. 
2026-03-20 15:10:03 - WORK - ... Calculating ... 
2026-03-20 15:10:07 - WORK - *** Finished: 3.877 seconds *** 
2026-03-20 15:10:07 - WORK - *** Starting auto-convergence ***  
2026-03-20 15:10:08 - WORK - Difference : 1e+99 at gridsize ('1.100e+01', '1.100e+01') 
2026-03-20 15:10:09 - WORK - Difference : 0.0405691343873108 at gridsize ('2.100e+01', '2.100e+01') 
2026-03-20 15:10:11 - WORK - Difference : 0.008350533620479272 at gridsize ('3.100e+01', '3.100e+01') 
2026-03-20 15:10:12 - WORK - Difference : 0.008464968578108356 at gridsize ('4.100e+01', '4.100e+01') 
2026-03-20 15:10:13 - WORK - Difference : 0.0027531237822962518 at gridsize ('5.100e+01', '5.100e+01') 
2026-03-20 15:10:14 - WORK - Difference : 0.0008179420554712546 at gridsize ('6.100e+01', '6.100e+01') 
2026-03-20 15:10:14 - RESULT - Found converged solution, gridsize: ('5.490e+02', '5.490e+02') 
2026-03-20 15:10:14 - WORK - *** Starting PO propagation *** 
2026-03-20 15:10:14 - WORK - Propagating JM_h on h_wo to e_wo, propagation mode: JMEH. 
2026-03-20 15:10:14 - WORK - Hardware: running 256 CUDA threads per block. 
2026-03-20 15:10:14 - WORK - ... Calculating ... 
2026-03-20 15:11:53 - WORK - *** Finished: 98.707 seconds *** 
2026-03-20 15:11:53 - WORK - *** Starting PO propagation *** 
2026-03-20 15:11:53 - WORK - Propagating JM_e on e_wo to focus2, propagation mode: EH. 
2026-03-20 15:11:53 - WORK - Hardware: running 256 CUDA threads per block. 
2026-03-20 15:11:53 - WORK - ... Calculating ... 
2026-03-20 15:11:54 - WORK - *** Finished: 1.336 seconds ***
Out[2]:
<PyPO.PyPOTypes.fields at 0x7f12a4e38e10>
In [3]:
s.plotBeam2D("Gauss1", FieldComponents.Ez, project=Projections.yz, vmin=-30)
s.plotBeam2D("EH_f2", FieldComponents.Ex, project=Projections.xy, vmin=-30)

E1, H1 = s.calcHPBW("Gauss1", FieldComponents.Ex)
E2, H2 = s.calcHPBW("EH_f2", FieldComponents.Ex)

mean_mag = (E2/E1 + H2/H1) / 2

print(f"Focus 1: E-HPBW = {E1:.2f}, H-HPBW = {H1:.2f} degrees")
print(f"Focus 2: E-HPBW = {E2:.2f}, H-HPBW = {H2:.2f} degrees")
print(f"Mean PO magnification: M = {mean_mag:.2f}")
print(f"Fractional error: dM/M = {(mean_mag - 1.89)/1.89:.2f}")
No description has been provided for this image
No description has been provided for this image 
Focus 1: E-HPBW = 5.88, H-HPBW = 5.88 degrees
Focus 2: E-HPBW = 11.23, H-HPBW = 11.23 degrees
Mean PO magnification: M = 1.91
Fractional error: dM/M = 0.01

The top plot shows the Gaussian beam in focus 1, i.e. the starting beam pattern. Below, we find the beam pattern in focus 2, i.e. the termination point of the simulation. It can be seen, qualitatively, that both beam patterns look similar.

Quantitatively, the HPBWs indicate that the Dragonian magnifies by a factor of 1.91. This is not the same as the design, but relatively close. In fact, the fractional error, simulated minus design magnification over design, is about 1%.