Bug fix for prism cross section

[oskooi]

Fixes NanoComp/meep#3311.

Symptom

The example in NanoComp/meep#3311 builds a mp.Prism with a 200 × 4 rectangular base in the xy plane, height=2.0, axis=ẑ, sidewall_angle=30°. The yz cross section rendered as a six-side polygon instead of the expected four-sided trapezoid.

Root cause

utils/geom.c:3139, in intersect_line_with_segment — the test for "does the query point q0 lie on the edge q1q2?":

if (fabs(dot) < (1.0 - THRESH) * q01 * q02) return NON_INTERSECTING;   // THRESH = 1e-5

It judged collinearity from the cosine of the angle subtended at q0. For a point at perpendicular distance d from the midpoint of an edge of length L, that cosine is ≈ 1 - 8d²/L², so the test accepted everything with

d < L · sqrt(THRESH/8) = 1.118e-3 · L

An angular tolerance used as a distance test — so the effective slop scales with edge length. node_in_or_on_polygon (geom.c:3204) treats that verdict as "on the boundary" and returns include_boundaries = 1, giving point_in_prism a false-positive halo around the entire region.

For this geometry L = 200 → halo = 0.2236. Measured against libctl by bisection, matching the closed form to 7 digits:

z	true-half-width	predicted true + halo	measured
0	2.000000	2.223607	2.2236074
2	0.845299	1.066324	1.0663248

Confirmed across aspect ratios (edge 200/40/8/4 → 0.22361/0.04472/0.00894/0.00447) and independent of the sidewall angle — identical at 0°, 5°, 30°, 45°, -30°, — and of height.

Why the sidewall angle exposes it, and why six sides

Meep only attributes a point to an object whose bounding box contains it. get_prism_bounding_box correctly unions the bottom and top polygons, giving y ∈ [-2, 2], so the rendered profile is min(2 - z·tan30° + halo(z), 2):

• z ∈ [0, 0.386] — clipped flat at y = ±2 → two spurious vertical walls
• z ∈ [0.386, 2] — follows the halo line 2.2236 — 0.5786·z

bottom + 2 vertical + 2 slanted + top = 6 sides. At sidewall_angle=0 the bounding box clips the halo uniformly and the bug is invisible — the taper doesn't cause it, it just stops the bounding box from hiding it.

The fix

Measure the actual perpendicular distance instead:

#define SEG_PERP_TOL 1.0e-10

double cross = M01 * RHSy - M11 * RHSx; /* (q1-q2) x (q1-q0) */
if (cross * cross > SEG_PERP_TOL * SEG_PERP_TOL * L2 * L2)
  return NON_INTERSECTING;

1e-10 matches the convention already used one line above for the near-parallel test on DetM. Boundary inclusion is perserved — a point exactly on an edge gives cross == 0.

Test

test_slanted_prism_cross_section() in utils/test-prism.c, run on a 200x4 and a 4x4 base — 45 assertions each. The 20 point probes are the discriminating ones: at five heights, a point 1e-3 inside the analytic wall w(z) = 2 - z·tan30° must be inside and one 1e-3 outside must be outside. 10/45 fail per geometry before the fix, 0 after.

Two existing octagon test points needed correcting (127.697→127.6965, 101.824→101.8242): they are 6-digit-rounded edge midpoints sitting ~1.5e-4 outside the polygon, and only read as on-boundary because of the old slop. Items 6 and 7 already used full precision for the same construction.