From f32ba924d419fb4037164cbe1b8e1e6c4dbf3ba6 Mon Sep 17 00:00:00 2001 From: Joseph Capriotti Date: Tue, 2 Dec 2025 10:42:12 -0700 Subject: [PATCH 1/3] Allows complex mu values in rTE kernel Updates the rTE kernel to accept complex values for mu, improving its applicability to a wider range of electromagnetic modeling scenarios. This change ensures that the mu parameter, representing permeability, can be complex, which is necessary for handling materials with magnetic losses or other complex electromagnetic properties. --- geoana/kernels/_extensions/_rTE.cpp | 15 ++++++++------- geoana/kernels/_extensions/_rTE.h | 8 ++++---- geoana/kernels/_extensions/rTE.pyx | 16 ++++++++-------- 3 files changed, 20 insertions(+), 19 deletions(-) diff --git a/geoana/kernels/_extensions/_rTE.cpp b/geoana/kernels/_extensions/_rTE.cpp index e0d905d9..e3d77f20 100644 --- a/geoana/kernels/_extensions/_rTE.cpp +++ b/geoana/kernels/_extensions/_rTE.cpp @@ -9,19 +9,20 @@ void funcs::rTE( double * frequencies, double * lambdas, complex_t * sigmas, - double * mus, + complex_t * mus, double * h, size_t n_frequency, size_t n_filter, size_t n_layers -){ +) noexcept(true){ double mu_0 = 4.0E-7*M_PI; complex_t j(0.0, 1.0); complex_t k2, u, Yh, Y, tanhuh, Y0; complex_t *sigi; - double *mui; - double omega, mu_c, l2; + complex_t *mui; + double omega, l2; + complex_t mu_c; for(size_t i_filt=0, i=0; i_filt Date: Tue, 2 Dec 2025 10:46:31 -0700 Subject: [PATCH 2/3] Explicitly cast mu to complex_t pointer --- geoana/kernels/_extensions/rTE.pyx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/geoana/kernels/_extensions/rTE.pyx b/geoana/kernels/_extensions/rTE.pyx index edbcf182..5024acac 100644 --- a/geoana/kernels/_extensions/rTE.pyx +++ b/geoana/kernels/_extensions/rTE.pyx @@ -109,12 +109,13 @@ def rTE_forward(frequencies, lamb, sigma, mu, thicknesses): # Types are same size and both pairs of numbers (r, i), so can cast one to the other complex_t *out_p = &out[0, 0] complex_t *sig_p = &sig[0, 0] + complex_t *mu_p = &c_mu[0, 0] REAL_t *h_p = NULL if n_layers > 1: h_p = &hs[0] with nogil: - rTE(out_p, &f[0], &lam[0], sig_p, &c_mu[0, 0], h_p, + rTE(out_p, &f[0], &lam[0], sig_p, mu_p, h_p, n_frequency, n_filter, n_layers) return np.array(out) @@ -194,6 +195,7 @@ def rTE_gradient(frequencies, lamb, sigma, mu, thicknesses): cdef: # Types are same size and both pairs of numbers (r, i), so can cast one to the other complex_t *sig_p = &sig[0, 0] + complex_t *mu_p = &c_mu[0, 0] complex_t *gsig_p = &gsig[0, 0, 0] complex_t *gmu_p = &gmu[0, 0, 0] complex_t *gh_p = NULL @@ -203,7 +205,7 @@ def rTE_gradient(frequencies, lamb, sigma, mu, thicknesses): gh_p = &gh[0, 0, 0] with nogil: - rTEgrad(gsig_p, gmu_p, gh_p, &f[0], &lam[0], sig_p, &c_mu[0, 0], h_p, + rTEgrad(gsig_p, gmu_p, gh_p, &f[0], &lam[0], sig_p, mu_p, h_p, n_frequency, n_filter, n_layers) return np.array(gsig), np.array(gh), np.array(gmu) From 1bea0df4901e3802477b2930d4e4b34c07f2a96b Mon Sep 17 00:00:00 2001 From: Joseph Capriotti Date: Tue, 2 Dec 2025 12:34:57 -0700 Subject: [PATCH 3/3] test kernels with complex conductivity and complex permeability --- tests/em/fdem/test_em_fdem_layered.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/em/fdem/test_em_fdem_layered.py b/tests/em/fdem/test_em_fdem_layered.py index 1f6bcf77..aced0961 100644 --- a/tests/em/fdem/test_em_fdem_layered.py +++ b/tests/em/fdem/test_em_fdem_layered.py @@ -294,9 +294,9 @@ def test_rTE_jacobian(self): thicknesses = np.ones(n_layer-1) lamb = np.logspace(0, 3, n_lambda) np.random.seed(0) - sigma = 1E-1*(1 + 0.1 * np.random.rand(n_layer, n_frequency)) - mu = mu_0 * (1 + 0.1 * np.random.rand(n_layer, n_frequency)) - dmu = mu_0 * 0.1 * np.random.rand(n_layer, n_frequency) + sigma = 1E-1*(1 + 1E-4j + 0.1 * np.random.rand(n_layer, n_frequency)) + mu = mu_0 * (1 + 1.0E-4j + 0.1 * np.random.rand(n_layer, n_frequency)) + dmu = mu_0 * (0.1 + 1.0E-5j) * np.random.rand(n_layer, n_frequency) def rte_sigma(x): sigma = x.reshape(n_layer, n_frequency) @@ -354,8 +354,8 @@ def setUp(self): thicknesses = np.ones(n_layer-1) lamb = np.logspace(-5, -3, n_lambda) np.random.seed(123) - sigma = 1E-1 * (1 + 1.0/(n_layer*n_frequency) * np.arange(n_layer*n_frequency).reshape(n_layer, n_frequency)) - mu = mu_0 * (1 + 1.0/(n_layer*n_frequency) * np.arange(n_layer*n_frequency).reshape(n_layer, n_frequency)) + sigma = 1E-1 * (1 + 1E-4j + 1.0/(n_layer*n_frequency) * np.arange(n_layer*n_frequency).reshape(n_layer, n_frequency)) + mu = mu_0 * (1 + 1.0E-4j + 1.0/(n_layer*n_frequency) * np.arange(n_layer*n_frequency).reshape(n_layer, n_frequency)) self.rTE1 = rTE_forward(frequencies, lamb, sigma, mu, thicknesses) self.rTE2 = _rTE_forward(frequencies, lamb, sigma, mu, thicknesses)