Alpha fold

[unalmis]

Summary

This PR takes the simple implementation of PlasmaControl#2220 and adds folded-alpha handling for Gamma_delta and Gamma_alpha.
The core idea is that bounce wells are found on long field lines spanning
multiple field periods. The new code
maps long-field-line wells back to effective single-period alpha samples, then
applies the Gamma_delta and Gamma_alpha reductions on that folded alpha
grid.

The implementation keeps the existing Gamma_c Velasco path unchanged.
Only Gamma_delta and Gamma_alpha call the folded-alpha path.

1. Folding Long-Field-Line Wells Into Alpha Samples

For each bounce well with endpoints (zeta_1, zeta_2), the code assigns the
well to a field-period segment using the midpoint approximation

zeta_m = (zeta_1 + zeta_2) / 2,
k = floor(NFP * zeta_m / (2*pi)).

This is implemented in _well_field_period(...).

For a base field-line label alpha_j, a well in field-period copy k
corresponds to the effective single-period label

alpha_eff = alpha_j + iota * k * (2*pi / NFP) mod 2*pi.

This follows the field-line map logic: moving along the long field line by one
field period shifts the equivalent single-period alpha label by
iota * 2*pi / NFP.

The folded representation is indexed conceptually as

(alpha_eff, pitch, local_well)

with optional surface or batch axes in front. A boolean mask records which
folded entries correspond to real wells.

2. Why Local Well Ranking Is Needed

The original well axis is sorted along the extended field line. After folding,
that is not the right branch index. We need the well index local to each
single-field-period copy.

Example:

global well axis:     w = 0    1    2    3    4    5
field period:         k = 0    0    1    1    2    2
local well in period: r = 0    1    0    1    0    1

_local_well_rank(...) computes this local index r. This keeps one physical
local branch coherent across the folded alpha grid. Without it, a single local
well branch could be split across different global long-line well slots, which
would make the alpha scans in Gamma_alpha inconsistent.

3. Gamma_delta Reduction

The threshold is stored in tangent space, so the outward candidate condition is
implemented without forming arctan:

radial > thresh * abs(poloidal).

For folded alpha data, the reduction computes

sum_well [ integral d alpha v_tau(alpha, well) ]
         * H(there exists an outward alpha candidate for this well).

In code:

weight = _alpha_weights(alpha, mask)
v_tau = (v_tau * weight).sum(-3)
has_alpha_out_candidate = alpha_out_candidate.any(-3)
return (v_tau * has_alpha_out_candidate).sum(-1)

The mask remains essential for alpha quadrature weights, but it is not needed
inside has_alpha_out_candidate: dummy wells have zero radial and zero
poloidal, so the strict candidate condition cannot be true for them.

4. Gamma_alpha Reduction

Gamma_alpha uses two signed branch scores:

alpha_out_candidate = radial - thresh * abs(poloidal),
alpha_in_candidate  = -radial - thresh * abs(poloidal).

Positive alpha_out_candidate means the radial drift is sufficiently outward.
Positive alpha_in_candidate means the radial drift is sufficiently inward.

The loss interval orientation depends on the sign of the poloidal drift:

if poloidal >= 0:
    interval starts at alpha_in and stops at alpha_out
else:
    interval starts at alpha_out and stops at alpha_in

The edge-case logic is:

if both alpha_out and alpha_in branches exist:
    use the interval returned by _LossCone
if alpha_out exists but alpha_in does not:
    count the whole alpha domain as lost
if alpha_out does not exist:
    count nothing

The folded path integrates

sum_well integral d alpha v_tau(alpha, well) * loss_cone(alpha, well).

5. Loss-Cone Indicator And Fractional Cells

_LossCone builds periodic intervals from two signed scores:

S_start(alpha), S_stop(alpha).

A start crossing occurs when S_start crosses from nonpositive to positive.
For adjacent samples with previous score p <= 0, current score s > 0, and
spacing h, linear interpolation gives

S(alpha_prev + t h) = p + t (s - p),
t = -p / (s - p).

Equivalently, the distance backward from the current sample to the root is

offset = h * s / (s - p).

That is the formula used by _root(...) and _root_nonuniform(...).

After roots are found, each start root is paired with the nearest downstream
stop root using forward periodic distance:

d(alpha_start, alpha_stop) = (alpha_stop - alpha_start) mod 2*pi.

For order=1, the output is a fractional generalized indicator, not just a
boolean mask. Each alpha sample receives

measure(cell_i intersect loss_interval) / measure(cell_i).

For folded nonuniform alpha samples, cell_i is a periodic Voronoi cell. The
neighbor distances and cell widths are computed by _periodic_voronoi_widths.

The final alpha contribution is

v_tau(alpha_i) * loss_fraction_i * alpha_weight_i.

This improves the placement of the discontinuous loss-cone boundary without
interpolating v_tau.