Add/general shaing callen formula

[GeorgGrassler]

The generalized Shaing-Callen bootstrap formula formula, where the inner most field line integral was truncated to get rid of the Bootstrap resonances, was derived here for near omnigenous configurations together with an alternative derivation of the Landreman-Catto formula for perfect omnigenous configurations.

Furthermore, the difference between these two formula was estimated analytically. The dominant contribution to the difference in case of omnigenity violation comes from a boundary term. Analytic estimates reveal a sqrt behavior with the violation of omnigenity and this is recaptured by numerical computations in the code.

Similarly, the total Shaing-Callen values are estimated to differ also by a sqrt scaling term, which again matches the analytical estimates.

These numerical experiments were done for a toy field inspired by Landreman-Paul QH of the form

$$ B = 1 + 0.15 \cos (\vartheta + 4\varphi) + \epsilon_1 \cos (\vartheta) ; . $$

While the numerics agree with the analytical estimates, the overall trend of sqrt(violation) is against the expectation that for good configurations like Landreman-Paul QH the truncated general formula and Landreman-Catto should well agree and for bigger violations, the general formula should be prefered over Landreman-Catto. For example, for the investigated Landreman-Paul QH configuration at half radius $s_\mathrm{tor}=0.25$ one finds that

$$ B \approx B_{00} \left(1 + 10^{-1} \cos (\vartheta + 4\varphi) + 3\cdot 10^{-4} \cos (\vartheta) \right); , $$

which results in an estimate for the relative difference of general vs. Landreman-Catto of

$$ (1 - 2 \cdot 10^{-1}) \sqrt{\frac{3\cdot 10^{-4}}{10^{-1} }} \approx 0.04 $$

which does agree with the actual found difference of $\sim 6$%. The problem is that

• (a) this difference is already visible in graphs and therefore significant and
• (b) the truncated formula does match worse than Landreman-Catto to full kinetic simulations like NEO-2

NEO-2 results (blue) and asymptotic model(orange) top: from generalized Shaing-Callen formula, bottom: from Landreman-Catto formula

While it seems the implementation is computationally sound, it is not yet clear how to interpret this result and if this generalized formula fits into the in rabe implemented asymptotic model. $\rightarrow$ Draft

[qodo-code-review]

You are nearing your monthly Qodo Merge usage quota. For more information, please visit here.

PR Compliance Guide 🔍

Below is a summary of compliance checks for this PR:

Security Compliance
🟢	No security concerns identified No security vulnerabilities detected by AI analysis. Human verification advised for critical code.
Ticket Compliance
⚪	🎫 No ticket provided Create ticket/issue
Codebase Duplication Compliance
⚪	Codebase context is not defined Follow the guide to enable codebase context checks.
Custom Compliance
🟢	Generic: Comprehensive Audit Trails Objective: To create a detailed and reliable record of critical system actions for security analysis and compliance. Status: Passed Learn more about managing compliance generic rules or creating your own custom rules
	Generic: Secure Error Handling Objective: To prevent the leakage of sensitive system information through error messages while providing sufficient detail for internal debugging. Status: Passed Learn more about managing compliance generic rules or creating your own custom rules
	Generic: Secure Logging Practices Objective: To ensure logs are useful for debugging and auditing without exposing sensitive information like PII, PHI, or cardholder data. Status: Passed Learn more about managing compliance generic rules or creating your own custom rules
🔴	Generic: Meaningful Naming and Self-Documenting Code Objective: Ensure all identifiers clearly express their purpose and intent, making code self-documenting Status: Misspelled output label: The new diagnostic string guess of non-omnigneous remainder contains a misspelling that reduces readability and searchability of logs/output. Referred Code print *, "guess of non-omnigneous remainder: ", remainder trapped_fraction_prime = calc_trapped_fraction_prime(field, fieldlines, n_eta) Learn more about managing compliance generic rules or creating your own custom rules
	Generic: Robust Error Handling and Edge Case Management Objective: Ensure comprehensive error handling that provides meaningful context and graceful degradation Status: Missing size validation: cumintegrate_over_phi_grid assumes size(phi_grid) >= 2 (uses t(2)), but does not validate this precondition, which can cause out-of-bounds access for small grids. Referred Code function cumintegrate_over_phi_grid(phi_grid, func) result(cumintegral) real(dp), dimension(:), intent(in) :: phi_grid procedure(func_i) :: func real(dp), dimension(size(phi_grid)) :: cumintegral integer :: n_phi real(dp) :: integrand_start, integrand_end real(dp) :: delta_phi, dt real(dp), dimension(size(phi_grid)) :: dphi_dt, t integer :: this n_phi = size(phi_grid) call linspace(0.0_dp, pi, n_phi, t) delta_phi = 0.5_dp*(phi_grid(1) - phi_grid(n_phi)) dphi_dt = -delta_phi*sin(t) dt = t(2) - t(1) cumintegral(1) = 0.0_dp ! as integral limits are same for first element do this = 2, n_phi integrand_start = func(phi_grid(this - 1))*dphi_dt(this - 1) integrand_end = func(phi_grid(this))*dphi_dt(this) ... (clipped 27 lines) Learn more about managing compliance generic rules or creating your own custom rules
⚪	Generic: Security-First Input Validation and Data Handling Objective: Ensure all data inputs are validated, sanitized, and handled securely to prevent vulnerabilities Status: No grid validation: The new integration routines do not validate basic assumptions about phi_grid (e.g., length, ordering/consistency with the t mapping), which may be acceptable for internal scientific usage but cannot be verified as safe from the diff alone. Referred Code function integrate_over_phi_grid(phi_grid, integrand) result(integral) real(dp), dimension(:), intent(in) :: phi_grid real(dp), dimension(:), intent(in) :: integrand real(dp) :: integral real(dp) :: delta_phi integer :: n_phi real(dp), dimension(size(phi_grid)) :: dphi_dt, t real(dp), dimension(size(phi_grid)) :: integrand_in_t n_phi = size(phi_grid) call linspace(0.0_dp, pi, n_phi, t) delta_phi = 0.5_dp*(phi_grid(1) - phi_grid(n_phi)) dphi_dt = -delta_phi*sin(t) integrand_in_t = integrand*dphi_dt if (ieee_is_nan(integrand(1))) then print *, "Warning in cumintegrate_over_phi_grid:" print *, "integrand is NaN at lower limit!" print *, "set to zero as weight dphi_dt=0 there!" integrand_in_t(1) = 0.0_dp end if ... (clipped 64 lines) Learn more about managing compliance generic rules or creating your own custom rules
Update

Compliance status legend

🟢 - Fully Compliant
🟡 - Partial Compliant
🔴 - Not Compliant
⚪ - Requires Further Human Verification
🏷️ - Compliance label

[qodo-code-review]

You are nearing your monthly Qodo Merge usage quota. For more information, please visit here.

PR Code Suggestions ✨

Explore these optional code suggestions:

Category	Suggestion	Impact
Possible issue	Fix potential runtime array shape mismatch To prevent a runtime error, move the allocation and deallocation of the antider_dBdtheta_over_B_cubed array inside the loop to ensure its size matches the phi_grid recalculated in each iteration. src/shaing_callen_remainder.f90 [141-152] res = 0.0_dp -phi_grid = get_phi_integration_grid(fieldlines(1)) -allocate (antider_dBdtheta_over_B_cubed(size(phi_grid))) do this = 1, n_fieldlines this_fieldline = fieldlines(this) phi_grid = get_phi_integration_grid(this_fieldline) + allocate (antider_dBdtheta_over_B_cubed(size(phi_grid))) antider_dBdtheta_over_B_cubed = cumintegrate_over_phi_grid(phi_grid, & wrapper_dBdtheta_over_B_cubed) res = res + integrate_over_phi_grid(phi_grid, antider_dBdtheta_over_B_cubed) + deallocate (antider_dBdtheta_over_B_cubed) + deallocate (phi_grid) this_fieldline = null_fieldline end do -deallocate (antider_dBdtheta_over_B_cubed) Apply / Chat Suggestion importance[1-10]: 9 __ Why: The suggestion correctly identifies a potential runtime error due to an array shape mismatch and provides a correct fix, preventing the program from crashing if fieldlines produce different grid sizes.	High
	Fix inconsistent array indexing in tests Fix a copy-paste error in the test by accessing M_pol, N_tor, and iota from the same fieldline element, qs_fieldlines(1), instead of from three different elements. test/helpers/test_shaing_callen_remainder.f90 [84-86] M_pol = qs_fieldlines(1)%M_pol -N_tor = qs_fieldlines(2)%N_tor -iota = qs_fieldlines(3)%iota +N_tor = qs_fieldlines(1)%N_tor +iota = qs_fieldlines(1)%iota Apply / Chat Suggestion importance[1-10]: 8 __ Why: This suggestion correctly identifies a clear copy-paste error in a test case where different indices are used to access properties that should belong to the same entity, which would lead to incorrect test validation.	Medium
High-level	Refactor global state from wrapper modules The shaing_callen_wrappers module uses module-level variables for state, creating hidden dependencies. This should be refactored to pass context explicitly, for instance by using a derived type to bundle the function pointer with its required data. Examples: src/shaing_callen_remainder.f90 [94-116] this_fieldline = fieldline phi_grid = get_phi_integration_grid(fieldline) n_phi = size(phi_grid) allocate (antider_dBdtheta_over_lambda_cubed(n_phi)) allocate (phi_integrand_F(n_phi)) do this = 1, size(eta_grid) this_eta = eta_grid(this) ... (clipped 13 lines) src/shaing_callen_wrappers.f90 [30-44] function wrapper_dBdtheta_over_lambda_cubed(phi) use fieldline_integrands, only: calc_lambda_squared real(dp), intent(in) :: phi real(dp) :: wrapper_dBdtheta_over_lambda_cubed real(dp) :: theta, B, dB_dx(3), dB_dtheta, lambda theta = this_fieldline%get_theta(phi) call this_field%compute_B_and_dB_dx(theta, phi, B, dB_dx) ... (clipped 5 lines) Solution Walkthrough: Before: ! In shaing_callen_remainder.f90 function calc_F_prime_for_fieldline(fieldline, eta_grid) ... this_fieldline = fieldline ! Set global state do this = 1, size(eta_grid) this_eta = eta_grid(this) ! Set global state antider_dBdtheta_over_lambda_cubed = cumintegrate_over_phi_grid(phi_grid, & wrapper_dBdtheta_over_lambda_cubed) ... this_eta = null_eta ! Unset global state end do this_fieldline = null_fieldline ! Unset global state ... end function ! In shaing_callen_wrappers.f90 function wrapper_dBdtheta_over_lambda_cubed(phi) ! Uses module variables: this_fieldline, this_field, this_eta theta = this_fieldline%get_theta(phi) call this_field%compute_B_and_dB_dx(...) lambda = sqrt(calc_lambda_squared(B, this_eta)) ... end function After: ! In shaing_callen_wrappers.f90 type :: integrand_context_t class(field_t), pointer :: field type(fieldline_t), pointer :: fieldline real(dp) :: eta end type ! In shaing_callen_remainder.f90 function calc_F_prime_for_fieldline(fieldline, eta_grid) type(integrand_context_t) :: context context%field => this_field context%fieldline => fieldline ... do this = 1, size(eta_grid) context%eta = eta_grid(this) antider_dBdtheta_over_lambda_cubed = cumintegrate_over_phi_grid(phi_grid, & wrapper_dBdtheta_over_lambda_cubed, & context) ... end do ... end function ! In shaing_callen_integration.f90 function cumintegrate_over_phi_grid(phi_grid, func, context) ! context is passed through to func end function Suggestion importance[1-10]: 8 __ Why: The suggestion correctly identifies a significant architectural flaw (use of module-level global state) that is heavily used and expanded by the new code, making it error-prone and not thread-safe.	Medium
General	Remove dead code from test function Remove the unused and overwritten assignment to res in the trial_func test function to eliminate dead code. test/unit/test_cumintegrate_over_phi_grid.f90 [63-69] function trial_func(phi) result(res) real(dp), intent(in) :: phi real(dp) :: res - res = sin(phi)/sqrt(1.0_dp - cos(phi)) res = sqrt(2.0_dp)*cos(0.5_dp*phi) end function trial_func Apply / Chat Suggestion importance[1-10]: 4 __ Why: The suggestion correctly identifies a line of dead code in a test function and recommends its removal, which improves code clarity and maintainability.	Low
Update

[krystophny]

Merged origin/main (boozer_sub/libneo restructure, flock refactor, granular targets). Adaptation decisions for your review:

• shaing_callen_remainder moved to src/shaing_callen/ and ported to the flock_of_fieldlines_t API; signatures follow main's (flock, field, n_eta) convention.
• Main had folded the covariant/dr_dAtheta scaling into get_non_omnigenous_remainder; kept that on top of your pitch+magnetic split, so callers pass dr_dAtheta as for calc_lambda_LC. The unscaled parts remain available as get_non_omnigenous_remainder_pitch/_magnetic.
• Added calc_total_lambda_SC(flock, field, n_eta, dr_dAtheta) wrapping your truncated formula; app/main.f90 writes total_lambda_SC_bB to rabe.nc and rabe.dat and NaN-fills it in unsafe mode.
• Kept cumintegrate_over_phi_grid; dropped cumint/defined_integral_i along with main's removal (their tests are gone on main).
• test/external/plot_deviation_landreman_paul_helical.f90 and plot_fieldlines_quasdex.f90 were deleted on main (they need the internal neo_field); accepted the deletion. Your plot_chi_max_over_xi_0 helper survives with explicit (fieldlines, M_pol, N_tor, iota_p) arguments since fieldline_t lost those members.
• Remainder tests live in test_shaing_callen_remainder as on your branch; for the tests main had also evolved I took main's bodies (M_pol/nfp analytic, get_global_B_max).
• test_get_non_omnigenous_remainders tolerances widened to main's fieldline pipeline (analytic label placement, machine-precision maxima): magnetic 1e-16 -> 1e-15 (measured 1.9e-16), pitch 1e-14 -> 1e-11 (measured 1.6e-12). Your scaling test passes unchanged.
• Golden record updated: your split evaluation removes the O(1/n_eta^2) cancellation bias of the monolithic formula, shifting remainder by up to 0.5% on the QH golden case; total_lambda_SC_bB added. Note the golden run now takes ~11 min locally (was much faster) because calc_trapped_fraction_prime runs per surface.
• Python bindings: shaing_callen_remainder added to the f90wrap sources so get_non_omnigenous_remainder stays on FlockOfFieldlines; calc_total_lambda_sc exported and covered by check_symbols.py.

Verification: Debug build clean; ctest quick 33/33 and slow scaling pass; golden run+compare passes against the updated record; f90wrap bindings build and check_symbols.py passes.

[GeorgGrassler]

@krystophny also leave this branch untouched. This is a research branch of sth that I was not sure about and left open for pure documentation. Maybe I should just archive this branch so that it is untouched