Skip to content
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
243 changes: 238 additions & 5 deletions docs/phase2-matrix-fill/math.md
Original file line number Diff line number Diff line change
Expand Up @@ -191,7 +191,7 @@ For the self-impedance element Z[m,m], the source and observation points are on

With the exact kernel, R_surf ≥ a > 0 — the singularity does not occur. However, when a is small, R_surf is small for s ≈ s' and the integrand is sharply peaked. Standard fixed-order Gauss-Legendre quadrature undersamples the peak and produces inaccurate results. This is the near-singular case.

### 5.2 Treatment — Singularity Extraction
### 5.2 Treatment — Near-Singular Quadrature (Exact Kernel)

The standard approach is singularity extraction. Split the self-impedance integral into a singular (or near-singular) part with a known analytic form, and a smooth remainder that is well-suited to Gauss-Legendre:

Expand All @@ -203,7 +203,97 @@ K_exact(s, s') = [K_exact(s, s') - K_singular(s, s')] + K_singular(s, s')

The analytic part K_singular is chosen to match the singular behavior of K_exact as s → s'. The smooth remainder is then integrated with standard Gauss-Legendre at modest quadrature order. The analytic part is integrated exactly.

The specific form of K_singular and its analytic integral are derived in terms of the segment geometry and wire radius. See the Fikioris references for the exact kernel self-impedance extraction.
**Conceptual clarification from Fikioris & Wu (2001).** That paper demonstrates that the
*approximate* (thin-wire) kernel K_ap diverges logarithmically at z = 0 — which is why the
approximate integral equation has no solution. The *exact* kernel K_ex is bounded everywhere.
For Arcanum's exact kernel, Z[m,m] has **no true singularity**. The challenge is
**near-singularity**: for thin wires (a << Δ), the integrand is sharply peaked near s = s'
because R_surf achieves its minimum value of a there. Standard fixed-order Gauss-Legendre
undersamples the peak. The correct treatment is smooth near-singular subtraction, not
classical singularity extraction.

**Straight-segment simplification.** For a straight segment, the observation point is on the
axis and the source surface point is at radius a. The distance R_surf is φ-independent:

```
R_surf = √(a² + (s-s')²) [φ drops out for straight wire]

K_exact(s, s') = e^{-jk√(a²+(s-s')²)} / √(a²+(s-s')²)
```

At s = s' this equals e^{-jka}/a — large for thin wires, but finite.

**Near-singular subtraction.** Define the near-singular part:

```
K_near(s, s') = 1 / √(a² + (s-s')²)
```

This matches the dominant behavior of K_exact near s = s'. The smooth remainder:

```
K_smooth(s, s') = K_exact(s, s') - K_near(s, s')
= [e^{-jk√(a²+(s-s')²)} - 1] / √(a²+(s-s')²)
```

As s → s': K_smooth → (e^{-jka} - 1)/a → -jk. Bounded and smooth everywhere.
Standard Gauss-Legendre integrates K_smooth accurately at order 8–16.

**Analytic double integral of K_near — T1 self-element.** Setting u = s - s_n, v = s' - s_n,
both in [0, Δ]:

```
I_near = ∫_0^Δ ∫_0^Δ 1/√(a²+(u-v)²) dv du
```

Inner integral: `∫_0^Δ 1/√(a²+(u-v)²) dv = sinh⁻¹((Δ-u)/a) + sinh⁻¹(u/a)`

Outer integral by symmetry (substituting u → Δ-u in the first term):

```
I_near = 2 ∫_0^Δ sinh⁻¹(u/a) du
= 2 [u sinh⁻¹(u/a) - √(a²+u²)]_0^Δ
= 2 [Δ sinh⁻¹(Δ/a) - √(a²+Δ²) + a]
```

This is the closed-form contribution from K_near to T1[m,m].

**Analytic single integral of K_near — T2 self-element.** The T2 term evaluates K_exact at
segment endpoints. For the self-element the observation point s sweeps the segment while the
source point is fixed at an endpoint s_end:

```
∫_{Δm} K_near(s, s_end) ds = sinh⁻¹((s_{m+1} - s_end)/a) - sinh⁻¹((s_m - s_end)/a)
```

T2[m,m] is bounded everywhere (R_surf ≥ a at any endpoint) and does not require the
near-singular subtraction. Standard Gauss-Legendre applies.

**Complete self-element evaluation:**

```
Z[m,m] = jωμ₀/(4π) × { I_near + ∫∫ K_smooth ds ds' }
- 1/(jωε₀4π) × ∫ [K_exact(s, r_m(s_{m+1})) - K_exact(s, r_m(s_m))] ds
```

where I_near is the closed form above and ∫∫ K_smooth is evaluated with adaptive
Gauss-Legendre.

**Thin-wire limit consistency check.** As a → 0, using sinh⁻¹(Δ/a) → ln(2Δ/a):

```
I_near → 2Δ [ln(2Δ/a) - 1] + O(a)
```

Multiplying by jωμ₀/(4π): `jωμ₀Δ/(2π) × [ln(2Δ/a) - 1]` — which is exactly the
thin-wire self-impedance formula in Section 5.3. ✓

**Curved-wire generalization.** For arc and helix segments, R_surf depends on φ because the
wire curves away from the observation point. The near-singular subtraction still applies in
the local tangent frame: near s ≈ s' the curved wire is locally straight, so K_near has the
same form with the same analytic integral. The curvature correction appears only in K_smooth,
which is smooth and integrable with standard methods. The precise K_exact for curved segments
requires the perpendicular frame (n̂_r, n̂_φ); see Section 7.4 for the construction.

### 5.3 Thin-Wire Limit of Self-Impedance

Expand All @@ -225,7 +315,77 @@ The logarithmic dependence on a/Δ confirms that the imaginary part of self-impe

For adjacent segments m and m+1 (sharing one endpoint), the minimum axis-to-axis distance R_axis can become very small as the integration variables s and s' approach the shared endpoint. With the exact kernel, R_surf ≥ a, so there is no true singularity — but for small a, the integrand is again sharply peaked near the shared endpoint.

The same singularity extraction approach applies. Identify the near-singular behavior analytically, subtract it, integrate the remainder with standard quadrature, add the analytic contribution.
The near-singular subtraction of Section 5.2 applies here too. The difference is geometry:
the near-singular region is concentrated near the **shared endpoint** P = r_m(s_{m+1}),
not spread over the full segments.

**Local coordinates near the shared endpoint.** Define:

```
ε = s_{m+1} - s [distance from P within segment m, ε ∈ [0, Δm]]
ε' = s' - s_{m+1} [distance from P within segment m+1, ε' ∈ [0, Δ_{m+1}]]
```

For straight segments meeting at angle θ_bend, near the shared endpoint:

```
R_surf ≈ √(a² + ε² + ε'² + 2εε' cos θ_bend)
```

For collinear segments (θ_bend = 0): R_surf ≈ √(a² + (ε + ε')²) — same form as
the self-element with the combined arc length (ε + ε').

**Near-singular subtraction for T1[m, m+1].** Define:

```
K_near(ε, ε') = (t̂_m · t̂_{m+1}) / √(a² + ε² + ε'² + 2εε' cos θ_bend)
```

where t̂_m · t̂_{m+1} = cos θ_bend. The smooth remainder is bounded for all ε, ε'.

**Analytic double integral of K_near.** For collinear segments (θ_bend = 0, cos θ_bend = 1),
the near-singular integral over the one-sided domain evaluates to:

```
∫_0^{Δm} ∫_0^{Δ_{m+1}} 1/√(a²+(ε+ε')²) dε' dε

= Δm sinh⁻¹(Δ_{m+1}/a) + Δ_{m+1} sinh⁻¹(Δm/a)
- √(a²+(Δm+Δ_{m+1})²) + √(a²+Δm²) + √(a²+Δ_{m+1}²) - a
```

Derivation: inner integral gives sinh⁻¹((Δ_{m+1}+ε)/a) - sinh⁻¹(ε/a); outer integral
applies the same antiderivative formula used in Section 5.2.

For a bent junction (θ_bend ≠ 0), the analytic integral acquires a correction factor from
cos θ_bend. For θ_bend small (gentle bends, typical in discretized curved wires), the
collinear formula is a good approximation. For sharp bends (θ_bend ≥ π/4), the correction
is non-negligible and the integral should be evaluated adaptively.

**Near-singular contribution to T2[m, m+1].** The T2 term for the near-neighbor element is:

```
T2[m, m+1] = ∫_{Δm} [K_exact(r_m(s), r_{m+1}(s_{m+2})) - K_exact(r_m(s), r_{m+1}(s_{m+1}))] ds
```

The second term K_exact(r_m(s), r_{m+1}(s_{m+1})) peaks as s → s_{m+1} (observation
approaches the shared endpoint P from within segment m, source fixed at P on segment m+1's
surface). Near-singular subtraction:

```
K_near_T2(s) = 1 / √(a² + (s_{m+1}-s)²)

Analytic integral: ∫ K_near_T2 ds = sinh⁻¹((s_{m+1}-s)/a) + C
```

The first term K_exact(r_m(s), r_{m+1}(s_{m+2})) is bounded away from any singularity and
requires no special treatment.

**`near_singular_distance_ratio` interpretation.** The design doc `MatrixFillConfig` field
`near_singular_distance_ratio` (default 3.0) defines how far from the shared endpoint the
near-singular treatment is applied, measured in multiples of the wire radius a. Concretely:
the near-singular subtraction is applied for ε < 3a and ε' < 3a. Beyond this distance the
integrand is smooth and standard Gauss-Legendre at order 8 is accurate. A value of 3.0 is
conservative; for most wire geometries 2.0 suffices.

### 6.2 Extent of Near-Neighbor Treatment

Expand Down Expand Up @@ -259,6 +419,77 @@ The default convergence threshold is ε_quad = 1×10⁻¹⁰ (10 significant fig

The azimuthal integral in K_exact is evaluated with a fixed-order Gauss-Legendre rule on [0, 2π]. For most wire radii (a ≤ Δ/2), order 16 (16 azimuthal points) achieves 8 significant figures in K_exact. For a/Δ > 0.5 (very thick wires), higher azimuthal order may be required.

### 7.4 Perpendicular Frame Computation for Exact Kernel

**Frame continuity is not required — proof.** The exact kernel azimuthal integral is:

```
K_exact(s, s') = 1/(2π) ∫₀^{2π} G₀(r_obs, r_axis + a[n̂_r cos φ + n̂_φ sin φ]) dφ
```

If the frame is rotated by angle α — replacing n̂_r with n̂_r cosα + n̂_φ sinα — the
integrand becomes G₀(..., φ + α). Since the integration is over a full period [0, 2π],
the result is identical for any α. Therefore:

**Any orthonormal pair {n̂_r, n̂_φ} perpendicular to t̂ produces the same K_exact.**
Each segment may use its own independently computed frame. Inter-segment continuity is
not required and need not be enforced. ✓

**Single construction for all curve types.** Since frame orientation is irrelevant, one
numerically stable algorithm serves all segment types. Given t̂ at any source point s':

```
Step 1 — choose reference vector v:
v = argmin over {x̂, ŷ, ẑ} of |v · t̂|
(the standard basis vector least aligned with t̂)

Step 2 — Gram-Schmidt orthogonalization:
n̂_r = (v - (v·t̂) t̂) / |v - (v·t̂) t̂|

Step 3 — right-hand perpendicular:
n̂_φ = t̂ × n̂_r
```

The minimum-component choice in Step 1 guarantees |v - (v·t̂)t̂| ≥ √(2/3) > 0
for all unit vectors t̂ — no special-case fallback needed.

**Verification by segment type.**

*Straight segment along ẑ:* t̂ = (0, 0, 1). Min-component basis vector: x̂ (or ŷ).
Taking v = x̂: n̂_r = x̂, n̂_φ = ẑ × x̂ = ŷ. Standard cylindrical frame. ✓

*Arc segment in xy-plane:* t̂(θ) = (−sin θ, cos θ, 0). The z-component of t̂ is zero,
so v = ẑ is the minimum-component choice. Since ẑ · t̂ = 0:

```
n̂_r = ẑ
n̂_φ = t̂ × ẑ = (cos θ, sin θ, 0) [outward radial direction]
```

The surface point r_surf = r_axis + a[ẑ cos φ + r̂ sin φ] traces the correct circular
cross-section perpendicular to t̂ at every θ. ✓

*Helix segment:* t̂ ∝ (−R sin θ, R cos θ, p/2π) where R is the helix radius and p is
the pitch. For a typical helix (R > p/2π) the z-component of t̂ is the smallest, so
v = ẑ again. For a steep helix (p/2π > R) the radial component is smallest. The
algorithm selects the correct v automatically in either case; no closed-form expression
for the frame is needed — t̂(s') is evaluated numerically from the parametric forms in
`docs/phase1-geometry/math.md` Sections 4–5, and the three-step construction above
completes in O(1) per quadrature point. ✓

**Relationship to Frenet-Serret.** For arc and helix segments, the Frenet-Serret
principal normal n̂_FS = (dt̂/ds)/|dt̂/ds| is also perpendicular to t̂. It would produce
the same K_exact by the dummy-variable argument above. The minimum-component
construction is preferred here because it requires only t̂ (already available from
phase1-geometry), not the curvature κ = |dt̂/ds|, and it is well-defined for straight
segments where κ = 0 and the Frenet-Serret frame is undefined.

**Implementation summary.** The function `evaluate_exact_kernel` (design.md Section 6.4)
receives t̂ as an input. It computes {n̂_r, n̂_φ} via the three-step construction above,
then evaluates the azimuthal quadrature. No segment-type branching is required inside
`evaluate_exact_kernel`; the curve-type specificity lives entirely in how t̂ is computed
by the caller.

---

## 8. Matrix Symmetry — Proof
Expand Down Expand Up @@ -345,8 +576,10 @@ For self and near-neighbor elements (|m-n| ≤ 1), singularity extraction is app
## 13. References

- Harrington, R.F. — *Field Computation by Moment Methods* (1968) — MoM formulation; Sections 3-4 directly relevant, cannot be overstated
- Fikioris, G. & Wu, T.T. — "On the Application of Numerical Methods to Hallén's Equation" — exact kernel self-impedance extraction essentially lifted
- Vande Ginste, D. et al. — conformal MoM basis function formulations for curved wires essentially lifted
- Fikioris, G. & Wu, T.T. — "On the Application of Numerical Methods to Hallén's Equation," IEEE TAP, 2001 — exact vs. approximate kernel behavior; confirms exact kernel is bounded (no true singularity); Toeplitz matrix element structure for straight wire
- Mystilidis, C., Papakanellos, P.J., Fikioris, G. & Mavrogordatos, T.K. — "On the Application of Numerical Methods to Hallén's Equation: The Case of a Lossy Medium," IEEE AWPL, 2020 — extends Fikioris & Wu to complex wavenumber (lossy surrounding medium); foundation for deferred Sommerfeld integral work (see Open Question 2)
- Champagne, N.J., Williams, J.T. & Wilton, D.R. — "The Use of Curved Segments for Numerically Modeling Thin Wire Antennas and Scatterers," IEEE TAP, 1992 — foundational CMoM curved-wire formulation
- Rogers, S.D. & Butler, C.M. — "An Efficient Curved-Wire Integral Equation Solution Technique," IEEE TAP, 2001 — exact kernel treatment for curved wires
- Burke & Poggio — *NEC-2 Method of Moments Code* (1981) — thin-wire kernel formulation (the baseline we improve upon!)
- `docs/phase1-geometry/math.md` — parametric forms r(τ), r'(τ), |r'(τ)| used in arc length elements
- `docs/phase2-matrix-fill/validation.md` — validation cases whose expected values are derived from this document
Expand Down
Loading