GeneralizedGradients.jl API Reference

const MAXTOT

Maximum total degree p + q of the monomials x^p y^q tabulated in tables/gg_coef_table.jl, and the highest derivative order nd the table carries for a(m,nd), b(m,nd) and bs(nd). This table was built with MAXTOT = 13.

It is a property of the table rather than a fit setting: it is the ceiling on GGFitInputParams.nd_max, which selects from what is tabulated here.

Covering a higher order means rebuilding the table, which is what programs/create_field_to_gg_coef_tables.jl is for – it writes both of the files in tables/ from one computation of the expansion:

MAXTOT=16 MMAX=17 julia programs/create_field_to_gg_coef_tables.jl

MAXTOT and MMAX are that program's only inputs. It reads each from the environment variable of the same name, falling back to the default built into the program, and writes both back into the table it generates – as this docstring and as the header comment.

Raising MAXTOT enlarges the table in every direction at once: how far the phi recurrence is carried, the (MAXTOT+1)(MAXTOT+2)/2 monomials per GG function, the nd range, and the internal sizes derived from it. Generation time grows steeply as a result, each step costing considerably more than the last. There is little reason to raise it for its own sake – a fit is limited by what the field grid supports long before it is limited by the table.

const MMAX

Maximum multipole order m of the a_m and b_m functions in tables/gg_coef_table.jl: the largest m appearing in a tabulated a(m,nd) or b(m,nd) key. This table was built with MMAX = 14.

It is the ceiling on GGFitInputParams.m_max: a larger m_max is clamped to MMAX, since that is the highest order the table can supply.

bs(nd) carries no multipole order – it is the derivative tower of a_0 – so it is unaffected by this limit.

Changing it means rebuilding the table; see MAXTOT for the command and for what the rebuild costs. MMAX may not exceed MAXTOT + 1: the seed series phi_0 and phi_1 are truncated at x^(MAXTOT+1), so a higher order would be dropped as the table was built, leaving it quietly incomplete rather than visibly wrong. create_field_to_gg_coef_tables.jl rejects that combination instead of writing such a table, so raising MMAX usually means raising MAXTOT with it.

mutable struct FieldGridTable{T}

Holds an electric and/or magnetic field sampled on a 3D grid.

magnetic and electric are 3D OffsetArrays whose elements are field 3-vectors: magnetic[ix,iy,iz] * scale == [Bx, By, Bz] (and likewise [Ex, Ey, Ez]). The grid indices (ix, iy, iz) need not start at 0 or 1; a grid point is at position r0 + dr .* (ix, iy, iz) relative to the anchor.

Eventually this structure may be used for other purposes than generalized gradient modeling. For generalized gradients fitting in particular, the field must be static (RF_frequency = 0).

Fields

• magnetic::OffsetArray{Vector{T}} — magnetic field 3-vectors [Bx, By, Bz] [T].
• electric::OffsetArray{Vector{T}} — electric field 3-vectors [Ex, Ey, Ez] [V/m].
• r0::Vector{T} — grid origin offset (x0, y0, z0) [m].
• dr::Vector{T} — grid spacing (dx, dy, dz) [m].
• g_ref::T — curvilinear-coordinate bending strength 1/bending_radius, in 1/m (0 for a straight reference curve).
• scale::T — Factor to scale the field by. Actual-field = struct-field * scale. Default is 1.0
• RF_frequency::T — RF frequency in Hz (0 for a static field).
• RF_phase::T — RF phase [rad].
• anchor_pt::GridAnchorPt.T — grid anchor point: Beginning, Center, or End. Default is GridGeometry.Center.
• geometry::GridGeometry.T — grid geometry: GridGeometry.XYZ (only possibility at present).

FieldGridTable() builds an empty table with T = Float64; read one from a file with read_field_grid_hdf5.

mutable struct GGFit

Holds the result of a gg_calc_fit fit: the fitted generalized-gradient (GG) coefficient functions sampled at the base planes plus per-plane diagnostics. Returned by gg_calc_fit and consumed by gg_show_fit_results and write_gg_fit.

Fields

• z_base::Vector{Float64} — z position of each base plane [m].
• params::Vector{Tuple{Symbol,Int,Int}} — list of fitted unknowns as (type, m, nd) tuples, where type is one of :a, :b, :bs (bs uses m = 0).
• a::Dict{Tuple{Int,Int},Vector{Float64}} — fitted a(m,nd) functions, (m,nd) => values_over_planes.
• b::Dict{Tuple{Int,Int},Vector{Float64}} — fitted b(m,nd) functions, (m,nd) => values_over_planes.
• bs::Dict{Int,Vector{Float64}} — fitted bs(nd) functions, nd => values_over_planes.
• rms_weighted_plane::Vector{Float64} — weighted RMS fit residual at each base plane, sqrt(Σ w·δ² / Σ w) over the points of that plane's fit region.
• rms_unweighted_plane::Vector{Float64} — RMS fit residual at each base plane over the same points as rms_weighted_plane but with all point weights set to
  1. Equal to rms_weighted_plane when core_weight = outer_plane_weight = 1.
• field_ave_plane::Vector{Float64} — average field magnitude |B| over the fitted grid points of each base plane [T]. Unweighted, and taken from the base plane alone (not the added planes), so it gives the field profile along z and a scale against which rms_weighted_plane can be judged.
• fit_radius_max::Float64 — radius about origin the fit was restricted to [m], or 0 if every grid point was used. Recorded because the residuals and the field contributions above cover that region only, so reading them, or measuring anything else against the fit, needs to know where the fit applies.
• m_max::Int — highest multipole order m retained by the fit. Lower than the cutoff the scan chose if pruning removed every function at the top orders.
• nd_max::Int — highest derivative order nd retained by the fit.
• scan::Vector{GGFitScanPoint} — one GGFitScanPoint per (m_max, nd_max) combination tried, in the order tried. Empty when no scan was requested.
• pruned::Vector{Tuple{Symbol,Int}} — GG functions dropped for producing negligible field, as (:a, m) / (:b, m) / (:bs, 0) tuples. Empty unless prune_ave_limit or prune_max_limit was set. A pruned function is absent from a/b/bs entirely — the dicts are sparse in m, and every consumer treats a missing key as an identically zero function. Not stored in the HDF5 file: which functions are present is already evident from the keys that were written.
• g_ref::Float64 — reference-frame bending strength = 1/bending_radius, in 1/m (0 for a straight reference frame).
• origin::Vector{Float64} — (x, y) line about which the GG coefficients are computed.
• dz_grid::Float64 — spacing between base planes [m].
• eval_plan::Union{Nothing,GGEvalPlan} — internal cache: the compiled GGEvalPlan used by field_and_potential_evaluate_at, built lazily on first evaluation. Not part of the fit data (not serialized); assumes the other fields are not mutated afterward.

mutable struct GGFitInputParams

Input parameters controlling a gg_calc_fit fit. See the documentation of gg_calc_fit for more documentation.

Fields

• origin::Vector{Float64} — Defines the line [x0, y0, z] about which the generalized gradient coefficients are computed. If g (1/bending_radius) is non-zero, origin must be [0, 0].

• n_planes_add::Int — Number of z-planes added to either side of the base z-plane to be used in the analysis of the derivatives at any given base z-plane. For example, for n_planes_add = 2, two planes would be added to either side of the base plane making the total number of planes used in the analysis equal to five.

• core_weight::Float64 — Merit function weight for "core" points (field table points whose transverse (x,y) position is near (0,0)). Default is 1.0 which gives an equal weight for all points of a given z-plane.

• outer_plane_weight::Float64 — Default is 1.0. Merit function weight for z-planes away from the base z-plane when n_planes_add is non-zero.

• m_max::Union{Int,AbstractVector{Int}} — Maximum multipole order m of the a_m/b_m functions used in a fit. m_max can be a single integer or a range. For example, m_max = 2:8 makes gg_calc_fit try each value in turn and keep the one that fits best (see fit_criterion). The default is 4:8. The bs(nd) (that is, a_0 derivative) unknowns carry no m and are never removed by m_max.

• nd_max::Union{Int,AbstractVector{Int}} — Maximum derivative order nd used in the fit. nd_max can be a single integer or a range. For example, nd_max = 2:8 makes gg_calc_fit try each value in turn and keep the one that fits best (see fit_criterion). The default is 3:7

• nd_max_for_m::Dict{Int,Int} — Per-multipole order override of nd_max, mapping a multipole order m to the highest derivative order kept for a_m and b_m. Key 0 sets the limit for b_s, which carries no multipole order. This is used to speed up fitting.

  For example: nd_max_for_m = Dict(4 => 2, 5 => 1) maps a_4/b_4 to max nd = 2 and a_5/b_5 to max nd = 1. Generally a good rule of thumb is that m + max nd should be roughly constant.

• fit_radius_max::Float64 — Fit only the field points that are within this radius transversely of the GG expansion axis origin (not necessarily (0, 0)). 0 meters (the default) uses every point of the field table. Vetoing points outside of some radius is useful if the field table is inaccurate at large radius or a fit at large radius is not needed since this is outside of where particles will travel.

• fit_criterion::Symbol — How a scan picks its winner. Each candidate model is given a score and the lowest score wins. Possible settings are:

    :aic    # Akaike information criterion
    :bic    # Bayesian, and the default

  See the gg_fit documentation for more details

• exclude_functions::Vector{Tuple{Symbol,Int}} — GG functions to leave out of the fit entirely, named as (:a, m), (:b, m) or (:bs, 0) tuples — the same form the pruned field of the result uses.

  exclude_functions = [(:a, 2), (:a, 4), (:b, 2), (:b, 4), (:bs, 0)]

  An excluded function is dropped when the list of unknowns is assembled, so it never gets a design-matrix column: it is not fitted, cannot influence the coefficients that are kept, and is absent from the result. b_s carries no multipole order, so the m of a :bs entry is ignored. Naming a function the model does not contain anyway is a harmless no-op.

  Use this where the answer is known in advance — a magnet whose symmetry forbids the even multipoles, say — since m_max can only cut at the top of the range while this removes orders from anywhere in it. Where the answer is not known in advance, prune_ave_limit/prune_max_limit below decide the same question from the fit itself.

• prune_ave_limit::Float64, prune_max_limit::Float64 for a given fit, these parameters veto GG functions that produce a negligible field. A function here is a whole a_m, a whole b_m, or b_s with all of the function's derivative orders/

  Both limits are fractions of the field table's mean |B|. prune_ave_limit is compared against the function's average contribution and prune_max_limit against its largest. A function is vetoed if it fails either test. A limit of 0 (the default for both) switches that test off. Example:

  • prune_ave_limit = 1e-4 # drop if the average contribution is under 0.01% of <|B|>
  • prune_max_limit = 1e-3 # ... and the largest contribution is under 0.1% of <|B|>

  Pruning is applied after a scan has picked its (m_max, nd_max) winner, and the surviving functions are then refit.

• output_file::String — Name of the output file.

struct GGFitScanPoint

One row of a gg_calc_fit (m_max, nd_max) scan: the model tried and how it scored. Collected in the scan field of the returned GGFit and printed by gg_show_fit_results.

Fields

• m_max::Int, nd_max::Int — the model this row is for.
• n_coef::Int — number of fitted coefficients per base plane.
• rms_weighted::Float64 — weighted RMS residual pooled over all base planes, sqrt(Σ w·δ² / Σ w) over every fitted field-component value.
• rms_weighted_comp::NTuple{3,Float64} — the same quantity restricted to one field component at a time, as the 3-tuple (Bx, By, Bs). Each entry is normalized by that component's own weight sum, so the three are directly comparable to each other and to rms_weighted, which is their weighted quadrature mean. A fit that is bad in only one component shows up here and nowhere else.
• rms_unweighted::Float64 — the same residual with all point weights set to 1.
• score::Float64 — value of the selection criterion; the scanned model with the lowest score is the one gg_calc_fit returns. It is N*ln(RSS/N) plus a penalty of 2 (:aic) or ln(N) (:bic) per fitted coefficient, where RSS is the pooled weighted sum of squared residuals and N the number of fitted field-component values. Only differences between candidates are meaningful — the absolute value is not. See the fit_criterion entry of GGFitInputParams.
• criterion::Symbol — which criterion score was computed with: :bic or :aic. Scores from different criteria are not comparable.

eval_plan(fit::GGFit) -> GGEvalPlan

Return fit's compiled GGEvalPlan, building it and caching it in fit.eval_plan on first use. Assumes fit is not mutated after the first evaluation (a stale plan is not detected).

The plan is the object every evaluator takes: build it once with eval_plan and pass it to field_and_potential_evaluate_at, potential_evaluate_at or field_evaluate_at (there is no GGFit method — evaluating from the plan is what keeps every call type-stable and allocation-free). The plan is also the GPU-ready, Adapt-able evaluation object: for GPU tracking pass eval_plan(fit) (not the fit itself, which holds Dicts and cannot cross to the device) as the field-evaluation parameters, then call the evaluators inside the kernel.

plan = eval_plan(fit)
A, dA = potential_evaluate_at(plan, x, y, s)

field_and_potential_evaluate(fit, ip::Integer, x::Real, y::Real) -> (B, A, dA)

Main entry point. Evaluate the field, vector potential and the Jacobian of A at grid plane ip and transverse position (x, y).

• fit — the GGFit struct returned by read_gg_fit.
• ip — 1-based plane index into fit.z_base.
• x, y — absolute transverse coordinates. fit.origin is subtracted internally to obtain the position relative to the GG expansion axis (the coordinate the expansion is written in). Pass an origin of (0,0) — or use the default — for axis-relative input.

Returns (B, A, dA) where

B  = [Bx, By, Bs]
A  = [Ax, Ay, As]
dA = 3x3 matrix, dA[i,j] = ∂A_i/∂u_j  with  (A_1,A_2,A_3) = (Ax,Ay,As)
     and (u_1,u_2,u_3) = (x,y,s).

field_and_potential_evaluate_at(plan::GGEvalPlan, x, y, s) -> (B, A, dA)

Evaluate the field, vector potential and Jacobian of A at an arbitrary (x, y, s) point, given the compiled evaluation plan.

Obtain the plan once from a fit with eval_plan and reuse it:

plan = eval_plan(fit)
B, A, dA = field_and_potential_evaluate_at(plan, x, y, s)

Taking the plan (rather than the fit) is what makes evaluation fast: this method is @inline, allocation-free (stack-resident SVector scratch), and generic over the coordinate type, so it can be called inside a GPU kernel on an Adapt.adapt-ed plan and with Float32 or ForwardDiff.Dual coordinates.

The GG coefficients are stored only at the grid planes fit.z_base, but the fit gives, at each plane, the whole derivative tower of every GG function: a(m,0..N), b(m,0..N), bs(0..N) with a(m,nd) = dⁿᵈaₘ/dsⁿᵈ and N the maximum order. So for an s between two planes z_L, z_R we have, for each function f, the value and its first N s-derivatives at both ends — 2(N+1) data — which fix a unique two-point Hermite polynomial H(s) of degree 2N+1. Each interpolated derivative is taken from the SAME polynomial, a(m,nd)(s) = H_aₘ⁽ⁿᵈ⁾(s), so the tower stays self-consistent: the interpolated a(m,1) is exactly d/ds of the interpolated a(m,0), etc. The plan compiles these Hermite polynomials once (see low_level.jl).

This is more accurate than independent per-order interpolation (error O(h^{2N+2}) for the base coefficient, using only the two straddling planes) and, because the orders are mutually consistent, the ∂A/∂s that field_and_potential_evaluate forms by bumping a(m,nd) → a(m,nd+1) equals the true s-derivative of the interpolated field. The curl identity B = ∇×A holds at s as before.

• plan — the compiled GGEvalPlan from eval_plan(fit).
• x, y — absolute transverse coordinates (the fit's origin subtracted internally).
• s — absolute longitudinal coordinate.

Returns (B, A, dA) with the same values as field_and_potential_evaluate, but as stack-allocated StaticArrays: B, A are SVector{3,T} and dA is an SMatrix{3,3,T} where T is the promoted coordinate type (Float64 for the usual Float64 inputs). They index like ordinary vectors/matrices (A[1], dA[1,3]).

field_coefficients_at_plane(fit, ip::Integer) -> (CBx, CBy, CBs)

Field-expansion coefficients at a grid plane.

• fit — the GGFit struct from read_gg_fit.
• ip — 1-based plane index into fit.z_base.

Returns (CBx, CBy, CBs); each is a matrix with CB[i+1, j+1] = CB_{c,i,j}, the coefficient of xⁱ yʲ in that field component at the plane.

field_coefficients_at_s(fit, s::Real) -> (CBx, CBy, CBs)

Field-expansion coefficients at an arbitrary s, via the same Hermite interpolation of the GG quantities used by field_and_potential_evaluate_at. Returns (CBx, CBy, CBs) where each CB is a matrix with CB[i+1, j+1] = CB_{c,i,j}, the coefficient of xⁱ yʲ in that field component at the plane.

field_evaluate_at(plan::GGEvalPlan, x, y, s) -> B

Like field_and_potential_evaluate_at but returns only the magnetic field B = [Bx, By, Bs] as an SVector{3,T} for the promoted coordinate type T, skipping the vector potential A and its Jacobian. B is identical to that of the full evaluator. Allocation-free, GPU-capable and type-generic; get plan from a fit with eval_plan(fit). See field_and_potential_evaluate_at for the (x, y, s) conventions.

gg_calc_fit(field::FieldGridTable, params::GGFitInputParams,
            fit_at::Union{Nothing,Tuple{Int,Int}} = nothing) -> GGFit

Fit a 3D DC magnetic field grid to generalized-gradient (GG) coefficients a_m(z), b_m(z), b_s(z) and their z-derivatives, plane by plane. A "plane" here is always a plane at constant z.

The returned GGFit holds the fitted coefficients and per-plane diagnostics. Use gg_show_fit_results to print a summary and write_gg_fit to save the result to an HDF5 file (readable by read_gg_fit).

See examples/run_gg_fit.jl for a complete, runnable example.

Arguments

• field::FieldGridTable — The field grid and associated parameters.
• params::GGFitInputParams - Input fit parameters.
• fit_at::Union{Nothing,Tuple{Int,Int}} — Optional (m_max, nd_max) override. When given, params.m_max and params.nd_max are ignored and the scan is done at this one point only. Everything else in params still applies. This is the way to pick a particular fit out of a scan that has already been run: run the scan once, read the scan table printed by gg_show_fit_results, then refit at whichever (m_max, nd_max) row is wanted without editing p.

How the fit works

The GG coefficients are computed at the equally spaced z-positions coincident with the input field-table planes. These planes are sometimes called "principal planes". The principal plane where a fit is being done is called the "base plane". The fit is done plane by plane: the coefficients at a given base plane are computed independently of the calculation of the coefficients at every other plane, and all coefficients of a given plane are solved for simultaneously by minimizing a merit function

Merit = Σₖ  weightₖ · (field_from_tableₖ - field_from_GG_coefsₖ)^2

The sum runs over all field points lying in a plane within n_planes_add of the base plane. For example, n_planes_add = 2 adds two principal planes on either side, so five planes are used in total. Near the ends of the table the count is reduced — a base plane at the very end of the table uses only three planes when n_planes_add = 2.

The merit function is a linear equation in the GG coefficients so each base plane is solved by weighted linear least squares fit over all field points lying within n_planes_add planes of the base plane.

Adding extra planes smooths the computed values between planes at the cost of making the fit at the principal planes worse. So adding more planes can give a worse fit.

Scanning over the fit size

Convention: m denotes the multipole order for GG functions a_m and b_m while for a(m, nd), b(m, nd), and bs(nd), the nd here denotes the derivative order.

A single fit to the field is done using GG coefficients up to some maximum multipole order and some maximum derivative order. A "scan" is a series of fits using differing maximum multipole order and differing maximum derivative order. As discussed below, The "best" fit is chosen based upon the goodness of fit and how many GG coefficients are needed for the fit.

The multipole orders used in a scan is set by p.m_max which can be an integer if only one order is to be used, or can be an integer array where all values of the array will be used in the scan. Similary, p.nd_max determines the range of maximum derivative orders. For example:

p.m_max  = 1:10       # try m_max = 1, 2, … 10
p.nd_max = 2:6        # crossed with nd_max = 2, 3, … 6

This runs 66 fits. Note: For a given fit, the same maximum multipole order and same maximum derivative order is used in the GG calculation for all planes.

To save time, the number of fits used in a scan can be fine tuned using nd_max_for_m which is a per-multipole override of nd_max, as a Dict mapping m to a derivative limit. For example: nd_max_for_m = Dict(4 => 2, 5 => 1) maps a_4/b_4 to max nd = 2 and a_5/b_5 to max nd = 1. Generally a good rule of thumb is that m + max nd should be roughly constant.

Excluding and pruning

A magnet with a symmetry may need only odd numbered multipole orders and not even. The exclude_functions parameter can be used to exclude GG functions that are not needed. Example: excludefunctions = [(:a, 2), (:a, 4), (:b, 2), (:b, 4), (:bs, 0)] Note `bshas no multipole order, so themof a:bs` entry is ignored.

prune_ave_limit = 1e-4      # against the function's average contribution
prune_max_limit = 1e-3      # against its largest contribution

Pruning is applied after a scan has picked its (m_max, nd_max) winner, and the surviving functions are then refit.

Weighting

The weight of a field point at (x, y) and plane offset dz (relative to the base plane) is the product of a transverse and a longitudinal factor:

weight(x,y,dz) = w_core(x,y) · w_plane(dz)

The transverse factor is determined by core_weight:

w_core(x,y) = core_weight · rmax^2 / (rmax^2 + r^2 · (core_weight - 1))

where r^2 = x^2 + y^2 and rmax is the maximum r over all points. core_weight = 1 (the default) makes w_core constant; core_weight > 1 favors the core (low-r) points at the expense of points farther out. A better core fit is usually desired since beam particles spend most of their time near the core.

The longitudinal factor is set by outer_plane_weight:

w_plane(dz) = 1 + (outer_plane_weight - 1) · |dz| / dz_max

where dz_max is the largest |dz| at the ends of the fit region. If n_planes_add = 0 (so dz_max = 0 and the expression is singular) w_plane is set to 1. outer_plane_weight = 1 (the default) makes w_plane constant; a value between 0 and 1 weights planes nearer the base plane more than the outer planes.

The fit_radius_max parameter is used to veto field points outside of the given radius in the transverse plane. The origin of the circle is the GG expansion axis origin (not necessarily (0, 0)). 0 meters (the default) uses every point of the field table. Vetoing points outside of some radius is useful if the field table is inaccurate at large radius or a fit at large radius is not needed since this is outside of where particles will travel.

Setting `fit_radius_max also re-scales core_weight, whose profile runs from core_weight on the axis to 1 at the outermost fitted point — with a radius set, that is the radius rather than the grid corner. field_ave_plane and the field contributions that drive prune_ave_limit/prune_max_limit likewise cover the fit region only, so pruning judges a GG function by the field it produces where the fit applies. gg_show_fit_residuals reports the residual split at this radius.

Setting fit_radius_max also moves three things onto the fit region, so that they keep describing the same set of points the residuals do: core_weight, whose profile runs from core_weight on the axis to 1 at the outermost fitted point; field_ave_plane; and the field contributions that prune_ave_limit / prune_max_limit act on, so a GG function is pruned on the field it produces where the fit applies rather than at a corner it was never fitted to.

Note that fit_radius_max bounds the merit function, not the model. The fitted GG functions still evaluate anywhere — including out at the corners, where they are now an extrapolation rather than a fit.

Choosing between fits (fit_criterion)

A fit using more GG coefficients will always have a lower RMS residual. A scan therefore cannot pick its winner by residual alone — it needs a rule that charges for coefficients. The input fit_criterion parameter selects that rule. Each scanned fit gets a score and the lowest score wins.

Possible settings of fit_criterion are:

• :aic score = N * ln(RSS/N) + 2 * k # Akaike information criterion
• :bic score = N * ln(RSS/N) + k * ln(N) # Bayesian, and the default

where

RSS = Σ_planes Σ_points  weight · (field_from_table - field_from_GG_coefs)^2
N   = total number of fitted field-component values, summed over base planes
      = Σ_planes · 3 · (points in that plane's fit region)
k   = total number of fitted coefficients
      = (coefficients per plane) · (number of principal planes)

RSS is the weighted residual — the merit function the fit actually minimized, pooled over every base plane. N counts Bx, By and Bs at each field point of each plane's fit region, so a point used by several base planes is counted once per base plane. k likewise counts the whole fit, since each base plane is solved for its own copy of the coefficients.

These are the standard information criteria for a least-squares fit with unknown error variance. The first term, N·ln(RSS/N), is (up to an additive constant that is the same for every candidate, and so cannot affect the ranking) minus twice the maximized Gaussian log-likelihood; it rewards a smaller residual. The second term is the penalty for fit size. They differ only in what one coefficient costs: 2 for AIC, ln(N) for BIC. Since ln(N) > 2 for any N > 7, BIC always penalizes more heavily than AIC and so selects equal or smaller fits. The scores are large negative numbers whose absolute value is meaningless — only differences between candidates matter.

To read the trade-off quantitatively: adding Δk coefficients improves the BIC score only if it lowers RSS by at least the fraction

1 - exp(-Δk·ln(N) / N)  ≈  Δk·ln(N) / N     (for Δk·ln(N) « N)

with 2 in place of ln(N) for AIC. This makes the practical weakness of these criteria explicit here: a field grid easily gives N of order 10^5, so each extra coefficient need only cut the residual sum of squares by of order ln(N)/N ~ 10^-4 to pay for itself under BIC, and 2/N under AIC. AIC in particular is often too weak to reject anything at that N, which is why :bic is the default.

There is a further caveat specific to this problem. AIC and BIC assume the residuals are independent draws from a common Gaussian. Here the residual is dominated by systematic truncation error — the field the fit cannot represent — which is smooth and strongly correlated from point to point, not noise. The effective number of independent measurements is therefore far below N, and the penalty far weaker than the theory intends. Treat the criteria as a defensible ranking heuristic rather than as a rigorous probability statement, and use the scan table's # coefs and RMS columns to make the final call: the honest signal is usually the point of diminishing returns, where the residual stops dropping appreciably as coefficients are added.

Other parameters

Other parameters in GGFitInputParams:

• origin = [x0, y0] — (x, y) line about which the GG coefficients are computed. If field.g_ref is non-zero, origin must be [0, 0]. Default [0.0, 0.0].
• output_file — name of the output HDF5 file. Default "gg_fit_results.h5".

Side note

In theory, the fitting does not require that the field table be a rectangular grid of equally spaced points. In fact, a set of randomly spaced field points would work. Also there is no fundamental requirement that the fit planes be evenly spaced. It is only for convenience that the gg_calc_fit function require a regularly spaced field table and that the output is at evenly spaced planes.

gg_coefficients_at_plane(fit, ip::Integer) -> (a, b, bs)

Generalized-gradient coefficients at a grid plane.

• fit — the GGFit struct from read_gg_fit.
• ip — 1-based plane index into fit.z_base.

Returns the three GG-function dicts of scalar values at the plane: a and b keyed by (m,nd) with a(m,nd) = dⁿᵈaₘ/dsⁿᵈ, b(m,nd) = dⁿᵈbₘ/dsⁿᵈ; and bs keyed by nd with bs(nd) = dⁿᵈ⁺¹a_0/dsⁿᵈ⁺¹ = dⁿᵈb_s/dsⁿᵈ.

gg_coefficients_at_s(fit, s::Real) -> (a, b, bs)

Generalized-gradient coefficients at an arbitrary s, Hermite-interpolated from the straddling grid planes (the same interpolation used by field_and_potential_evaluate_at). Returns the three GG-function dicts of scalar values, as in gg_coefficients_at_plane.

gg_make_fit_residual_table(gg_fit::GGFit, field::FieldGridTable, plane::Integer;
                           dplane::Integer = 0) -> NamedTuple

Field table minus GG fit over the transverse grid of one plane — the data behind one plane's RMS residual, and what to plot to see the shape of a bad fit.

• plane — 1-based index into gg_fit.z_base.
• dplane — grid planes away from that base plane, so dplane = 0 (the default) is the base plane itself. A non-zero offset evaluates the fit by the same Taylor extrapolation gg_calc_fit used when it fitted that neighbouring plane, so the table is comparable with the stored residual for any n_planes_add.

Returns (; x, y, z, plane, dplane, origin, r_fit, B_table, B_fit, dB). x and y are the absolute grid coordinates (vectors, whatever the grid's own index range) and z the plane's longitudinal position. The three field arrays are indexed [ix, iy, component] with component 1, 2, 3 = Bx, By, Bs, and dB = B_table - B_fit. r_fit is the radius the diagnostics treat as the edge of the fit region: gg_fit.fit_radius_max if the fit set one, otherwise the largest circle inside the grid.

The table covers the whole grid either way — a residual outside the fit region is worth seeing, it just is not something the fit was asked to make small.

r = gg_make_fit_residual_table(gg_fit, field, 12)
surface(r.x, r.y, r.dB[:, :, 1])       # Bx residual over the plane

Note that with n_planes_add > 0 a base plane's stored rms_weighted_plane covers its neighbouring planes too, so it will not equal the RMS of this one table.

gg_show_fit_residuals(gg_fit::GGFit, field::FieldGridTable;
                      planes = eachindex(gg_fit.z_base), detail = Int[],
                      mmax::Integer = 8)

Print what each plane's residual is made of, to separate a fit that is too small from a field table a GG expansion cannot represent.

The residual is field table - GG fit over the transverse grid of a base plane (gg_make_fit_residual_table). planes selects which base planes to report and detail those to additionally print a full harmonic table for. mmax is the highest azimuthal harmonic examined.

The measurements

Where on the plane the residual is. A field grid is rectangular and the GG expansion is a series in r, so the grid's corners stick out well beyond the largest circle the expansion is well posed on — √2 further out on a square grid, where every multipole is at its largest. in-rms is the residual inside that circle (or inside fit_radius_max, when the fit set one) and out% the share of the squared residual coming from outside it. Read this column first: an out% near 100 means the plane's RMS is a statement about the corners and not about the fit, and no amount of model adjustment will fix it because it is the series expansion itself running out of convergence there. That is what fit_radius_max is for — with it set, out% describes points the fit was never asked about, and in-rms is the residual that was actually minimized.

Rough versus smooth. The residual is split into a point-to-point irregular part (rough, estimated from second differences, see _resid_rough) and what is left over (smooth, from rms² = rough² + smooth²), and rough% is the first as a percentage of the plane's rms. Interpolation noise in the field table is rough; a GG term the fit does not have is smooth. A residual that is mostly rough will not improve no matter how many GG terms are added; a residual that is mostly smooth is a model that is missing something. Both are measured over the largest centred block inside the inscribed circle, so the corners cannot dominate the answer.

Which term is missing. For the smooth case, the residual is decomposed into azimuthal harmonics on circles about the GG axis, with the transverse part resolved into B_r and B_θ (see _azimuthal_harmonics). A missing multipole of order m shows up as harmonic m growing as r^(m-1) in B_r/B_θ and r^m in B_s. The top harmonic column names the largest one and the exponent p measured for it; when p matches that expectation, the residual really is a multipole the fit does not have, and raising m_max (or lifting an nd_max_for_m cap, or un-excluding a function) will remove it. An exponent near 0 means the harmonic does not grow with radius and is not a multipole at all.

Whether the table can be fitted at all. ∇·B and ∇×B measure the field table against Maxwell's equations by centred differences, with no reference to the fit. A GG expansion is Maxwellian by construction, so whatever the table violates by is error no model can remove — and a large ∇×B in particular means current inside the grid, which the expansion cannot represent at all.

The catch is that those same differences truncate at O(h²·∂³B), so a rapidly varying but perfectly Maxwellian field also shows a violation, and the raw number means nothing. The columns are therefore ratios to the same statistic computed on the fitted GG field tabulated over the same grid — a field that satisfies Maxwell exactly and has much the same structure, so what it shows is the truncation floor for this grid and this field. 1.0 means the table is as Maxwellian as the differencing can tell; a value of several means a violation the grid is fine enough to see.

d²B/dz² is the second difference of the table along z in tesla, read down the column rather than in absolute terms: a plane whose value stands orders of magnitude above its neighbours' does not belong with them, which is a defect in the map and not something to fit.

Reading them together: a high out% is a fit region the expansion cannot cover, whatever else the row says; rough, with ∇·B/∇×B well above 1, is a table-quality problem; smooth with a harmonic whose exponent matches its order is a model that needs that term; a d²B/dz² far above the neighbouring planes' is a seam or a bad plane in the map.

The check none of this replaces is the direct one: refit with a larger m_max and see whether the plane's residual actually falls. A residual that is missing GG terms drops; one that is not saturates, and the saturated level is the honest floor for that plane.

gg_show_fit_results(gg_fit::GGFit, field::FieldGridTable, params::GGFitInputParams)

Print a human-readable summary of a gg_calc_fit result gg_fit: the fit settings, the (m_max, nd_max) scan table if a scan was run, a per-plane table of the weighted and unweighted RMS residuals alongside the average field magnitude of the plane, and the leading multipoles at the central plane as a quick sanity check.

The scan table breaks the weighted residual out by field component (wRMS Bx, wRMS By, wRMS Bs) next to the pooled wRMS resid. When a fit is poor, the split says whether all three components are equally bad — pointing at a model too small, or at a grid the GG expansion cannot represent — or whether one component alone carries the error.

potential_evaluate_at(plan::GGEvalPlan, x, y, s) -> (A, dA)

Like field_and_potential_evaluate_at but returns only the vector potential A and its Jacobian dA, skipping the magnetic field B. This is the entry point used for tracking; get plan from a fit with eval_plan(fit) and reuse it.

For tracking, only A and dA are needed. The B field is the majority of the per-call work (its monomial expansion has more terms than A's), so skipping it is roughly 1.8x faster than the full evaluator while returning identical A, dA. A is an SVector{3,T} and dA an SMatrix{3,3,T} for the promoted coordinate type T. Allocation-free, GPU-capable and type-generic. See field_and_potential_evaluate_at for the (x, y, s) conventions.

read_field_grid_hdf5(path; index = 1) -> FieldGridTable

Read a Bmad/openPMD field_grid HDF5 file (as written by write_field_grid_hdf5, or by Bmad itself) into a FieldGridTable. The magnetic/electric OffsetArrays are indexed (ix_lo:ix_hi, …) with each element a [Bx,By,Bz] 3-vector; the grid index ranges come from gridLowerBound/gridSize (the grid is not assumed to start at zero). r0 is gridOriginOffset (so a point (ix,iy,iz) is at dr .* (ix,iy,iz) + r0 relative to the anchor). An absent field type is left at the struct default. index selects which grid under /ExternalFieldMesh/ to read (default 1).

read_gg_fit(path::AbstractString) -> (fit::GGFit, meta::NamedTuple)

Load a gg_calc_fit result HDF5 file (written by write_gg_fit). Returns a two-tuple whose first component is a GGFit struct holding the GG coefficient dictionaries a, b, bs (and z_base, m_max, nd_max, rms_weighted_plane, rms_unweighted_plane, field_ave_plane, fit_radius_max, g_ref, origin, dz_grid), and whose second component is a NamedTuple of the associated fit-control metadata (n_planes_add, core_weight, outer_plane_weight). The params field of the returned struct is empty (the unknown list is not stored in the file).

fit, meta = read_gg_fit(path)
fit.a            # Dict{(m,nd) => values_over_planes}
fit.g_ref        # reference curvature

write_bmad_field_grid_element(field::Union{AbstractString,FieldGridTable};
                           ele_name::AbstractString = "fieldmap_ele",
                           output_base::AbstractString = ele_name,
                           field_scale::Real = 1.0,
                           hdf5::Bool = true) -> (ele_file, grid_file)

Create a file with a single Bmad lattice element that uses the field grid table for Runge-Kutta tracking. Also write the field grid to a file in HDF5 or ASCII format.

Two files are written: <output_base>.bmad (the lattice element) and the field grid (<output_base>_grid.h5 or _grid.bmad). The reference-curve bending strength is taken from the field grid's g_ref (= 1/bend_radius): if non-zero the element is written as an sbend, otherwise an em_field. The grid is anchored at the entrance of the element (ele_anchor_pt = beginning) with length L = dz*(nz-1); the field-grid r0 keeps the transverse offset (x0, y0) of the input grid and shifts z so the first grid plane sits at the element entrance.

Input:

• field - Field table or name of an HDF5 field table file.
• ele_name — name of the Bmad lattice element. Default "fieldmap_ele".
• output_base — base path for the two output files. Default ele_name.
• field_scale — overall field scale factor written to the field grid. Default 1.
• hdf5 — if true (the default), write the field grid as an openPMD HDF5 file (<output_base>_grid.h5) instead of a plain-text block.

Output

• (ele_file, grid_file) - Tuple of file names.

write_bmad_gg_fit(fit::GGFit; ele_name, output_base, cutoff) -> lattice_file_path
write_bmad_gg_fit(input::AbstractString; output_base, cutoff) -> lattice_file_path

Convert generalized-gradient (GG) coefficients produced by gg_calc_fit into Bmad gen_grad_map format, producing a Bmad lattice element with the GG map attached. Returns the path of the lattice-element file.

The GG fit is supplied either as a loaded GGFit struct (fit, as returned by read_gg_fit) or as the path to a ggcalcfit HDF5 file (output of write_gg_fit), which is read with read_gg_fit.

Usage

using GeneralizedGradients
write_bmad_gg_fit("gg_fit_result.h5")

From the shell (see programs/run_write_bmad_gg_fit.jl):

julia programs/run_write_bmad_gg_fit.jl <gg_fit_result.h5> [output_base] [cutoff]

Keyword arguments:

• ele_name — name of the Bmad lattice element. Default "gen_grad_ele" (or, for a file input, the output_base basename).
• output_base — base name for the two output files: <output_base>.bmad (the lattice element) and <output_base>_gg.bmad (the attached gen_grad_map). Defaults to the input file name without extension, or ele_name for an in-memory fit.
• cutoff — relative magnitude cutoff for pruning negligible multipole curves. A curve is dropped if its peak |GG| is below cutoff * (largest peak |GG| of any curve). Default 0 (keep every non-zero curve).

The reference-coordinates bending "strength" 1/bend_radius [1/m] is taken from fit.g_ref; non-zero => the element is an sbend with curved_ref_frame = T, otherwise an em_field.

Background: the two GG conventions

This project (Van der Schueren / Sagan) characterizes the field by midplane- derivative generalized gradients a_m(s), b_m(s), b_s(s):

B_x(x,0,s) = Σ_{m≥1} a_m(s) x^{m-1}/(m-1)!     (skew / "cos" family)
B_y(x,0,s) = Σ_{m≥1} b_m(s) x^{m-1}/(m-1)!     (normal / "sin" family)
B_s(0,0,s) = b_s(s)                            (solenoidal, m = 0)

Bmad's gen_grad_map (Venturini-Dragt) instead uses azimuthal-harmonic gradients C_{m,α}(z), α ∈ {sin, cos}, where the field is (Sagan, IPAC23 Eq. 4)

B_ρ = Σ_{m≥1,n} f(m,n)(2n+m) ρ^{2n+m-1}[C^{[2n]}_{m,s} sin mθ + C^{[2n]}_{m,c} cos mθ]
      + Σ_{n≥1} f(0,n)(2n) ρ^{2n-1} C^{[2n]}_{0,c}
B_θ = Σ_{m≥1,n} f(m,n) m ρ^{2n+m-1}[C^{[2n]}_{m,s} cos mθ - C^{[2n]}_{m,c} sin mθ]
B_z = Σ_{m≥0,n} f(m,n) ρ^{2n+m}[C^{[2n+1]}_{m,s} sin mθ + C^{[2n+1]}_{m,c} cos mθ]

with f(m,n) = (-1)^n m!/(4^n n!(n+m)!), sin = normal, cos = skew.

Equating the two on the midplane gives the exact relations used here. For each azimuthal m and derivative order j (with k ≡ m):

C^{[j]}_{m,s} = (1/m!)[ b^{[j]}_m - (m-1)! Σ_{n≥1, m-2n≥1} Wn(m,n) C^{[j+2n]}_{m-2n,s} ]
C^{[j]}_{m,c} = (1/m!)[ a^{[j]}_m - (m-1)! Σ_{n≥1, m-2n≥1} Wc(m,n) C^{[j+2n]}_{m-2n,c}
                                 - (m even) (m-1)! Us(m) b_s^{[m+j-1]} ]
C^{[j]}_{0,c} = b_s^{[j-1]}                                            (j ≥ 1)

Wn(m,n) = (-1)^n (m-2n)!(m-2n)/(4^n n!(m-n)!)     (normal radial mixing)
Wc(m,n) = (-1)^n (m-2n)! m   /(4^n n!(m-n)!)       (skew radial mixing)
Us(m)   = (-1)^{m/2} m /(4^{m/2} ((m/2)!)^2)        (skew↔solenoid coupling)

where x^{[j]} ≡ dʲx/dsʲ is supplied directly by the fit (a[(m,j)], b[(m,j)], bs[j]). These recursions are solved in order of increasing m, reusing the lower-m towers. Truncation at the fit's maximum derivative order nd_max bounds the radial-correction sums exactly as the fit itself is bounded, so the resulting gen_grad_map reproduces the project field to machine precision.

Output

A Bmad gen_grad_map (field_type = magnetic) attached to a lattice element. As with grid fields, the map is anchored at the entrance of the element (ele_anchor_pt = beginning), z-positions run 0, dz, 2dz, … and the element length is L = (n_planes - 1) * dz. The transverse anchor r0 is the GG expansion axis (origin). For a curved reference (g_ref ≠ 0) the element is an sbend with g = g_ref and curved_ref_frame = T; otherwise it is an em_field.

write_field_grid(output_file::AbstractString, fg::FieldGridTable)

Write a FieldGridTable. If output_file ends in .h5/.hdf5 it is written as a Bmad openPMD field_grid HDF5 file (write_field_grid_hdf5, readable by Bmad); otherwise it is written as a Julia source file (like ags-snakes/wsnk_fieldmap.jl) that defines fg when included.

write_field_grid_hdf5(hdf5_output_file::AbstractString, fg::FieldGridTable)

Write a FieldGridTable as an openPMD HDF5 field_grid file matching Bmad's hdf5_write_grid_field with geometry = xyz.

Input

• hdf5_output_file – Output file name.
• fg – Field grid table.

write_gg_fit(gg_fit::GGFit, field::FieldGridTable, params::GGFitInputParams) -> output_file_path

Write a gg_calc_fit result gg_fit to an HDF5 file (readable by read_gg_fit).

Stores the fitted GG coefficients plus enough metadata to reproduce and interpret the fit later. The (large) input field table is NOT stored. The file is written to params.output_file and its path is returned.

HDF5 schema

root datasets   : z_base, rms_weighted_plane, rms_unweighted_plane,
                  field_ave_plane, origin                     (Float64[])
root attributes : g_ref, dz_grid (Float64); m_max, nd_max, n_planes_add (Int);
                  fit_radius_max, core_weight, outer_plane_weight (Float64)
groups a, b     : m (Int[]), nd (Int[]), values (Float64[nkeys, nplanes])
                  -- reconstruct Dict{(m,nd) => values[i,:]}
group  bs       : nd (Int[]), values (Float64[nkeys, nplanes])
                  -- reconstruct Dict{nd => values[i,:]}

GGEvalPlan{VF,TWS,CPS,NG,NP,NQ}

Compiled, type-stable evaluation plan built once per fit (see lowlevel.jl). Holds the interpolation towers (an NTuple of [`Tower](@ref)) and the per-component monomial term listscomps(orderBx By Bs Ax Ay As dAx dAy dAs`).

The scratch-sizing constants are carried as Val fields (ng = ngvals, np = pmax+1, nq = qmax+1) rather than plain integers so that (a) they are compile-time constants inside the evaluator — sizing the stack-resident SVectors with no heap allocation — and (b) they pass through Adapt.adapt unchanged. Every array field is generic over its backing type, and the whole plan is Adapt.@adapt_structured, so adapt(CuArray, plan) (or whatever backend Adapt targets) yields a plan whose evaluation runs inside a GPU kernel. The GG value getters (_interp_gvals) and component evaluators (_comp_value / _comp_full) are all allocation-free and generic over the coordinate type, so the same plan tracks in Float64, Float32, or ForwardDiff.Dual.

Fields

• origin::NTuple{2,Float64} — (x, y) line the GG expansion is written about.
• z::VF — base-plane positions [m], in increasing order.
• towers::TWS — the interpolation towers, an NTuple{NT,_Tower} (one per multipole m present, plus the single bs tower).
• comps::CPS — per-component monomial term lists, an NTuple{9,_CompTerms} in the order Bx By Bs Ax Ay As dAx dAy dAs.
• ng::Val{NG} — number of gvals slots, as a compile-time constant.
• np::Val{NP} — pmax + 1, the x power-table length.
• nq::Val{NQ} — qmax + 1, the y power-table length.

_CompTerms{VI,VF}

One output component's monomial terms. Evaluating the component accumulates value += Σ w[t] * gvals[slot[t]] * x^p[t] * y^q[t] over all terms t.

Generic over its backing-array types (VI for the integer arrays, VF for the weights) so the same struct is Vector-backed on the host and device-array backed after Adapt.adapt — see GGEvalPlan.

Fields

• slot::VI — index into gvals of the GG value each term multiplies.
• w::VF — coefficient of each term.
• p::VI — power of x of each term.
• q::VI — power of y of each term.

_Tower{VI,VF,MF}

One GG derivative tower (a fixed multipole m, or the single bs tower). poly[d+1, pair] is the coefficient of u^d (with u = s - zref[pair]) of the interpolant on plane-pair pair; interpolating gives H⁽ⁿᵈ⁾(s) for the tower's orders nd = 0..N, scattered into gvals at slots[nd+1]. Non-contiguous orders (nd > N) are taken from the nearest (left) plane via extra_vals[e, pair].

Generic over its backing-array types (VI/VF/MF for the integer vectors, float vectors and float matrices) so it survives Adapt.adapt to the GPU.

Fields

• N::Int — highest contiguous derivative order the tower interpolates.
• deg::Int — polynomial degree of the interpolant: 2N+1 for Hermite, N for Taylor.
• slots::VI — gvals slot for each order nd = 0..N.
• poly::MF — interpolant coefficients, (deg+1) x npairs.
• zref::VF — left-plane position of each plane-pair [m], length npairs.
• extra_slots::VI — gvals slots for the non-contiguous orders nd > N (rare).
• extra_vals::MF — value of each extra order at each plane, n_extra x P.

_accum(tdict, valfun, g_ref) -> K

Coefficient-array builder. K[p+1,q+1] = coefficient of xᵖ yᵠ. valfun(key) returns the GG function value multiplying that table entry.

_anchor_from_str(s) -> GridAnchorPt.T

Parse an openPMD eleAnchorPt string into a GridAnchorPt enum value.

_anchor_to_str(a::GridAnchorPt.T) -> String

Map a GridAnchorPt enum value to its openPMD eleAnchorPt string.

_attr(obj, name, default)

Read an attribute if present, else return default.

_azimuthal_harmonics(map, mmax; nth = 64, nrad = 5) -> NamedTuple or nothing

Azimuthal Fourier decomposition of a residual map, on circles centred on the GG expansion axis.

The residual's transverse part is resolved into radial and azimuthal components B_r, B_θ before the transform, because those are what a single multipole makes clean. A missing multipole of order m has a scalar potential going as r^m·sin(mθ), hence

B_r, B_θ  ~  r^(m-1) · {sin,cos}(mθ)        B_s  ~  r^m · {sin,cos}(mθ)

— one azimuthal harmonic m, with a definite radial power. Noise has no preferred harmonic and no radial growth. So both the harmonic and its measured radial exponent have to line up before a residual can be blamed on a missing GG term, and the exponent is the part that is hard to fake.

Returns (; radii, amp, expo), where amp[k][m+1, ir] is the amplitude [T] of harmonic m of kind k (1, 2, 3 = B_r, B_θ, B_s) on the circle of radius radii[ir], and expo[k][m+1] is the exponent of the power law fitted through amp by least squares in log-log, over the outer half of the radii. Returns nothing when no circle centred on the axis fits inside the grid.

The circles stop at map.r_fit, so this sees only the part of the plane the expansion is well posed on; _radial_split covers what happens outside it.

_build_comp(Ta, Tb, Tbs, bump, g_ref, slot_a, slot_b, slot_bs) -> _CompTerms

Flatten one output component's (a, b, bs) monomial tables into a flat term list, folding g_ref^k into each weight and resolving every (m,nd)/nd key to a gvals slot. bump shifts the GG derivative order by one (used for ∂A/∂s). Terms whose GG value is structurally zero (missing key) are dropped.

_build_eval_plan(fit::GGFit) -> GGEvalPlan

Compile fit into a GGEvalPlan: assign a dense gvals slot to every GG value, build the interpolation towers, and flatten the nine output components' monomial tables into term lists. Called once per fit and cached by eval_plan.

_build_towers_mnd!(towers, z, d, slotmap, next) -> next

Assign gvals slots for the keys of an (m,nd)-keyed dict d (a or b), build one _Tower per multipole m, and append them to towers. Returns the next free slot index.

_build_towers_nd!(towers, z, d, slotmap, next) -> next

Like _build_towers_mnd! but for the nd-keyed bs dict: a single _Tower.

_coefsum(terms, x::Float64, y::Float64, g_ref)

CB coefficient sum: Σ coeff·g_ref^k·x^p·y^q over the table entries for one (component, function) — one entry of the CB grids built in gg_calc_fit.

_comp_array(Ta, Tb, Tbs, aval, bval, bsval, g_ref) -> K

Combined coefficient array of a component: sum of its a, b and bs parts. Ta/Tb are keyed by (m,nd) and Tbs by nd.

_comp_full(ct, gvals, xp, yq) -> (val, dvx, dvy)

Like _comp_value but also returns the x and y partial derivatives of the component.

_comp_value(ct, gvals, xp, yq) -> val

Evaluate a component's value Σ w * gvals[slot] * x^p * y^q, where xp[i+1] = x^i and yq[j+1] = y^j are the power tables. Generic over the coordinate type.

_component_dataset(field, c)

Lay component c of a (ix, iy, iz) OffsetArray of 3-vectors out as a 1-based (nx, ny, nz) complex array. HDF5.jl reverses dims on write, so the dataset lands on disk exactly like Bmad's own Fortran writer (H5Screate_simple_f with Fortran dims [nx,ny,nz]): Bmad's reader gets data_dim = (nx,ny,nz) and, with data_order "F", reads the column-major buffer back into pt[ix,iy,iz] correctly.

_contiguous_order(orders) -> N

Largest N such that orders 0,1,…,N are all present in orders (sorted).

_eval_scratch(plan, x, y, s) -> (gvals, xp, yq)

Prepare the per-call scratch for evaluating plan at (x, y, s): interpolate every GG tower onto s (gvals) and build the x/y power tables (xp[i+1] = x^i, yq[j+1] = y^j, with plan.origin already subtracted). All three are stack-resident SVectors, so a call allocates nothing on the heap and runs unchanged on the GPU. Generic over the coordinate type T (so Float64, Float32, and ForwardDiff.Dual all work). Shared by field_and_potential_evaluate_at and potential_evaluate_at.

_expand_coefs(params_list, theta) -> (a, b, bs)

Scatter a solved coefficient matrix (theta[col, plane], rows following params_list) into the a, b and bs dictionaries of a GGFit. Only the unknowns in params_list get keys, so a pruned function is simply absent.

_field_CB(fit, ip::Integer) -> (CBx, CBy, CBs)

Field-expansion coefficients B_c(x,y,s) = Σ_{i,j} CB_{c,i,j}(s) xⁱ yʲ. Returns full _NMAX×_NMAX arrays summed over the a, b, bs parts.

_fit_candidates(spec, table_max) -> Vector{Int}

Resolve an m_max/nd_max input-parameter setting into the list of values gg_calc_fit should try. An Int pins that value; a vector or range is taken as-is. Values are clamped to table_max (the largest the coefficient table supports) and deduplicated, so an over-wide request such as 0:100 does not refit the same top fit repeatedly.

_fit_over_planes(cand_cols, geom) -> NamedTuple

Solve the weighted least-squares fit at every base plane, for every candidate column subset in cand_cols, from one shared design matrix per plane. geom carries the grid, weighting and basis data that does not depend on which columns are fitted (see the call site in gg_calc_fit), so this can be re-run on a reduced column set — as the pruning pass does — without rebuilding any of it.

Returns (; thetas, rmsw_c, rmsu_c, rmsw_comp_c, nrow_pl, wsum_pl, wsum_comp_pl, field_ave_plane): per candidate the fitted coefficients (thetas[c][col, plane], rows following that candidate's column subset), the weighted and unweighted per-plane residuals, and the weighted residual split by field component; then the per-plane row counts and weight sums, and the average |B| of each base plane.

_fit_score(criterion, rss, ndata, nparam) -> Float64

Selection score for one scanned fit; the lowest score wins. rss is the pooled weighted sum of squared residuals over ndata field-component values, and nparam the total number of fitted coefficients (per-plane count times the number of base planes). With RSS = rss, N = ndata and k = nparam:

:aic    N*log(RSS/N) + 2k              Akaike information criterion
:bic    N*log(RSS/N) + k*log(N)        Bayesian information criterion

N in the log term is the count of data values entering the Gaussian log-likelihood, not a weight total.

:aic and :bic share the leading term — minus twice the maximized Gaussian log-likelihood, dropping an additive constant that is common to every candidate and so cannot change the ranking — and differ only in the cost of one coefficient, 2 versus log(N). See the gg_calc_fit docstring for how to read the trade-off and for the caveats that apply when the residual is dominated by systematic truncation error rather than by noise.

_geometry_from_str(s) -> GridGeometry.T

Parse an openPMD gridGeometry string into a GridGeometry enum value.

_geometry_to_str(::GridGeometry.T) -> String

Map a GridGeometry enum value to its openPMD gridGeometry string (only XYZ is supported).

_gg_exclude_set(exclude) -> Set{Tuple{Symbol,Int}}

Validate and normalize a user exclude_functions list into the group form used throughout. b_s carries no multipole order, so any (:bs, m) normalizes to (:bs, 0). Naming a function the model does not contain is a harmless no-op; a type other than :a, :b or :bs is a typo and raises.

_gg_field_contributions(gg_fit::GGFit, field::FieldGridTable) -> (rows, b_ave)

How much field each fitted GG function actually produces, measured over every fitted transverse grid point of every base plane — that is, over gg_fit.fit_radius_max when one is set, and over the whole grid otherwise.

rows holds one (typ, m, ave, max) per GG function — (:a, m) and (:b, m) for each multipole order m present in the fit, plus (:bs, 0) — where ave and max are the mean and the largest |B| [T] that function generates with every other GG coefficient set to zero. All derivative orders nd of the function are included, since they too contribute to the field at the plane. The field is linear in the GG coefficients, so a row is exactly what dropping that function from the fit would remove from the modeled field. b_ave is the mean |B| of the field table itself, the scale the rows are to be read against.

Both gg_show_fit_results and the pruning pass of gg_calc_fit read these rows — which is why they carry the (typ, m) group rather than a printable label.

This is the useful form of "how big is this coefficient": the raw a/b values are not comparable across m, since the basis function each multiplies carries a different power of r and so a different size over the grid.

_gg_group(param) -> (typ, m)

The GG function a (typ, m, nd) unknown belongs to: all derivative orders of one a_m or b_m share a group, and every bs unknown lands in (:bs, 0). This is the granularity at which contributions are measured and functions are pruned.

_gg_label(typ, m) -> String

Printable name of a GG function group: "a_3", "b_5", "b_s".

_gg_nd_caps(nd_max_for_m) -> Dict{Int,Int}

Validate a user nd_max_for_m setting, the per-multipole override of nd_max. Keys are multipole orders (0 being b_s, which carries no order) and values the highest derivative order kept for that multipole. Naming an m the coefficient table does not contain is a harmless no-op; a negative order or a negative limit is a mistake and raises.

_gg_num(x::Real) -> String

Lossless, compact Float64 text: repr emits the shortest string that parses back to the identical Float64 (Bmad's Fortran reader accepts the e-notation). Without this, cancellation in B_s (which is a small difference of larger terms) magnifies the rounding of a fixed-precision format.

_gg_taylor_getters(gg_fit, ip, dz) -> (aval, bval, bsval)

GG coefficient getters for base plane ip, Taylor-shifted to the longitudinal offset dz: aval(m, nd) returns a_m^[nd] at z_base[ip] + dz as the fit models it,

Σ_{j ≥ nd}  a(m,j)[ip] · dz^(j-nd) / (j-nd)!

and likewise for b and b_s. At dz = 0 this is just the stored coefficient.

This is the same Taylor extrapolation the gg_calc_fit design matrix uses to carry a base plane's coefficients onto its neighbouring planes, so a residual built from these getters is the residual the fit actually minimized.

_hermite_derivs(zL, zR, fL, fR, sq) -> Vector

Two-point Hermite tower: fL[j+1] = f⁽ʲ⁾(zL), fR[j+1] = f⁽ʲ⁾(zR), j = 0..N. Returns [H⁽ⁿᵈ⁾(sq) for nd=0..N] where H is the degree-(2N+1) Hermite interpolant. Built via confluent Newton divided differences in the local coordinate u = s - zL (nodes: 0 with multiplicity N+1, hstep with multiplicity N+1).

_hermite_poly(zL, zR, fL, fR) -> Vector{Float64}

Monomial coefficients (in u = s - zL) of the degree-(2N+1) two-point Hermite interpolant with fL[j+1] = f⁽ʲ⁾(zL), fR[j+1] = f⁽ʲ⁾(zR), j = 0..N. Same confluent-Newton construction as _hermite_derivs, but returns the polynomial rather than evaluating it.

_interp2(v, x, y, px, py) -> Float64

Bilinear interpolation of v[ix, iy] sampled on the evenly spaced coordinate vectors x and y, at the point (px, py). Points outside the grid are clamped to its edge; the callers only ask for points inside it.

_interp_gg_fit(fit, s::Real) -> fit::GGFit

Take GG fit results fit which give the GG functions at a set of planes and return a similar GGFit but with one plane: the GG coefficients for that plane are the interpolated GG coefficients at the given s-position.

• fit — GG coefficients for all planes.

Builds a single virtual plane at s by Hermite-interpolating every GG derivative tower from the two straddling grid planes (one-plane Taylor if only one plane).

_interp_gvals(plan, s) -> gvals::SVector

Interpolate every tower onto s, returning the dense GG value vector as a stack-resident SVector{ngvals, T} (T = typeof(s)). Allocation-free and GPU-safe: the accumulator is an immutable SVector updated with Base.setindex (no MArray, no heap), and the power u^(d-nd) is carried incrementally so no scratch buffer is needed. Iterating the NTuple of towers is unrolled by the compiler.

_interp_mnd_dict(d, iL, iR, zL, zR, sq, single) -> Dict

Interpolate an (m,nd)-keyed dict (a, b): build one Hermite per multipole m.

_interp_nd_dict(d, iL, iR, zL, zR, sq, single) -> Dict

Interpolate an nd-keyed dict (bs): a single Hermite tower.

_interp_tower(fL, fR, zL, zR, sq, single) -> Vector

Interpolate one GG function's derivative tower onto sq (Hermite, or Taylor if single).

_is_hdf5_path(path)

True if path should be treated as an HDF5 file (".h5" or ".hdf5" suffix).

_make_dA(Axx, Axy, dAxv, Ayx, Ayy, dAyv, Asx, Asy, dAsv) -> SMatrix{3,3}

Assemble the 3x3 Jacobian dA[i,j] = ∂A_i/∂u_j (rows Ax,Ay,As; columns x,y,s) as a stack-allocated SMatrix (no heap allocation; eltype inferred from the arguments). Arguments are given row-major; the SMatrix constructor takes them column-major.

_make_tower(z, Fm, slots, extra_slots, extra_planevals) -> _Tower

Build one tower from a per-m value matrix Fm[j+1, plane] = f⁽ʲ⁾ at that plane, its slots, and any extras. z are the base planes. Precomputes the interpolant's monomial coefficients on every straddling plane-pair (or the single-plane Taylor series when there is one plane).

extra_planevals[e] is the per-plane value vector of extra order e; it is flattened into the n_extra x P matrix _Tower.extra_vals so the tower holds only plain arrays (adaptable to the GPU; a vector-of-vectors is not).

_maxwell_stats(Bm, B0, Bp, r0, dr, g_ref) -> (divB, curlB, zjump)

Measure a three-plane stack of field values against Maxwell's equations. Bm, B0, Bp are the planes at iz-1, iz, iz+1, each indexed [ix, iy, component], r0/dr the grid origin and spacing.

divB and curlB are the RMS over the middle plane's interior points of |∇·B| and |∇×B| as a fraction of the RMS of the field's own gradient there (sqrt(Σ_ij (∂B_i/∂u_j)²)). Dividing by the gradient makes the numbers dimensionless and comparable from plane to plane, but the absolute level means nothing on its own: the same centred differences that measure the violation truncate at O(h²·∂³B), so a perfectly Maxwellian field that varies rapidly across a cell produces a large ratio too. It has to be read against the same statistic computed on a field known to be Maxwellian — which is what gg_show_fit_residuals does, running this on the fitted GG field tabulated over the same grid and reporting the table's level as a multiple of it.

zjump is the RMS over the plane of |B[iz-1] - 2B[iz] + B[iz+1]| [T], the second difference along z. It is left unnormalized because it is read down the column: a smoothly varying map gives h²·∂²B on every plane, so a plane whose value stands orders of magnitude above its neighbours' is a seam between two maps, a wrong z step, or a dropped plane.

The derivatives are taken in the curvilinear frame set by g_ref (scale factors 1, 1, g with g = 1 + g_ref·x):

∇·B    = (1/g)[∂(g·Bx)/∂x + ∂(g·By)/∂y + ∂Bs/∂s]
(∇×B)x = ∂Bs/∂y - (1/g)·∂By/∂s
(∇×B)y = (1/g)·∂Bx/∂s - (1/g)(g_ref·Bs + g·∂Bs/∂x)
(∇×B)s = ∂By/∂x - ∂Bx/∂y

In a current-free region both vanish for any physical field. A GG expansion is Maxwellian by construction and so can only ever fit the part of a table that obeys them: whatever the table violates them by is error no model can remove. A large ∇×B in particular says the region carries current — a coil inside the grid — which the expansion has no way to represent.

_peak(d, m) -> Float64

Peak |value| of a derivative tower's value column (j = 0), used for cutoffs.

_polyval(K, x, y) -> (val, dvx, dvy)

Value and (x,y) partials of the plain polynomial Σ K[i,j] xⁱ yʲ.

_pow_table(::Val{N}, x) -> SVector{N}

Stack-resident power table xp[i+1] = x^i for i = 0..N-1, as an immutable SVector (no heap, GPU-safe).

_radial_split(map) -> (r_fit, rms, rms_in, out_frac)

Split a residual map by radius about the GG expansion axis, at map.r_fit.

Field grids are rectangular and the GG expansion is a series in r, so the two do not match: the corners of a square grid sit at r = √2 · r_in, where every multipole is at its largest and the series is at its least convergent — a corner point of a ±a grid weighs a term of order m by 2^(m/2) against the same term at r = a. Those points are also the majority of a rectangular grid's area, so a plane's RMS residual can be almost entirely made of them while the fit is good everywhere the expansion is meant to be used. That is the case fit_radius_max exists to remove, and when the fit set one, this splits there instead.

Returns the split radius, the RMS residual over the whole plane, the RMS over the points inside it, and out_frac, the share of the total squared residual contributed by the points outside it.

_read_field_group(g1, name, lb, nx, ny, nz)

Read a field group ("magneticField"/"electricField") into an (ix, iy, iz) OffsetArray of [Bx,By,Bz] 3-vectors indexed from lb, or nothing if absent.

In a Bmad field_grid file each component dataset is written Fortran-order (logical dims [nx,ny,nz]; on-disk C-dims (nz,ny,nx)). HDF5.jl reverses dims on read, so it hands back a 1-based (nx, ny, nz) array that is already the field – no transpose needed.

_resid_rough(d) -> Float64

Scale of the point-to-point irregular part of a residual map d[ix, iy].

For a residual that is white noise of standard deviation σ, the second difference d[i-1] - 2d[i] + d[i+1] has variance 6σ², so mean(D²)/6 recovers σ². For a smooth residual the same second difference is h²·∂²d and is small by two powers of the grid spacing. The estimate is taken along x and along y and the smaller kept, so a residual that varies rapidly in one direction only is not mistaken for noise.

Smooth structure does leak in — h²·∂²d is not zero — so this is an upper bound on the noise, which is the safe direction: it can only understate how much of the residual a larger model could remove.

_taylor_derivs(z0, f0, sq) -> Vector

Single-point Taylor tower: from f and its derivatives at z0, return [P⁽ⁿᵈ⁾(sq) for nd=0..N] with P the Taylor series, i.e. f extrapolated to sq.

_taylor_poly(f0) -> Vector{Float64}

Monomial coefficients (in u = s - z0) of the single-plane Taylor series: poly[d+1] = f0[d+1] / d!, so that P⁽ⁿᵈ⁾(u) matches _taylor_derivs.

_trim3(CBx, CBy, CBs) -> (CBx, CBy, CBs)

Trim three coefficient arrays to the smallest (x,y) extent holding every nonzero entry, so the returned matrices are indexed CB[i+1, j+1] = CB_{c,i,j}.

_val(::Val{N}) -> N

Recover the wrapped integer of a Val field as a compile-time constant (used to size the stack-resident SVector scratch).

_write_coef_group(parent, name, d; single::Bool = false)

Write a Dict keyed by (m,nd) (or by nd, if single) as index arrays + matrix.

_write_field_component_jl(io, name, field)

Write the fg.<name> OffsetArray of [Bx,By,Bz] 3-vectors as include-able Julia.

_write_field_grid_text(path, mag, r0, dr, is_bend, field_scale)

Write the plain-text field-grid block from an (ix, iy, iz) OffsetArray of [Bx,By,Bz] 3-vectors, using the grid's own indices (origin r0, spacing dr, anchor = beginning).

_write_field_group(g1, name, field, unit_dim, unit_sym)

Write one field group ("magneticField"/"electricField") from an (ix,iy,iz) OffsetArray of 3-vectors.

_write_fixed_str_array(parent, name, strs::AbstractVector{<:AbstractString})

Write a fixed-length (null-terminated, ASCII) string-array attribute, matching Bmad's hdf5_write_attribute_string rank-1. HDF5.jl writes String arrays as variable-length strings by default, which Bmad's reader cannot convert into its fixed character buffers (it aborts on axisLabels).

gg_to_bmad_curves(fit) -> (cs, cc, c0c, nplanes, nd_max, kmax)

Compute the Bmad azimuthal-harmonic GG derivative towers from a loaded gg_calc_fit result (fit, the GGFit struct returned by read_gg_fit). Returns

cs[(m,j)]  :: Vector  -- C^{[j]}_{m,sin}(plane)  (normal multipole m)
cc[(m,j)]  :: Vector  -- C^{[j]}_{m,cos}(plane)  (skew multipole m)
c0c[j]     :: Vector  -- C^{[j]}_{0,cos}(plane)  (solenoid, j ≥ 1)

each a per-plane vector, for j = 0 … nd_max.