AddWaveGuidePort for circular WG port, Python/Windows

[g4dbn]

Before I go too deep into this Python 3.11 version 0.0.36 in Windows 11 under VSCode, I'd like to known whether that build version will support creating circular waveguide ports using AddWaveGuidePort() with suitable E and H expressions? I don't see a round guide class, and when I try to create a round guide with TE11, I get hundred of errors like this:

CSPropExcitation::GetWeightedExcitation: Error evaluating the weighting function (ID: 2, n=1):
CSPropExcitation::GetWeightedExcitation: Error evaluating the weighting function (ID: 2, n=0):

I suspect I'm doing something wrong with the expressions. Are there any examples that show the correct way to express a simple TE11 in terms of the E and H expressions? I've tried this, but without some solid examples of what they are supposed to look like - and confirmation that the function will work with v0.0.36, python 3.11 and Windows, I'm a bit stuck and tempted to ditch Python and go back to Octave, or maybe just ditch Windows.

E_rho = 'sin(alpha)'
E_phi = 'cos(alpha)'
E_z   = '0'

H_rho = 'cos(alpha)'
H_phi = '-1 * sin(alpha)'
H_z   = '0'

E_func = [E_rho, E_phi, E_z]
H_func = [H_rho, H_phi, H_z]

ports = []
ports.append( FDTD.AddWaveGuidePort(1, start_port1, stop_port1, 'z', E_func, H_func, kc_val, excite=1) )

[thliebig]

Where does "alpha" come from? I think in openEMS the angle around the z-axis is just "a".
These are all the available "coordinate" names: "x,y,z,rho,a,r,t"

[g4dbn]

That was from an attempt to debug what was going on after I misunderstood a note about coordinate names that I found on a random forum sorry!

I've tried to simplify the example to tease out what is going on. It feels like the E_func and H_func values being passed to AddWaveGuidePort() are not being parsed correctly because I've got something wrong.

I thought it might be divide by zero for rho, so added an epsilon offset variable, I thought it might be the excitation overlapping the metal of the waveguide, so I made the square excitation smaller so it would fit inside the guide, I tried formatting the the kc variable as a float, as an exponent, as a literal, I tried various formats for the other values, I removed all the numbers and set zero, tried zero for kc, tried literals, removed the divisor /, put it back, tried shrinking the excitation to minimal size, put it back, changed the grid settings, tried setting different lengths of the excitation from 0 upwards to 4 cells, and zero failed of course, but nothing worked.

I'd love to see a minimal example of a round waveguide simulation in the Python/Windows environment so I can try to understand what on earth I'm doing wrong. Is it something related to my use of CSX = ContinuousStructure() perhaps? I've made so many edits, I can't remember why I ended up setting zero for some of the H_func elements. Just seeing a working example would be a huge help, or some confirmation that it should work if I get the H and E functions correct.

When I call appCSXCAD, it shows the structure and ports correctly. I've left the small-sized square excitations in place, although I've never had a problem with the excitation port overlapping metal in Octave, so I don't know whether that matters.

My reason for wanting to migrate to the Python stuff is to try out some genetic algorithms for optimising feedhorn sidelobe levels, and that seemed to be easier with Python than Octave. My reason for using Windows is that my video editing machine is the most powerful processor I have, and I want to integrate the output field animations via Blender into Davinci Resolve for Youtube, and use Autodesk Fusion to create the solids for import, without having to buy another hugely powerful computer. Am I setting myself up to fail by sticking with Windows?
Current terminal output with the error and versions is here:

C:/Python311/python.exe c:/opt/openEMS/python/Tutorials/circ_waveguide_test.py
E_func: ['(j1(0.024880861892*rho)/(0.024880861892*rho+1.0e-12))*cos(a)', '(j1(0.024880861892*rho)/(0.024880861892*rho+1.0e-12))*sin(a)', '0']
H_func: ['0', '0', '(j1(0.024880861892*rho)/(0.024880861892*rho+1.0e-12))*sin(a)']
kc_val 0.02488086189189189
 ---------------------------------------------------------------------- 
 | openEMS 64bit -- version v0.0.36
 | (C) 2010-2023 Thorsten Liebig <thorsten.liebig@gmx.de>  GPL license
 ---------------------------------------------------------------------- 
        Used external libraries:
                CSXCAD -- Version: v0.6.3
                hdf5   -- Version: 1.12.0
                          compiled against: HDF5 library version: 1.12.0
                tinyxml -- compiled against: 2.6.2
                fparser
                boost  -- compiled against: 1_72
                vtk -- Version: 8.2.0
                       compiled against: 8.2.0

Create FDTD operator (compressed SSE + multi-threading)
Multithreaded operator using 8 threads.
Operator::CalcTimestep_Var3: Smallest timestep (2.22376e-09s) found at position: 0 : 0;0;0
Compressing the FDTD operator... this may take a while...
Excitation::buildExcitationSignal: Unknown error... excitation setup failed!!

My much-amended code without any postprocessing is here:

import os
import numpy as np
from openEMS import openEMS
from openEMS.physical_constants import *
from CSXCAD import ContinuousStructure
import tempfile

# --- Simulation parameters ---
SimBox = [200, 200, 400]
radius = 75
thickness = 2

# --- FDTD setup ---
FDTD = openEMS(EndCriteria=1e-5)
f_start = 1e9
f_stop  = 4e9
FDTD.SetGaussExcite((f_start+f_stop)/2, (f_stop-f_start)/2)
FDTD.SetBoundaryCond(['PML_8']*6)

# --- Geometry and mesh ---
CSX = ContinuousStructure()
FDTD.SetCSX(CSX)
mesh = CSX.GetGrid()
x_lines = np.linspace(-SimBox[0]/2, SimBox[0]/2, int(SimBox[0]/1)+1)
y_lines = np.linspace(-SimBox[1]/2, SimBox[1]/2, int(SimBox[1]/1)+1)
z_lines = np.linspace(0, SimBox[2], int(SimBox[2]/2)+1)
CSX.DefineGrid({'x':x_lines, 'y':y_lines, 'z':z_lines}, 1)

metal = CSX.AddMetal('metal')
wg_air = CSX.AddMaterial('wg_air')
start_geom = [0, 0, 0]
stop_geom  = [0, 0, SimBox[2]]
metal.AddCylinder(start_geom, stop_geom, radius)
wg_air.AddCylinder(start_geom, stop_geom, radius-thickness, priority=10)

# --- Port setup (TE11 mode) ---
p11_p = 1.84118378
kc_val = p11_p / (radius-thickness/2)
kc = kc_val
epsilon = 1e-12

radius_air = radius - thickness
margin = 2
port_xy = (radius_air - margin) / np.sqrt(2)
start_port1 = [-port_xy, -port_xy, 20]
stop_port1  = [ port_xy,  port_xy, 22]
start_port2 = [-port_xy, -port_xy, 378]
stop_port2  = [ port_xy,  port_xy, 380]

# TE11 mode field expressions (openEMS expects 'j1', 'rho', 'a')
E_func = [
    "(j1({:.12f}*rho)/({:.12f}*rho+{:.1e}))*cos(a)".format(kc, kc, epsilon),
    "(j1({:.12f}*rho)/({:.12f}*rho+{:.1e}))*sin(a)".format(kc, kc, epsilon),
    "0"
]
H_func = [
    "0",
    "0",
    "(j1({:.12f}*rho)/({:.12f}*rho+{:.1e}))*sin(a)".format(kc, kc, epsilon)
]
print("E_func:", E_func)
print("H_func:", H_func)
print("kc_val", kc_val)

# --- Add waveguide ports ---
ports = []
ports.append(FDTD.AddWaveGuidePort(1, start_port1, stop_port1, 'z', E_func, H_func, kc_val, excite=1))
ports.append(FDTD.AddWaveGuidePort(2, start_port2, stop_port2, 'z', E_func, H_func, kc_val))

# --- Write CSX for debugging ---
Sim_Path = os.path.join(tempfile.gettempdir(), 'Circ_WG_SM')
CSX_file = os.path.join(Sim_Path, 'circ_wg.xml')
if not os.path.exists(Sim_Path):
    os.mkdir(Sim_Path)
CSX.Write2XML(CSX_file)

# --- Run simulation ---
FDTD.Run(Sim_Path, verbose=2)

[g4dbn]

Dave points out that the error is generated if GetNyquistNum() in excitation.cpp returns 0, that is the error code I am seeing. He has tested on a Linux machine with Python and it fails as well , so is there something wrong with my script connected with the SetGaussExcite() perhaps?

My code has

f_start = 1e9
f_stop  = 4e9
FDTD.SetGaussExcite((f_start+f_stop)/2, (f_stop-f_start)/2)

and in excitation.cpp there is

void Excitation::SetupGaussianPulse(double f0, double fc)
{
	m_Excit_Type = Excitation::GaissianPulse;
	m_f0 = f0;
	m_fc = fc;
	m_f_max = f0+fc;
	m_SignalPeriod = 0;
}

and in
Excitation::buildExcitationSignal(unsigned int maxTS)
there is

switch (m_Excit_Type)
	{
	case Excitation::GaissianPulse:
		CalcGaussianPulsExcitation(m_f0,m_fc,maxTS);
		break;
...
        if (GetNyquistNum() == 0)
	{
		cerr << "Excitation::buildExcitationSignal: Unknown error... excitation setup failed!!" << endl;
		return false;
	}

Is there anything I need to change to avoid that error?