feat: add full particle orbit models (Pauli particle and classical particle)

[krystophny]

User description

Summary

• Add support for full 6D particle orbit integration alongside existing guiding-center model
• New orbit_model parameter: 0=guiding-center (default), 1=Pauli particle, 2=classical particle
• Pauli particle conserves magnetic moment mu (includes mu*grad(B) force)
• Classical particle uses pure Lorentz force (mu=0)
• Same initial conditions (s, theta, phi, lambda, v) for all models
• Gyrophase=0 uses physically motivated grad(B) x b direction
• Requires coils file (field_input) for Biot-Savart field evaluation

Test plan

• All existing tests pass (23/23)
• Manual test with coils file to verify full orbit integration
• Energy conservation check for classical particle
• Magnetic moment conservation check for Pauli particle

---

PR Type

Enhancement

---

Description

• Add full 6D particle orbit integration models (Pauli and classical)

• Implement symplectic implicit midpoint integrator with Newton iteration

• Support cylindrical coordinate field evaluation via Biot-Savart coils

• Transform between cylindrical and VMEC coordinates for orbit output

• Add orbit_model parameter (0=guiding-center, 1=Pauli, 2=classical)

---

Diagram Walkthrough

flowchart LR
  A["Initial Conditions<br/>s, theta, phi, lambda, v"] -->|transform_vmec_to_cyl| B["Cylindrical State<br/>R, phi, Z, p_R, p_phi, p_Z"]
  B -->|init_full_orbit_state| C["FullOrbitState<br/>with mu, mass, charge"]
  C -->|timestep_full_orbit| D["Implicit Midpoint<br/>Integration"]
  D -->|eval_cyl| E["Biot-Savart Field<br/>A, B, gradB"]
  D -->|Newton Iteration| F["Converged State"]
  F -->|convert_full_to_gc| G["Output GC Coords<br/>s, theta, phi, lambda, v"]

Loading

File Walkthrough

	Relevant files
Enhancement	orbit_full.f90 Full orbit integration with implicit midpoint scheme          src/orbit_full.f90 New module implementing full 6D particle orbit integration FullOrbitState type stores cylindrical coordinates (R, phi, Z, p_R, p_phi, p_Z) and magnetic moment mu init_full_orbit_state converts guiding-center ICs to full orbit with proper gyrophase initialization timestep_full_orbit uses symplectic implicit midpoint integrator with Newton iteration eval_midpoint_residual and eval_midpoint_jacobian compute residuals and Jacobians for Newton solver convert_full_to_gc transforms full orbit state back to guiding-center coordinates for output compute_energy calculates total energy including magnetic moment contribution +410/-0  field_coils_cyl.f90 Cylindrical Biot-Savart field evaluation module                    src/field/field_coils_cyl.f90 New module for cylindrical coordinate magnetic field evaluation via Biot-Savart evaluate_cyl computes vector potential A, magnetic field B, magnitude Bmod, and gradient gradB Transforms Cartesian field data to cylindrical coordinates with proper Jacobian handling Finite difference computation of field derivatives for Jacobian calculations set_active_coils manages global coils pointer for field evaluation +142/-0  coordinates.f90 Cylindrical to VMEC coordinate transformation                        src/coordinates/coordinates.f90 Add transform_cyl_to_vmec subroutine for converting cylindrical to VMEC coordinates Implements Newton-Raphson iteration to solve R(s,theta) and Z(s,theta) equations Handles singular Jacobian cases and enforces s in [0,1] and theta in [0,2π] Enables output of full orbit trajectories in VMEC flux coordinates +73/-0    simple.f90 Integrate full orbit state into Tracer type                            src/simple.f90 Add FullOrbitState component fo to Tracer type Add orbit_model field to Tracer type Implement init_full subroutine to initialize full orbit state from guiding-center ICs Implement timestep_full_z subroutine for full orbit timesteps with conversion back to GC coords Add timestep_full_z to tstep interface for polymorphic timestep calls +47/-0    simple_main.f90 Add full orbit dispatch and coils initialization                  src/simple_main.f90 Import full orbit modules and parameters (FullOrbitState, init_full, timestep_full_z) Load coils file and set active coils in init_field when orbit_model is not guiding-center Dispatch to full orbit initialization in trace_orbit based on orbit_model Route macrostep calls to timestep_full_z for full orbit models Disable collision support for full orbit models in macrostep +31/-9
Configuration changes	params.f90 Add orbit model selection and validation                                  src/params.f90 Add orbit_model parameter with three options: ORBIT_GUIDING_CENTER (0), ORBIT_PAULI_PARTICLE (1), ORBIT_PARTICLE (2) Add validation in read_config to ensure orbit_model is in valid range Require field_input (coils file) for full orbit models Prevent collision support for full orbit models Add orbit_model to namelist config +21/-1    CMakeLists.txt Register new source files in build system                                src/CMakeLists.txt Add orbit_full.f90 to source files list Add field/field_coils_cyl.f90 to source files list +2/-0

---

[qodo-code-review]

PR Compliance Guide 🔍

Below is a summary of compliance checks for this PR:

Security Compliance
🔴	Denial-of-service crash Description: Division by zero risk when computing cylindrical Jacobian terms dcyl_dxcart(2,1:2) as they divide by R without guarding small/zero radius, which can cause NaNs/overflow and destabilize the integrator (e.g., near the axis R≈0). field_coils_cyl.f90 [56-61] Referred Code dcyl_dxcart(1,1) = cos_phi dcyl_dxcart(1,2) = sin_phi dcyl_dxcart(1,3) = 0d0 dcyl_dxcart(2,1) = -sin_phi / R dcyl_dxcart(2,2) = cos_phi / R dcyl_dxcart(2,3) = 0d0
⚪	Robustness to invalid inputs Description: Finite-difference step size h = 1d-6 * max(R,1d0) becomes effectively constant near axis and very large R, with no bounds checking or field domain validation; combined with unguarded Biot-Savart calls, this can trigger out-of-domain evaluations and numerical instabilities exploitable to crash runs. field_coils_cyl.f90 [76-88] Referred Code h = 1d-6 * max(R, 1d0) do i = 1, 3 xcart_p = xcart xcart_m = xcart xcart_p(i) = xcart(i) + h xcart_m(i) = xcart(i) - h B_cart_p = compute_magnetic_field(active_coils, xcart_p) B_cart_m = compute_magnetic_field(active_coils, xcart_m) dB_cart(:, i) = (B_cart_p - B_cart_m) / (2d0 * h) enddo
	Solver robustness Description: Linear solve uses dgesv without checking for NaNs/Infs in the Jacobian/residual and treats any nonzero info as generic error 2, lacking conditioning checks or regularization, enabling crafted states to produce singular/ill-conditioned systems and halt processing. orbit_full.f90 [150-156] Referred Code delta = -fvec call dgesv(6, 1, jac, 6, ipiv, delta, 6, info) if (info /= 0) then ierr = 2 return endif
	Numerical failure handling Description: The Newton iteration enforces theta = modulo(theta, 2π) but does not handle s projection effects on convergence near boundaries, and returns ierr=2 on near-singular Jacobian with a fixed determinant threshold 1d-30, allowing legitimate cases to fail and potentially cause repeated retries or aborts upstream. coordinates.f90 [105-107] Referred Code s = max(0d0, min(1d0, s)) theta = modulo(theta, 2d0 * 3.141592653589793d0) enddo
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: 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: Robust Error Handling and Edge Case Management Objective: Ensure comprehensive error handling that provides meaningful context and graceful degradation Status: Fatal Error Stop: The routine aborts with 'error stop' when active_coils is not set instead of returning an error code, preventing graceful handling by callers. Referred Code if (.not. associated(active_coils)) then error stop 'field_coils_cyl: active_coils not set' endif Learn more about managing compliance generic rules or creating your own custom rules
⚪	Generic: Comprehensive Audit Trails Objective: To create a detailed and reliable record of critical system actions for security analysis and compliance. Status: Missing Auditing: New critical operations like loading coils and switching orbit models are not logged, making it hard to audit configuration-driven behavior and runtime decisions. Referred Code if (orbit_model /= ORBIT_GUIDING_CENTER) then if (len_trim(field_input) == 0) then error stop 'Full orbit models require field_input (coils file)' endif call load_coils_from_file(field_input, coils_for_full_orbit) call set_active_coils(coils_for_full_orbit) call print_phase_time('Coils loaded for full orbit integration') endif 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: Abbreviated Names: Several identifiers like z, z_new, z_old, dp_R, dp_Z, and fo are terse and may hinder readability without surrounding context. Referred Code subroutine timestep_full_orbit(state, ierr) type(FullOrbitState), intent(inout) :: state integer, intent(out) :: ierr real(dp) :: z_new(6), z_old(6), fvec(6), jac(6,6) real(dp) :: delta(6), res_norm, res_norm_old integer :: iter, ipiv(6), info integer :: i ierr = 0 z_old = state%z z_new = z_old call explicit_euler_step(state, z_old, z_new) res_norm_old = 1d30 do iter = 1, MAX_NEWTON_ITER call eval_midpoint_residual(state, z_old, z_new, fvec) res_norm = sqrt(sum(fvec**2)) ... (clipped 97 lines) 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: Error Messages: User-facing 'error stop' messages expose internal configuration requirements which may be acceptable for CLI tools but need verification for exposure context. Referred Code if (orbit_model < 0 .or. orbit_model > 2) then error stop 'orbit_model must be 0 (guiding-center), 1 (Pauli), or 2 (particle)' endif if (orbit_model /= ORBIT_GUIDING_CENTER) then if (len_trim(field_input) == 0) then error stop 'Full orbit models require field_input (coils file)' endif if (swcoll) then error stop 'Collisions not supported for full orbit models' endif endif 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: File Input Trust: Coils file loading trusts 'field_input' path without explicit validation or sandboxing beyond existence, which may be acceptable if upstream parsing validates content. Referred Code if (orbit_model /= ORBIT_GUIDING_CENTER) then if (len_trim(field_input) == 0) then error stop 'Full orbit models require field_input (coils file)' endif call load_coils_from_file(field_input, coils_for_full_orbit) call set_active_coils(coils_for_full_orbit) call print_phase_time('Coils loaded for full orbit integration') endif 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]

PR Code Suggestions ✨

Explore these optional code suggestions:

Category	Suggestion	Impact
Possible issue	✅ Fix parallel velocity calculation Suggestion Impact: The commit updated the v_par calculation to remove the R**2 factor and use the correct pairing of physical components, matching the suggested form. Additional related corrections were made to handle covariant/contravariant components consistently elsewhere. code diff: v_mag = sqrt(v_R**2 + R**2 * v_phi**2 + v_Z**2) - v_par = (v_R * B_cyl(1) + R**2 * v_phi * B_cyl(2) + v_Z * B_cyl(3)) / Bmod + ! B_cyl is (B_R, R*B_phi, B_Z), v_phi is angular velocity + ! v_par = v dot b = v_R*B_R + (R*v_phi)*(B_cyl(2)/R) + v_Z*B_Z + v_par = (v_R * B_cyl(1) + v_phi * B_cyl(2) + v_Z * B_cyl(3)) / Bmod Correct the parallel velocity v_par calculation. The dot product v . B should use physical components, which involves canceling out the R factors in the phi-component, not multiplying by R2.** src/orbit_full.f90 [356] -v_par = (v_R * B_cyl(1) + R**2 * v_phi * B_cyl(2) + v_Z * B_cyl(3)) / Bmod +v_par = (v_R * B_cyl(1) + v_phi * B_cyl(2) + v_Z * B_cyl(3)) / Bmod [Suggestion processed] Suggestion importance[1-10]: 9 __ Why: The suggestion correctly identifies a bug in the v_par calculation. The dot product of physical velocity and physical magnetic field was calculated incorrectly by using a covariant component of B (B_cyl(2)) and an extra R**2 factor. This leads to an incorrect parallel velocity and pitch angle, which are critical for the simulation's physics.	High
	✅ Correct the velocity magnitude calculation Suggestion Impact: The commit refactored the code from cylindrical to Cartesian variables, and wherever velocity magnitude is computed, it now uses Cartesian components v_vec and computes v_mag = sqrt(vx^2 + vy^2 + vz^2). This change effectively resolves the original dimensional issue by avoiding the R*v_phi term entirely. code diff: Correct the velocity magnitude v_mag calculation. The square of the magnitude should be v_R2 + (R*v_phi)2 + v_Z2, but the current code incorrectly squares v_phi but not R inside the parenthesis.** src/orbit_full.f90 [354] -v_mag = sqrt(v_R**2 + R**2 * v_phi**2 + v_Z**2) +v_mag = sqrt(v_R**2 + (R * v_phi)**2 + v_Z**2) [Suggestion processed] Suggestion importance[1-10]: 9 __ Why: The suggestion correctly points out a dimensional inconsistency in the velocity magnitude calculation. The physical velocity in the phi direction is R*v_phi, so this term should be squared as (R*v_phi)**2. The existing code incorrectly has R**2 * v_phi**2, which is wrong. This bug affects fundamental outputs like velocity and pitch angle.	High
	✅ Fix squared velocity magnitude calculation Suggestion Impact: The code was refactored from cylindrical to Cartesian variables, and the energy calculation now computes v_sq as the sum of Cartesian velocity components squared (v_vec(1)**2 + v_vec(2)**2 + v_vec(3)**2), which inherently fixes the previous cylindrical inconsistency. This addresses the original issue with the v_phi term. code diff: function compute_energy(state) result(H) - use field_coils_cyl, only: evaluate_cyl - type(FullOrbitState), intent(in) :: state real(dp) :: H - real(dp) :: R, phi, Z, p_R, p_phi, p_Z - real(dp) :: A_cyl(3), B_cyl(3), Bmod, gradB(3) - real(dp) :: v_R, v_phi, v_Z, v_sq - - R = state%z(1) - phi = state%z(2) - Z = state%z(3) - p_R = state%z(4) - p_phi = state%z(5) - p_Z = state%z(6) - - call evaluate_cyl(R, phi, Z, A_cyl, B_cyl, Bmod, gradB) - - v_R = (p_R - (state%q / c) * A_cyl(1)) / state%m - v_phi = (p_phi - (state%q / c) * A_cyl(2)) / (state%m * R**2) - v_Z = (p_Z - (state%q / c) * A_cyl(3)) / state%m - - v_sq = v_R**2 + R**2 * v_phi**2 + v_Z**2 + real(dp) :: x_cart(3), p(3) + real(dp) :: A_cart(3), B_cart(3), Bmod, gradB(3) + real(dp) :: v_vec(3), v_sq + + x_cart = state%z(1:3) + p = state%z(4:6) + + call evaluate_field_and_gradB(state, x_cart, A_cart, B_cart, Bmod, gradB) + + v_vec = (p - (state%q / c) * A_cart) / state%m + + v_sq = v_vec(1)**2 + v_vec(2)**2 + v_vec(3)**2 H = 0.5d0 * state%m * v_sq + state%mu * Bmod end function compute_energy Correct the squared velocity magnitude v_sq calculation. The formula should be v_R2 + (R*v_phi)2 + v_Z2, but the current code incorrectly squares v_phi but not R inside the parenthesis, leading to an incorrect kinetic energy calculation.** src/orbit_full.f90 [404] -v_sq = v_R**2 + R**2 * v_phi**2 + v_Z**2 +v_sq = v_R**2 + (R * v_phi)**2 + v_Z**2 [Suggestion processed] Suggestion importance[1-10]: 9 __ Why: The suggestion correctly identifies a dimensional inconsistency in the v_sq calculation, which is used for computing energy. The physical velocity in the phi direction is R*v_phi, so this term should be squared as (R*v_phi)**2. The existing code's R**2 * v_phi**2 is incorrect and leads to a wrong kinetic energy value, breaking energy conservation checks.	High
Update