Merge pull request #76 from symbiotic-engineering/radial-eigenfunction-fix

R_1n radial eigenfunction fix

Author: Y-Bimali

.gitignore

@@ -206,4 +206,12 @@ html/
 *.mat
 
 # vscode configuration files
-.vscode/
+.vscode/
+
+.calkit/overleaf
+
+# Calkit build artifacts
+.calkit/env-locks/
+.calkit/logs/
+.calkit/notebooks/executed/
+.calkit/notebooks/cleaned/

package/src/openflash/multi_equations.py

@@ -195,18 +195,11 @@ def R_1n_vectorized(n, r, i, h, d, a):
     # ---------------------------------------------------------------------------------
 
     cond_n_is_zero = (n == 0)
-    cond_r_at_boundary = (r == scale(a)[i])
+    cond_r_at_scale = (r == scale(a)[i])
 
     # --- Define the outcomes for each condition ---
-    # FIX: n=0 logic must match scalar version
-    if i == 0:
-        outcome_for_n_zero = np.full_like(r, 0.5)
-    else:
-        # Annulus: 1.0 + 0.5 * log(r / outer)
-        # Use np.divide to handle r=0 gracefully just in case
-        with np.errstate(divide='ignore', invalid='ignore'):
-            val = 1.0 + 0.5 * np.log(np.divide(r, scale(a)[i]))
-        outcome_for_n_zero = val
+    # Outcome 1: n=0 eigenfunction (constant)
+    outcome_for_n_zero = np.full_like(r, 0.5)
 
     # Outcome 2: If n>=1 and r is at the boundary, the value is 1.0.
     outcome_for_r_boundary = 1.0
@@ -221,7 +214,7 @@ def R_1n_vectorized(n, r, i, h, d, a):
                   exp(safe_lambda * (r - scale(a)[i]))
 
     # --- Apply the logic using nested np.where ---
-    result_if_n_not_zero = np.where(cond_r_at_boundary, outcome_for_r_boundary, bessel_term)
+    result_if_n_not_zero = np.where(cond_r_at_scale, outcome_for_r_boundary, bessel_term)
 
     return np.where(cond_n_is_zero, outcome_for_n_zero, result_if_n_not_zero)
 
@@ -237,11 +230,8 @@ def diff_R_1n_vectorized(n, r, i, h, d, a):
 
     condition = (n == 0)
 
-    if i == 0:
-        value_if_true = np.zeros_like(r)
-    else:
-        # Derivative is 1/(2r)
-        value_if_true = np.divide(1.0, 2 * r, out=np.full_like(r, np.inf), where=(r!=0))
+    # n=0 eigenfunction is constant -> derivative is zero
+    value_if_true = np.zeros_like(r)
 
     # --- Calculation for when n > 0 ---
     lambda_val = lambda_ni(n, i, h, d)
@@ -276,31 +266,31 @@ def R_2n_vectorized(n, r, i, a, h, d):
     # -------------------------------------------
 
     # LEGACY: Use Outer Radius
-    outer_r = scale(a)[i]
+    scale_r = scale(a)[i]
 
     cond_n_is_zero = (n == 0)
-    cond_r_at_boundary = (r == outer_r)
+    cond_r_at_scale = (r == scale_r)
 
     # Case 1: n = 0
     with np.errstate(divide='ignore', invalid='ignore'):
-         outcome_for_n_zero = 0.5 * np.log(np.divide(r, outer_r))
+         outcome_for_n_zero = 0.5 * np.log(np.divide(r, scale_r))
 
     # Case 2: n > 0 and r is at the boundary
-    outcome_for_r_boundary = 1.0
+    outcome_for_r_scale = 1.0
 
     # Case 3: n > 0 and r is not at the boundary
     lambda_val = lambda_ni(n, i, h, d)
 
     # Mask input where n=0 to prevent 'inf' errors
     lambda_safe = np.where(cond_n_is_zero, 1.0, lambda_val)
 
-    # LEGACY: Denom uses OUTER radius
-    denom = besselke(0, lambda_safe * outer_r)
+    # LEGACY: Denom uses scaled radius
+    denom = besselke(0, lambda_safe * scale_r)
     denom = np.where(np.abs(denom) < 1e-12, np.nan, denom)
 
-    bessel_term = (besselke(0, lambda_safe * r) / denom) * exp(lambda_safe * (outer_r - r))
+    bessel_term = (besselke(0, lambda_safe * r) / denom) * exp(lambda_safe * (scale_r - r))
 
-    result_if_n_not_zero = np.where(cond_r_at_boundary, outcome_for_r_boundary, bessel_term)
+    result_if_n_not_zero = np.where(cond_r_at_scale, outcome_for_r_scale, bessel_term)
 
     return np.where(cond_n_is_zero, outcome_for_n_zero, result_if_n_not_zero)
 
@@ -316,19 +306,19 @@ def diff_R_2n_vectorized(n, r, i, h, d, a):
 
     # Case n > 0
     lambda_val = lambda_ni(n, i, h, d)
-    outer_r = scale(a)[i] # LEGACY ANCHOR
+    scale_r = scale(a)[i] # LEGACY ANCHOR
 
     lambda_safe = np.where(n == 0, 1.0, lambda_val)
 
-    # LEGACY: Denom uses OUTER radius
-    denom = besselke(0, lambda_safe * outer_r)
+    # LEGACY: Denom uses scale radius
+    denom = besselke(0, lambda_safe * scale_r)
     safe_denom = np.where(np.abs(denom) < 1e-10, 1e-10, denom)
 
     with np.errstate(divide='ignore', invalid='ignore'):
         numerator = -lambda_safe * besselke(1, lambda_safe * r)
         ratio = numerator / safe_denom
         # LEGACY: Exponential decay from OUTER radius
-        exp_term = exp(lambda_safe * (outer_r - r))
+        exp_term = exp(lambda_safe * (scale_r - r))
         value_if_false = ratio * exp_term
 
     return np.where(n == 0, value_if_true, value_if_false)
@@ -386,9 +376,9 @@ def Lambda_k_vectorized(k, r, m0, a, m_k_arr):
     # -------------------------------------------
 
     cond_k_is_zero = (k == 0)
-    cond_r_at_boundary = (r == scale(a)[-1])
+    cond_r_at_scale = (r == scale(a)[-1])
 
-    outcome_boundary = 1.0
+    outcome_scale = 1.0
 
     # --- Case 2: k = 0 ---
     if m0 == inf:
@@ -426,13 +416,13 @@ def Lambda_k_vectorized(k, r, m0, a, m_k_arr):
                                  where=np.isfinite(denom_k_nonzero) & (denom_k_nonzero != 0))
     outcome_k_nonzero = bessel_ratio * exp(safe_m_k * (scale(a)[-1] - r))
 
-    result_if_not_boundary = np.where(cond_k_is_zero,
+    result_if_not_scale = np.where(cond_k_is_zero,
                                       outcome_k_zero,
                                       outcome_k_nonzero)
 
-    return np.where(cond_r_at_boundary,
-                    outcome_boundary,
-                    result_if_not_boundary)
+    return np.where(cond_r_at_scale,
+                    outcome_scale,
+                    result_if_not_scale)
 
 # Differentiate wrt r 
 def diff_Lambda_k_vectorized(k, r, m0, a, m_k_arr):
@@ -567,27 +557,11 @@ def diff_Z_k_e_vectorized(k, z, m0, h, m_k_arr, N_k_arr):
 #integrating R_1n * r
 # Integration
 def int_R_1n(i, n, a, h, d):
-    if n == 0:
+    if n == 0: # Integral of 0.5 * r
         if i == 0:
-            # Central cylinder: Integral of 0.5 * r
             return a[i]**2/4 
         else:
-            # Annulus: Integral of r * (1.0 + 0.5 * ln(r/outer_r))
-            outer_r = scale(a)[i]
-            inner_r = a[i-1]
-            
-            cyl_term = (outer_r**2 - inner_r**2) / 2.0
-            
-            def log_indefinite_int(r):
-                # Integral of 0.5 * r * ln(r/outer)
-                # = 0.5 * [ (r^2/2)*ln(r/outer) - r^2/4 ]
-                log_val = np.log(r/outer_r) if r > 0 else 0
-                return 0.5 * ((r**2 / 2.0) * log_val - (r**2 / 4.0))
-
-            val_outer = log_indefinite_int(outer_r) # log(1)=0, so just -r^2/4 term
-            val_inner = log_indefinite_int(inner_r)
-            
-            return cyl_term + (val_outer - val_inner)
+            return (a[i]**2 - a[i-1]**2)/4
     else:
         # Standard Bessel I integral (Unchanged)
         lambda0 = lambda_ni(n, i, h, d)
@@ -608,8 +582,9 @@ def int_R_2n(i, n, a, h, d):
         raise ValueError("i cannot be 0")
 
     # LEGACY: Use Outer Radius
-    outer_r = scale(a)[i]
+    outer_r =a[i]
     inner_r = a[i-1] # Previous radius is inner
+    scale_r = scale(a)[i]
 
     if n == 0:
         # Integral of r * (0.5 * ln(r/outer_r))
@@ -618,7 +593,7 @@ def int_R_2n(i, n, a, h, d):
         def indefinite(r):
             if r <= 0: return 0
             # ln(r/outer_r) is 0 when r=outer_r
-            term_log = np.log(r / outer_r)
+            term_log = np.log(r / scale_r)
             return 0.5 * ((r**2 / 2) * term_log - (r**2 / 4))
 
         val_outer = indefinite(outer_r)
@@ -629,13 +604,14 @@ def indefinite(r):
         lambda0 = lambda_ni(n, i, h, d)
 
         term_outer = (outer_r * besselke(1, lambda0 * outer_r)) 
-        # exp(la - la) = 1, so no exp factor needed for outer term.
+        # exp(la - la) = 1, for scale = outer_r, so no exp factor needed for outer term in that case.
+        term_outer *= np.exp(lambda0 * (scale_r - outer_r))
 
         term_inner = (inner_r * besselke(1, lambda0 * inner_r))
         # exp(la - li). We need to multiply by this.
-        term_inner *= np.exp(lambda0 * (outer_r - inner_r))
+        term_inner *= np.exp(lambda0 * (scale_r - inner_r))
 
-        denom = lambda0 * besselke(0, lambda0 * outer_r)
+        denom = lambda0 * besselke(0, lambda0 * scale_r)
 
         return (term_inner - term_outer) / denom
 #integrating phi_p_i * d_phi_p_i/dz * r *d_r at z=d[i]
@@ -657,8 +633,8 @@ def z_n_d(n):
 #############################################
 def excitation_phase(x, NMK, m0, a): # x-vector of unknown coefficients
     coeff = x[-NMK[-1]] # first coefficient of e-region expansion
-    local_scale = scale(a)
-    return -(pi/2) + np.angle(coeff) - np.angle(besselh(0, m0 * local_scale[-1]))
+    exterior_scale = scale(a)[-1]
+    return -(pi/2) + np.angle(coeff) - np.angle(besselh(0, m0 * exterior_scale))
 
 def excitation_force(damping, m0, h):
     # --- FIX: Handle infinite m0 ---

package/test/test_multi_equations.py

@@ -301,23 +301,12 @@ def test_int_R_1n_n0_i0(a, h, d):
     assert np.isclose(int_R_1n(0, 0, a, h, d), expected)
 
 def test_int_R_1n_n0_i_positive(a, h, d):
-    # Fixed expected value: Integral of r * (1.0 + 0.5 * log(r/outer))
+    # Fixed expected value: Integral of r * 0.5
     i = 1
     outer_r = a[i]
     inner_r = a[i-1]
 
-    # 1. Cylinder term (integral of r*1.0)
-    cyl_term = (outer_r**2 - inner_r**2) / 2.0
-    
-    # 2. Log term helper
-    def log_indefinite_int(r):
-        log_val = np.log(r/outer_r) if r > 0 else 0
-        return 0.5 * ((r**2 / 2.0) * log_val - (r**2 / 4.0))
-
-    val_outer = log_indefinite_int(outer_r) 
-    val_inner = log_indefinite_int(inner_r)
-    
-    expected = cyl_term + (val_outer - val_inner)
+    expected = (outer_r**2 - inner_r**2) / 4
     assert np.isclose(int_R_1n(1, 0, a, h, d), expected)
 
 def test_int_R_1n_n_positive(test_n, test_i, a, h, d):

sea-lab-utils

@@ -0,0 +1 @@
+Subproject commit f9d69427a24304292c4a140b272e63122df9b6cf