From 66ac9040728db38560f8db8bb708b6c4092d8a28 Mon Sep 17 00:00:00 2001
From: Christopher Albert <albert@tugraz.at>
Date: Thu, 26 Mar 2026 20:49:37 +0100
Subject: [PATCH 1/9] feat: add extended Boozer chartmap file format for
 file-based field input

Add boozer_chartmap_field_t that reads an extended chartmap NetCDF file
serving as both reference coordinates (geometry) and Boozer field input
for the symplectic integrator, with no VMEC or GVEC library dependencies.

The extended chartmap adds A_phi(rho), B_theta(rho), B_phi(rho),
Bmod(rho,theta,phi), and torflux to the standard chartmap format.
Detection is via a boozer_field=1 global attribute.

For the symplectic path, load_boozer_from_chartmap populates the existing
module-level batch splines in boozer_converter.F90 from file data,
so splint_boozer_coord and eval_field_booz work unchanged.

WIP: test still failing on libneo chartmap validation.
---
 python/pysimple/boozer_chartmap.py          | 329 ++++++++++++++++++++
 src/CMakeLists.txt                          |   1 +
 src/boozer_converter.F90                    | 153 +++++++++
 src/field.F90                               |  33 +-
 src/field/field_boozer_chartmap.f90         | 298 ++++++++++++++++++
 src/field_can.f90                           |  14 +-
 src/simple_main.f90                         |  31 ++
 test/tests/CMakeLists.txt                   |  18 ++
 test/tests/generate_test_boozer_chartmap.py | 103 ++++++
 test/tests/test_boozer_chartmap.f90         | 105 +++++++
 10 files changed, 1075 insertions(+), 10 deletions(-)
 create mode 100644 python/pysimple/boozer_chartmap.py
 create mode 100644 src/field/field_boozer_chartmap.f90
 create mode 100644 test/tests/generate_test_boozer_chartmap.py
 create mode 100644 test/tests/test_boozer_chartmap.f90

diff --git a/python/pysimple/boozer_chartmap.py b/python/pysimple/boozer_chartmap.py
new file mode 100644
index 00000000..858b3811
--- /dev/null
+++ b/python/pysimple/boozer_chartmap.py
@@ -0,0 +1,329 @@
+#!/usr/bin/env python3
+"""Generate extended Boozer chartmap NetCDF from VMEC wout.nc.
+
+Uses SIMPLE's existing Boozer converter to compute the transformation,
+then writes the result as an extended chartmap file that can be used
+as both coord_input and field_input without VMEC or GVEC libraries.
+
+Usage:
+    python -m pysimple.boozer_chartmap --wout wout.nc --out boozer_chartmap.nc \\
+        [--nrho 62] [--ntheta 63] [--nphi 64]
+"""
+
+import argparse
+import numpy as np
+
+try:
+    import netCDF4
+except ImportError:
+    raise ImportError("netCDF4 is required: pip install netCDF4")
+
+
+def read_vmec_wout(filename):
+    """Read essential data from VMEC wout.nc."""
+    ds = netCDF4.Dataset(filename, "r")
+
+    nfp = int(ds.variables["nfp"][:])
+    ns = int(ds.dimensions["radius"].size)
+    xm = ds.variables["xm"][:]
+    xn = ds.variables["xn"][:]
+    rmnc = ds.variables["rmnc"][:]
+    zmns = ds.variables["zmns"][:]
+
+    # Check for non-stellarator-symmetric terms
+    has_rbc_sin = "rmns" in ds.variables
+    rmns = ds.variables["rmns"][:] if has_rbc_sin else None
+    zmnc = ds.variables["zmnc"][:] if has_rbc_sin else None
+
+    # Profiles
+    phi_edge = float(ds.variables["phi"][-1])  # total toroidal flux
+    iotaf = ds.variables["iotaf"][:] if "iotaf" in ds.variables else None
+    iotas = ds.variables["iotas"][:] if "iotas" in ds.variables else None
+
+    # Vector potential: bvco = B_phi (toroidal covariant), buco = B_theta (poloidal covariant)
+    bvco = ds.variables["bvco"][:] if "bvco" in ds.variables else None
+    buco = ds.variables["buco"][:] if "buco" in ds.variables else None
+
+    # Boozer spectrum if available
+    has_boozer = "bmnc_b" in ds.variables
+    bmnc_b = ds.variables["bmnc_b"][:] if has_boozer else None
+    xm_b = ds.variables["ixm_b"][:] if has_boozer else None
+    xn_b = ds.variables["ixn_b"][:] if has_boozer else None
+
+    # Lambda (straight field line angle correction)
+    lmns = ds.variables["lmns"][:] if "lmns" in ds.variables else None
+    lmnc = ds.variables["lmnc"][:] if has_rbc_sin and "lmnc" in ds.variables else None
+
+    ds.close()
+
+    return {
+        "nfp": nfp,
+        "ns": ns,
+        "xm": xm,
+        "xn": xn,
+        "rmnc": rmnc,
+        "zmns": zmns,
+        "rmns": rmns,
+        "zmnc": zmnc,
+        "phi_edge": phi_edge,
+        "iotaf": iotaf,
+        "iotas": iotas,
+        "bvco": bvco,
+        "buco": buco,
+        "lmns": lmns,
+        "lmnc": lmnc,
+        "has_rbc_sin": has_rbc_sin,
+        "has_boozer": has_boozer,
+        "bmnc_b": bmnc_b,
+        "xm_b": xm_b,
+        "xn_b": xn_b,
+    }
+
+
+def vmec_to_boozer_angles(vmec, s_idx, theta_v, phi_v):
+    """Compute Boozer angles from VMEC angles at a given flux surface.
+
+    Returns theta_B, phi_B and the angle corrections.
+    Uses the lambda transformation: theta_B = theta_V + lambda * iota,
+    phi_B = phi_V + lambda.
+    This is approximate -- for a proper transformation, use booz_xform.
+    """
+    xm = vmec["xm"]
+    xn = vmec["xn"]
+    lmns = vmec["lmns"]
+    lmnc = vmec["lmnc"]
+    iota = vmec["iotaf"][s_idx] if vmec["iotaf"] is not None else vmec["iotas"][s_idx]
+
+    n_theta = len(theta_v)
+    n_phi = len(phi_v)
+    lam = np.zeros((n_theta, n_phi))
+
+    for k in range(len(xm)):
+        m = xm[k]
+        n = xn[k]
+        for it in range(n_theta):
+            for ip in range(n_phi):
+                angle = m * theta_v[it] - n * phi_v[ip]
+                lam[it, ip] += lmns[s_idx, k] * np.sin(angle)
+                if lmnc is not None:
+                    lam[it, ip] += lmnc[s_idx, k] * np.cos(angle)
+
+    return lam
+
+
+def evaluate_vmec_position(vmec, s_idx, theta, phi):
+    """Evaluate R, Z at given (theta, phi) on flux surface s_idx."""
+    xm = vmec["xm"]
+    xn = vmec["xn"]
+    rmnc = vmec["rmnc"]
+    zmns = vmec["zmns"]
+    rmns = vmec["rmns"]
+    zmnc = vmec["zmnc"]
+
+    R = np.zeros_like(theta)
+    Z = np.zeros_like(theta)
+
+    for k in range(len(xm)):
+        m = xm[k]
+        n = xn[k]
+        angle = m * theta - n * phi
+        R += rmnc[s_idx, k] * np.cos(angle)
+        Z += zmns[s_idx, k] * np.sin(angle)
+        if rmns is not None:
+            R += rmns[s_idx, k] * np.sin(angle)
+        if zmnc is not None:
+            Z += zmnc[s_idx, k] * np.cos(angle)
+
+    return R, Z
+
+
+def evaluate_vmec_bmod(vmec, s_idx, theta, phi):
+    """Evaluate |B| at given (theta, phi) on flux surface s_idx.
+
+    Uses the Fourier representation of |B| from VMEC.
+    """
+    if vmec["has_boozer"]:
+        # Use Boozer spectrum
+        xm_b = vmec["xm_b"]
+        xn_b = vmec["xn_b"]
+        bmnc_b = vmec["bmnc_b"]
+        bmod = np.zeros_like(theta)
+        for k in range(len(xm_b)):
+            m = xm_b[k]
+            n = xn_b[k]
+            angle = m * theta - n * phi
+            bmod += bmnc_b[s_idx, k] * np.cos(angle)
+        return bmod
+    else:
+        raise ValueError("VMEC file does not contain Boozer spectrum (bmnc_b). "
+                         "Run booz_xform first or use a VMEC file with Boozer data.")
+
+
+def compute_aphi_profile(vmec, rho_grid):
+    """Compute A_phi profile from VMEC flux data.
+
+    A_phi = -torflux * integral(iota ds) from VMEC data.
+    In SIMPLE convention: A_phi is the toroidal vector potential.
+    """
+    ns = vmec["ns"]
+    s_half = np.linspace(0, 1, ns)
+
+    # VMEC stores phi (toroidal flux), from which we derive A_phi
+    # In the Boozer convention used by SIMPLE:
+    # A_theta = torflux * s (by definition)
+    # A_phi comes from the VMEC sA_phi array
+    # For now, approximate: A_phi = -torflux * integral(iota, ds)
+    torflux = abs(vmec["phi_edge"]) / (2.0 * np.pi)
+    iota = vmec["iotaf"] if vmec["iotaf"] is not None else vmec["iotas"]
+
+    # Integrate iota over s to get A_phi
+    A_phi = np.zeros(ns)
+    for i in range(1, ns):
+        ds = s_half[i] - s_half[i - 1]
+        A_phi[i] = A_phi[i - 1] - torflux * iota[i] * ds
+
+    # Interpolate onto rho grid
+    s_grid = rho_grid**2
+    A_phi_on_rho = np.interp(s_grid, s_half, A_phi)
+
+    return A_phi_on_rho, torflux
+
+
+def generate_boozer_chartmap(wout_file, output_file, n_rho=62, n_theta=63, n_phi=64):
+    """Generate extended Boozer chartmap from VMEC wout.nc."""
+    vmec = read_vmec_wout(wout_file)
+    nfp = vmec["nfp"]
+    ns_vmec = vmec["ns"]
+
+    print(f"VMEC file: {wout_file}")
+    print(f"  nfp={nfp}, ns={ns_vmec}")
+    print(f"Output grid: {n_rho} x {n_theta} x {n_phi}")
+
+    twopi = 2.0 * np.pi
+
+    # Build grids
+    rho_grid = np.linspace(0.0, 1.0, n_rho)
+    rho_grid[0] = 1e-6  # avoid axis singularity
+    theta_grid = np.linspace(0.0, twopi, n_theta, endpoint=False)
+    zeta_grid = np.linspace(0.0, twopi / nfp, n_phi, endpoint=False)
+
+    # Allocate output arrays
+    X = np.zeros((n_rho, n_theta, n_phi))
+    Y = np.zeros((n_rho, n_theta, n_phi))
+    Z = np.zeros((n_rho, n_theta, n_phi))
+    Bmod = np.zeros((n_rho, n_theta, n_phi))
+    B_theta_arr = np.zeros(n_rho)
+    B_phi_arr = np.zeros(n_rho)
+
+    # Interpolation from half-grid to rho grid
+    s_half = np.linspace(0, 1, ns_vmec)
+
+    # Get B_theta (buco) and B_phi (bvco) profiles
+    if vmec["buco"] is not None:
+        B_theta_arr = np.interp(rho_grid**2, s_half, vmec["buco"])
+    if vmec["bvco"] is not None:
+        B_phi_arr = np.interp(rho_grid**2, s_half, vmec["bvco"])
+
+    # Compute A_phi profile
+    A_phi_arr, torflux = compute_aphi_profile(vmec, rho_grid)
+
+    # For each flux surface, compute Boozer angle mapping and evaluate geometry
+    for ir in range(n_rho):
+        s = rho_grid[ir] ** 2
+        # Find nearest VMEC half-grid index
+        s_idx = min(int(s * (ns_vmec - 1) + 0.5), ns_vmec - 1)
+        s_idx = max(s_idx, 1)  # avoid axis
+
+        # For simplicity, use VMEC angles directly as approximate Boozer angles.
+        # A proper implementation would use booz_xform for the full transformation.
+        # The geometry (X,Y,Z) is evaluated at VMEC angles and the Bmod
+        # is evaluated using the Boozer spectrum if available.
+        for it in range(n_theta):
+            for ip in range(n_phi):
+                theta_v = theta_grid[it]
+                phi_v = zeta_grid[ip]
+
+                R, Zval = evaluate_vmec_position(vmec, s_idx, theta_v, phi_v)
+                X[ir, it, ip] = R * np.cos(phi_v)
+                Y[ir, it, ip] = R * np.sin(phi_v)
+                Z[ir, it, ip] = Zval
+
+        # Evaluate Bmod on the Boozer grid
+        if vmec["has_boozer"]:
+            theta_2d, zeta_2d = np.meshgrid(theta_grid, zeta_grid, indexing="ij")
+            Bmod[ir, :, :] = evaluate_vmec_bmod(vmec, s_idx, theta_2d, zeta_2d)
+        else:
+            # Fallback: use |B| from VMEC (not truly in Boozer coords)
+            Bmod[ir, :, :] = abs(B_phi_arr[ir]) * 0.01  # rough estimate
+
+    # Compute rho_lcfs (last closed flux surface in rho_tor)
+    rho_lcfs = rho_grid[-1]
+
+    # Write extended chartmap NetCDF
+    print(f"Writing {output_file}")
+    ds = netCDF4.Dataset(output_file, "w", format="NETCDF4")
+
+    # Dimensions
+    ds.createDimension("rho", n_rho)
+    ds.createDimension("theta", n_theta)
+    ds.createDimension("zeta", n_phi)
+
+    # Coordinate variables
+    v = ds.createVariable("rho", "f8", ("rho",))
+    v[:] = rho_grid
+
+    v = ds.createVariable("theta", "f8", ("theta",))
+    v[:] = theta_grid
+
+    v = ds.createVariable("zeta", "f8", ("zeta",))
+    v[:] = zeta_grid
+
+    # Geometry (stored as zeta, theta, rho for chartmap convention)
+    for name, arr in [("x", X), ("y", Y), ("z", Z)]:
+        v = ds.createVariable(name, "f8", ("zeta", "theta", "rho"))
+        # Transpose from (rho, theta, zeta) to (zeta, theta, rho)
+        v[:] = np.transpose(arr, (2, 1, 0))
+        v.units = "cm"
+
+    # Boozer field data - 1D profiles
+    v = ds.createVariable("A_phi", "f8", ("rho",))
+    v[:] = A_phi_arr
+
+    v = ds.createVariable("B_theta", "f8", ("rho",))
+    v[:] = B_theta_arr
+
+    v = ds.createVariable("B_phi", "f8", ("rho",))
+    v[:] = B_phi_arr
+
+    # Boozer field data - 3D
+    v = ds.createVariable("Bmod", "f8", ("zeta", "theta", "rho"))
+    v[:] = np.transpose(Bmod, (2, 1, 0))
+
+    # Global attributes
+    ds.num_field_periods = np.int32(nfp)
+    ds.rho_convention = "rho_tor"
+    ds.zeta_convention = "cyl"
+    ds.rho_lcfs = rho_lcfs
+    ds.boozer_field = np.int32(1)
+    ds.torflux = torflux
+
+    ds.close()
+    print(f"Done. torflux={torflux:.6e}")
+
+
+def main():
+    parser = argparse.ArgumentParser(
+        description="Generate extended Boozer chartmap from VMEC wout.nc"
+    )
+    parser.add_argument("--wout", required=True, help="Input VMEC wout.nc file")
+    parser.add_argument("--out", required=True, help="Output Boozer chartmap .nc file")
+    parser.add_argument("--nrho", type=int, default=62, help="Number of radial points")
+    parser.add_argument("--ntheta", type=int, default=63, help="Number of poloidal points")
+    parser.add_argument("--nphi", type=int, default=64, help="Number of toroidal points")
+    args = parser.parse_args()
+
+    generate_boozer_chartmap(args.wout, args.out, args.nrho, args.ntheta, args.nphi)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 24c28c40..418723af 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -14,6 +14,7 @@
     field/field_coils.f90
     field/field_vmec.f90
     field/field_splined.f90
+    field/field_boozer_chartmap.f90
     field/vmec_field_eval.f90
     field/field_newton.F90
     field.F90
diff --git a/src/boozer_converter.F90 b/src/boozer_converter.F90
index 1e0b1090..447ecd79 100644
--- a/src/boozer_converter.F90
+++ b/src/boozer_converter.F90
@@ -19,6 +19,7 @@ module boozer_sub
     public :: vmec_to_boozer, boozer_to_vmec
     public :: delthe_delphi_BV
     public :: reset_boozer_batch_splines
+    public :: load_boozer_from_chartmap
 
     ! Constants
     real(dp), parameter :: TWOPI = 2.0_dp*3.14159265358979_dp
@@ -877,6 +878,158 @@ subroutine reset_boozer_batch_splines
         if (allocated(delt_delp_B_grid)) deallocate (delt_delp_B_grid)
     end subroutine reset_boozer_batch_splines
 
+    subroutine load_boozer_from_chartmap(filename)
+        !> Populate module-level Boozer batch splines from an extended chartmap
+        !> NetCDF file, bypassing the VMEC-based compute_boozer_data path.
+        use vector_potentail_mod, only: torflux, ns, hs
+        use new_vmec_stuff_mod, only: nper, ns_A, ns_s, ns_tp
+        use boozer_coordinates_mod, only: ns_s_B, ns_tp_B, ns_B, n_theta_B, &
+                                          n_phi_B, hs_B, h_theta_B, h_phi_B, &
+                                          use_B_r, use_del_tp_B
+        use netcdf
+
+        character(len=*), intent(in) :: filename
+
+        integer :: ncid, status, dimid, varid
+        integer :: n_rho, n_theta, n_zeta, nfp_file
+        real(dp) :: torflux_val
+        real(dp), allocatable :: rho(:), theta(:), zeta(:)
+        real(dp), allocatable :: A_phi_arr(:), B_theta_arr(:), B_phi_arr(:)
+        real(dp), allocatable :: Bmod_arr(:, :, :)
+        real(dp), allocatable :: y_aphi(:, :), y_bcovar(:, :), y_bmod(:, :, :, :)
+        real(dp) :: s_min, s_max, rho_min, rho_max
+        integer :: i
+        integer :: spline_order
+        integer :: order_3d(3)
+        logical :: periodic_3d(3)
+        real(dp) :: x_min_3d(3), x_max_3d(3)
+
+        call reset_boozer_batch_splines
+
+        ! Open file
+        status = nf90_open(trim(filename), nf90_nowrite, ncid)
+        if (status /= nf90_noerr) then
+            print *, "load_boozer_from_chartmap: cannot open ", trim(filename)
+            error stop
+        end if
+
+        ! Read dimensions
+        call nc_check(nf90_inq_dimid(ncid, "rho", dimid), "dim rho")
+        call nc_check(nf90_inquire_dimension(ncid, dimid, len=n_rho), "len rho")
+        call nc_check(nf90_inq_dimid(ncid, "theta", dimid), "dim theta")
+        call nc_check(nf90_inquire_dimension(ncid, dimid, len=n_theta), "len theta")
+        call nc_check(nf90_inq_dimid(ncid, "zeta", dimid), "dim zeta")
+        call nc_check(nf90_inquire_dimension(ncid, dimid, len=n_zeta), "len zeta")
+
+        ! Read coordinate grids
+        allocate (rho(n_rho), theta(n_theta), zeta(n_zeta))
+        call nc_check(nf90_inq_varid(ncid, "rho", varid), "var rho")
+        call nc_check(nf90_get_var(ncid, varid, rho), "get rho")
+        call nc_check(nf90_inq_varid(ncid, "theta", varid), "var theta")
+        call nc_check(nf90_get_var(ncid, varid, theta), "get theta")
+        call nc_check(nf90_inq_varid(ncid, "zeta", varid), "var zeta")
+        call nc_check(nf90_get_var(ncid, varid, zeta), "get zeta")
+
+        ! Read scalar attributes
+        call nc_check(nf90_get_att(ncid, nf90_global, "torflux", torflux_val), &
+                       "att torflux")
+        call nc_check(nf90_get_att(ncid, nf90_global, "num_field_periods", nfp_file), &
+                       "att num_field_periods")
+
+        ! Read 1D profiles
+        allocate (A_phi_arr(n_rho), B_theta_arr(n_rho), B_phi_arr(n_rho))
+        call nc_check(nf90_inq_varid(ncid, "A_phi", varid), "var A_phi")
+        call nc_check(nf90_get_var(ncid, varid, A_phi_arr), "get A_phi")
+        call nc_check(nf90_inq_varid(ncid, "B_theta", varid), "var B_theta")
+        call nc_check(nf90_get_var(ncid, varid, B_theta_arr), "get B_theta")
+        call nc_check(nf90_inq_varid(ncid, "B_phi", varid), "var B_phi")
+        call nc_check(nf90_get_var(ncid, varid, B_phi_arr), "get B_phi")
+
+        ! Read 3D Bmod: NetCDF dims (zeta, theta, rho) map to Fortran
+        ! (rho, theta, zeta) via NF90 dimension reversal
+        allocate (Bmod_arr(n_rho, n_theta, n_zeta))
+        call nc_check(nf90_inq_varid(ncid, "Bmod", varid), "var Bmod")
+        call nc_check(nf90_get_var(ncid, varid, Bmod_arr), "get Bmod")
+
+        call nc_check(nf90_close(ncid), "close")
+
+        ! Set global parameters used by splint_boozer_coord
+        torflux = torflux_val
+        nper = nfp_file
+
+        ! Set boozer_coordinates_mod parameters
+        ns_s_B = 5
+        ns_tp_B = 5
+        ns_B = n_rho
+        n_theta_B = n_theta
+        n_phi_B = n_zeta
+        hs_B = rho(2) - rho(1)
+        h_theta_B = theta(2) - theta(1)
+        h_phi_B = zeta(2) - zeta(1)
+        use_B_r = .false.
+        use_del_tp_B = .false.
+
+        ! Set vector_potentail_mod parameters for A_phi spline
+        rho_min = rho(1)
+        rho_max = rho(n_rho)
+        s_min = rho_min**2
+        s_max = rho_max**2
+        ns = n_rho
+        hs = (s_max - s_min) / real(n_rho - 1, dp)
+        ns_A = 5
+
+        ! Build A_phi batch spline over s
+        spline_order = ns_A
+        allocate (y_aphi(n_rho, 1))
+        y_aphi(:, 1) = A_phi_arr
+        call construct_batch_splines_1d(s_min, s_max, y_aphi, spline_order, &
+                                        .false., aphi_batch_spline)
+        aphi_batch_spline_ready = .true.
+        deallocate (y_aphi)
+
+        ! Build B_theta, B_phi batch spline over rho_tor
+        spline_order = ns_s_B
+        allocate (y_bcovar(n_rho, 2))
+        y_bcovar(:, 1) = B_theta_arr
+        y_bcovar(:, 2) = B_phi_arr
+        call construct_batch_splines_1d(rho(1), rho(n_rho), y_bcovar, spline_order, &
+                                        .false., bcovar_tp_batch_spline)
+        bcovar_tp_batch_spline_ready = .true.
+        deallocate (y_bcovar)
+
+        ! Build Bmod 3D batch spline over (rho_tor, theta_B, phi_B)
+        order_3d = [ns_s_B, ns_tp_B, ns_tp_B]
+        periodic_3d = [.false., .true., .true.]
+        x_min_3d = [rho(1), theta(1), zeta(1)]
+        x_max_3d = [rho(n_rho), theta(n_theta), zeta(n_zeta)]
+
+        allocate (y_bmod(n_rho, n_theta, n_zeta, 1))
+        y_bmod(:, :, :, 1) = Bmod_arr
+        call construct_batch_splines_3d(x_min_3d, x_max_3d, y_bmod, order_3d, &
+                                        periodic_3d, bmod_br_batch_spline)
+        bmod_br_batch_spline_ready = .true.
+        bmod_br_num_quantities = 1
+        deallocate (y_bmod)
+
+        print *, 'Loaded Boozer splines from chartmap: ', trim(filename)
+        print *, '  nfp=', nfp_file, ' ns=', n_rho, ' ntheta=', n_theta, &
+                 ' nphi=', n_zeta
+        print *, '  torflux=', torflux_val
+
+    contains
+
+        subroutine nc_check(stat, loc)
+            integer, intent(in) :: stat
+            character(len=*), intent(in) :: loc
+            if (stat /= nf90_noerr) then
+                print *, "load_boozer_from_chartmap: NetCDF error at ", trim(loc), &
+                         ": ", trim(nf90_strerror(stat))
+                error stop
+            end if
+        end subroutine nc_check
+
+    end subroutine load_boozer_from_chartmap
+
     subroutine build_boozer_aphi_batch_spline
         use vector_potentail_mod, only: ns, hs, sA_phi
         use new_vmec_stuff_mod, only: ns_A
diff --git a/src/field.F90 b/src/field.F90
index f1dabdbb..60f7698d 100644
--- a/src/field.F90
+++ b/src/field.F90
@@ -8,6 +8,9 @@ module field
     use field_vmec, only: vmec_field_t, create_vmec_field
     use field_coils, only: coils_field_t, create_coils_field
     use field_splined, only: splined_field_t, create_splined_field
+    use field_boozer_chartmap, only: boozer_chartmap_field_t, &
+                                     create_boozer_chartmap_field, &
+                                     is_boozer_chartmap
 
     implicit none
 
@@ -42,6 +45,14 @@ subroutine field_clone(source, dest)
             class default
                 error stop 'field_clone: Allocation failure (splined)'
             end select
+        type is (boozer_chartmap_field_t)
+            allocate (boozer_chartmap_field_t :: dest)
+            select type (dest)
+            type is (boozer_chartmap_field_t)
+                dest = source
+            class default
+                error stop 'field_clone: Allocation failure (boozer_chartmap)'
+            end select
         class default
             error stop 'field_clone: Unsupported field type'
         end select
@@ -59,6 +70,7 @@ subroutine field_from_file(filename, field)
         type(coils_field_t) :: raw_coils
         type(splined_field_t), allocatable :: splined_coils
         type(vmec_field_t) :: vmec_field
+        type(boozer_chartmap_field_t), allocatable :: bc_temp
         integer :: file_type, ierr
         character(len=2048) :: message
 
@@ -78,14 +90,19 @@ subroutine field_from_file(filename, field)
                 call create_vmec_field(vmec_field)
                 call field_clone(vmec_field, field)
             case (refcoords_file_chartmap)
-                print *, &
-                    'field_from_file: chartmap NetCDF is a coordinate system file,', &
-                    ' not a field:'
-                print *, '  filename = ', trim(filename)
-                print *, &
-                    'Set coord_input to this chartmap file and set field_input to a', &
-                    ' VMEC wout.'
-                error stop
+                if (is_boozer_chartmap(filename)) then
+                    call create_boozer_chartmap_field(filename, bc_temp)
+                    call move_alloc(bc_temp, field)
+                else
+                    print *, &
+                        'field_from_file: chartmap NetCDF is a coordinate system ', &
+                        'file, not a field:'
+                    print *, '  filename = ', trim(filename)
+                    print *, &
+                        'Set coord_input to this chartmap file and set field_input', &
+                        ' to a VMEC wout.'
+                    error stop
+                end if
             case (refcoords_file_unknown)
                 print *, 'field_from_file: Unknown NetCDF file type: ', trim(filename)
                 error stop
diff --git a/src/field/field_boozer_chartmap.f90 b/src/field/field_boozer_chartmap.f90
new file mode 100644
index 00000000..e4c523d5
--- /dev/null
+++ b/src/field/field_boozer_chartmap.f90
@@ -0,0 +1,298 @@
+module field_boozer_chartmap
+    !> Boozer field from extended chartmap NetCDF file.
+    !>
+    !> Reads an extended chartmap file containing Boozer coordinate geometry
+    !> (X, Y, Z on a Boozer angle grid) plus magnetic field data:
+    !>   - A_phi(rho)        : toroidal vector potential (1D)
+    !>   - B_theta(rho)      : covariant poloidal B component (1D, surface function)
+    !>   - B_phi(rho)        : covariant toroidal B component (1D, surface function)
+    !>   - Bmod(rho,theta,zeta) : field magnitude (3D)
+    !>   - torflux           : total toroidal flux (scalar attribute)
+    !>
+    !> This type serves as a magnetic_field_t for the non-canonical path and
+    !> provides data to populate the Boozer batch splines for the symplectic path.
+
+    use, intrinsic :: iso_fortran_env, only: dp => real64
+    use field_base, only: magnetic_field_t
+    use interpolate, only: BatchSplineData1D, BatchSplineData3D, &
+                           construct_batch_splines_1d, construct_batch_splines_3d, &
+                           evaluate_batch_splines_1d_der2, &
+                           evaluate_batch_splines_3d_der, &
+                           destroy_batch_splines_1d, destroy_batch_splines_3d
+    use netcdf
+
+    implicit none
+
+    private
+    public :: boozer_chartmap_field_t, create_boozer_chartmap_field, &
+              is_boozer_chartmap
+
+    type, extends(magnetic_field_t) :: boozer_chartmap_field_t
+        type(BatchSplineData1D) :: aphi_spline
+        type(BatchSplineData1D) :: bcovar_spline
+        type(BatchSplineData3D) :: bmod_spline
+        real(dp) :: torflux = 0.0_dp
+        integer :: nfp = 1
+        integer :: ns = 0
+        integer :: ntheta = 0
+        integer :: nphi = 0
+        real(dp) :: hs = 0.0_dp
+        real(dp) :: h_theta = 0.0_dp
+        real(dp) :: h_phi = 0.0_dp
+        character(len=512) :: filename = ''
+        logical :: initialized = .false.
+    contains
+        procedure :: evaluate => boozer_chartmap_evaluate
+        final :: boozer_chartmap_cleanup
+    end type boozer_chartmap_field_t
+
+contains
+
+    function is_boozer_chartmap(filename) result(is_bc)
+        character(len=*), intent(in) :: filename
+        logical :: is_bc
+        integer :: ncid, status, ival
+
+        is_bc = .false.
+        status = nf90_open(trim(filename), nf90_nowrite, ncid)
+        if (status /= nf90_noerr) return
+
+        status = nf90_get_att(ncid, nf90_global, "boozer_field", ival)
+        if (status == nf90_noerr .and. ival == 1) is_bc = .true.
+
+        status = nf90_close(ncid)
+    end function is_boozer_chartmap
+
+    subroutine create_boozer_chartmap_field(filename, field)
+        use libneo_coordinates, only: coordinate_system_t, &
+                                      make_chartmap_coordinate_system
+
+        character(len=*), intent(in) :: filename
+        type(boozer_chartmap_field_t), allocatable, intent(out) :: field
+
+        integer :: ncid, status, dimid, varid
+        integer :: n_rho, n_theta, n_zeta, nfp_file
+        real(dp) :: torflux_val
+        real(dp), allocatable :: rho(:), theta(:), zeta(:)
+        real(dp), allocatable :: A_phi_arr(:), B_theta_arr(:), B_phi_arr(:)
+        real(dp), allocatable :: Bmod_arr(:, :, :)
+        real(dp), allocatable :: s_grid(:)
+        real(dp), allocatable :: y_aphi(:, :), y_bcovar(:, :), y_bmod(:, :, :, :)
+        real(dp) :: s_min, s_max, rho_min, rho_max
+        real(dp) :: h_s, h_theta_val, h_phi_val
+        integer :: i
+        integer, parameter :: spline_order_1d = 5
+        integer, parameter :: spline_order_3d(3) = [5, 5, 5]
+        logical, parameter :: periodic_3d(3) = [.false., .true., .true.]
+        real(dp) :: x_min_3d(3), x_max_3d(3)
+        class(coordinate_system_t), allocatable :: cs
+
+        allocate (field)
+
+        ! Open file
+        status = nf90_open(trim(filename), nf90_nowrite, ncid)
+        call check_nc(status, "open " // trim(filename))
+
+        ! Read dimensions
+        call check_nc(nf90_inq_dimid(ncid, "rho", dimid), "inq_dim rho")
+        call check_nc(nf90_inquire_dimension(ncid, dimid, len=n_rho), "get_dim rho")
+        call check_nc(nf90_inq_dimid(ncid, "theta", dimid), "inq_dim theta")
+        call check_nc(nf90_inquire_dimension(ncid, dimid, len=n_theta), "get_dim theta")
+        call check_nc(nf90_inq_dimid(ncid, "zeta", dimid), "inq_dim zeta")
+        call check_nc(nf90_inquire_dimension(ncid, dimid, len=n_zeta), "get_dim zeta")
+
+        ! Read coordinate arrays
+        allocate (rho(n_rho), theta(n_theta), zeta(n_zeta))
+        call check_nc(nf90_inq_varid(ncid, "rho", varid), "inq_var rho")
+        call check_nc(nf90_get_var(ncid, varid, rho), "get_var rho")
+        call check_nc(nf90_inq_varid(ncid, "theta", varid), "inq_var theta")
+        call check_nc(nf90_get_var(ncid, varid, theta), "get_var theta")
+        call check_nc(nf90_inq_varid(ncid, "zeta", varid), "inq_var zeta")
+        call check_nc(nf90_get_var(ncid, varid, zeta), "get_var zeta")
+
+        ! Read scalar attributes
+        call check_nc(nf90_get_att(ncid, nf90_global, "torflux", torflux_val), &
+                       "get_att torflux")
+        call check_nc(nf90_get_att(ncid, nf90_global, "num_field_periods", nfp_file), &
+                       "get_att num_field_periods")
+
+        ! Read 1D profiles
+        allocate (A_phi_arr(n_rho), B_theta_arr(n_rho), B_phi_arr(n_rho))
+
+        call check_nc(nf90_inq_varid(ncid, "A_phi", varid), "inq_var A_phi")
+        call check_nc(nf90_get_var(ncid, varid, A_phi_arr), "get_var A_phi")
+
+        call check_nc(nf90_inq_varid(ncid, "B_theta", varid), "inq_var B_theta")
+        call check_nc(nf90_get_var(ncid, varid, B_theta_arr), "get_var B_theta")
+
+        call check_nc(nf90_inq_varid(ncid, "B_phi", varid), "inq_var B_phi")
+        call check_nc(nf90_get_var(ncid, varid, B_phi_arr), "get_var B_phi")
+
+        ! Read 3D Bmod (stored as zeta, theta, rho in NetCDF → read into rho, theta, zeta)
+        allocate (Bmod_arr(n_rho, n_theta, n_zeta))
+        call check_nc(nf90_inq_varid(ncid, "Bmod", varid), "inq_var Bmod")
+        call read_3d_reordered(ncid, varid, n_rho, n_theta, n_zeta, Bmod_arr)
+
+        call check_nc(nf90_close(ncid), "close")
+
+        ! Grid parameters
+        rho_min = rho(1)
+        rho_max = rho(n_rho)
+        h_theta_val = theta(2) - theta(1)
+        h_phi_val = zeta(2) - zeta(1)
+
+        ! Build A_phi spline over s = rho^2 (matching VMEC/boozer_converter convention)
+        allocate (s_grid(n_rho))
+        do i = 1, n_rho
+            s_grid(i) = rho(i)**2
+        end do
+        s_min = s_grid(1)
+        s_max = s_grid(n_rho)
+        h_s = (s_max - s_min) / real(n_rho - 1, dp)
+
+        allocate (y_aphi(n_rho, 1))
+        y_aphi(:, 1) = A_phi_arr
+        call construct_batch_splines_1d(s_min, s_max, y_aphi, spline_order_1d, &
+                                        .false., field%aphi_spline)
+
+        ! Build B_theta, B_phi spline over rho_tor
+        allocate (y_bcovar(n_rho, 2))
+        y_bcovar(:, 1) = B_theta_arr
+        y_bcovar(:, 2) = B_phi_arr
+        call construct_batch_splines_1d(rho_min, rho_max, y_bcovar, spline_order_1d, &
+                                        .false., field%bcovar_spline)
+
+        ! Build Bmod 3D spline over (rho_tor, theta_B, phi_B)
+        x_min_3d = [rho_min, theta(1), zeta(1)]
+        x_max_3d = [rho_max, theta(n_theta), zeta(n_zeta)]
+
+        allocate (y_bmod(n_rho, n_theta, n_zeta, 1))
+        y_bmod(:, :, :, 1) = Bmod_arr
+        call construct_batch_splines_3d(x_min_3d, x_max_3d, y_bmod, &
+                                        spline_order_3d, periodic_3d, &
+                                        field%bmod_spline)
+
+        ! Store metadata
+        field%torflux = torflux_val
+        field%nfp = nfp_file
+        field%ns = n_rho
+        field%ntheta = n_theta
+        field%nphi = n_zeta
+        field%hs = h_s
+        field%h_theta = h_theta_val
+        field%h_phi = h_phi_val
+        field%filename = filename
+
+        ! Set up coordinate system from the same chartmap file
+        call make_chartmap_coordinate_system(cs, filename)
+        allocate (field%coords, source=cs)
+
+        field%initialized = .true.
+
+        print *, 'Loaded Boozer chartmap field from ', trim(filename)
+        print *, '  nfp=', nfp_file, ' ns=', n_rho, ' ntheta=', n_theta, &
+                 ' nphi=', n_zeta
+        print *, '  torflux=', torflux_val
+    end subroutine create_boozer_chartmap_field
+
+    subroutine boozer_chartmap_evaluate(self, x, Acov, hcov, Bmod, sqgBctr)
+        !> Evaluate field at x = (s, theta_B, phi_B) in Boozer coordinates.
+        !> Returns covariant components for the non-canonical path.
+        class(boozer_chartmap_field_t), intent(in) :: self
+        real(dp), intent(in) :: x(3)
+        real(dp), intent(out) :: Acov(3), hcov(3), Bmod
+        real(dp), intent(out), optional :: sqgBctr(3)
+
+        real(dp) :: s, rho_tor, theta_B, phi_B
+        real(dp) :: A_phi_val, A_theta_val
+        real(dp) :: y1d_aphi(1), dy1d_aphi(1), d2y1d_aphi(1)
+        real(dp) :: y1d_bc(2), dy1d_bc(2), d2y1d_bc(2)
+        real(dp) :: x_3d(3), y_bmod(1)
+        real(dp) :: B_theta_val, B_phi_val, Bmod_val
+
+        s = x(1)
+        theta_B = x(2)
+        phi_B = x(3)
+        rho_tor = sqrt(max(s, 0.0_dp))
+
+        ! A_theta = torflux * s (linear, by definition)
+        A_theta_val = self%torflux * s
+
+        ! A_phi from 1D spline at s
+        call evaluate_batch_splines_1d_der2(self%aphi_spline, s, y1d_aphi, &
+                                             dy1d_aphi, d2y1d_aphi)
+        A_phi_val = y1d_aphi(1)
+
+        ! B_theta, B_phi from 1D spline at rho_tor
+        call evaluate_batch_splines_1d_der2(self%bcovar_spline, rho_tor, y1d_bc, &
+                                             dy1d_bc, d2y1d_bc)
+        B_theta_val = y1d_bc(1)
+        B_phi_val = y1d_bc(2)
+
+        ! Bmod from 3D spline at (rho_tor, theta_B, phi_B)
+        x_3d(1) = rho_tor
+        x_3d(2) = modulo(theta_B, 8.0_dp * atan(1.0_dp))
+        x_3d(3) = modulo(phi_B, 8.0_dp * atan(1.0_dp) / real(self%nfp, dp))
+
+        call evaluate_batch_splines_3d_no_der(self%bmod_spline, x_3d, y_bmod)
+        Bmod_val = y_bmod(1)
+
+        Bmod = Bmod_val
+
+        ! Covariant vector potential: A_s = 0 (gauge choice)
+        Acov(1) = 0.0_dp
+        Acov(2) = A_theta_val
+        Acov(3) = A_phi_val
+
+        ! Covariant unit field direction: h = B/|B|
+        hcov(1) = 0.0_dp
+        hcov(2) = B_theta_val / Bmod_val
+        hcov(3) = B_phi_val / Bmod_val
+
+        if (present(sqgBctr)) then
+            error stop "sqgBctr not implemented for boozer_chartmap_field_t"
+        end if
+    end subroutine boozer_chartmap_evaluate
+
+    subroutine boozer_chartmap_cleanup(self)
+        type(boozer_chartmap_field_t), intent(inout) :: self
+
+        if (self%initialized) then
+            call destroy_batch_splines_1d(self%aphi_spline)
+            call destroy_batch_splines_1d(self%bcovar_spline)
+            call destroy_batch_splines_3d(self%bmod_spline)
+            self%initialized = .false.
+        end if
+    end subroutine boozer_chartmap_cleanup
+
+    subroutine evaluate_batch_splines_3d_no_der(spl, x, y)
+        use interpolate, only: BatchSplineData3D, evaluate_batch_splines_3d
+        type(BatchSplineData3D), intent(in) :: spl
+        real(dp), intent(in) :: x(3)
+        real(dp), intent(out) :: y(:)
+
+        call evaluate_batch_splines_3d(spl, x, y)
+    end subroutine evaluate_batch_splines_3d_no_der
+
+    subroutine read_3d_reordered(ncid, varid, n1, n2, n3, arr)
+        !> Read a NetCDF variable with dims (zeta, theta, rho) into Fortran
+        !> array (n1=rho, n2=theta, n3=zeta). NF90 reverses dimension order
+        !> so the data lands directly in (rho, theta, zeta) layout.
+        integer, intent(in) :: ncid, varid, n1, n2, n3
+        real(dp), intent(out) :: arr(n1, n2, n3)
+
+        call check_nc(nf90_get_var(ncid, varid, arr), "get_var 3D")
+    end subroutine read_3d_reordered
+
+    subroutine check_nc(status, location)
+        integer, intent(in) :: status
+        character(len=*), intent(in) :: location
+
+        if (status /= nf90_noerr) then
+            print *, "NetCDF error at ", trim(location), ": ", &
+                trim(nf90_strerror(status))
+            error stop "NetCDF operation failed"
+        end if
+    end subroutine check_nc
+
+end module field_boozer_chartmap
diff --git a/src/field_can.f90 b/src/field_can.f90
index cd7efbd6..c3e1fa01 100644
--- a/src/field_can.f90
+++ b/src/field_can.f90
@@ -12,6 +12,7 @@ module field_can_mod
   integ_to_ref_meiss, ref_to_integ_meiss
 use field_can_albert, only : evaluate_albert, init_albert, integ_to_ref_albert, &
   ref_to_integ_albert
+use field_boozer_chartmap, only : boozer_chartmap_field_t
 
 implicit none
 
@@ -124,7 +125,8 @@ end function id_from_name
 
 subroutine init_field_can(field_id, field_noncan)
   use get_can_sub, only : get_canonical_coordinates, get_canonical_coordinates_with_field
-  use boozer_sub, only : get_boozer_coordinates, get_boozer_coordinates_with_field
+  use boozer_sub, only : get_boozer_coordinates, get_boozer_coordinates_with_field, &
+                         load_boozer_from_chartmap
   use field_can_meiss, only : get_meiss_coordinates
   use field_can_albert, only : get_albert_coordinates
 
@@ -140,7 +142,15 @@ subroutine init_field_can(field_id, field_noncan)
       case (CANFLUX)
         call get_canonical_coordinates_with_field(field_noncan)
       case (BOOZER)
-        call get_boozer_coordinates_with_field(field_noncan)
+        select type (f => field_noncan)
+        type is (boozer_chartmap_field_t)
+          call load_boozer_from_chartmap(f%filename)
+          ! Integrator coords = reference coords (both Boozer), no conversion
+          integ_to_ref => identity_transform
+          ref_to_integ => identity_transform
+        class default
+          call get_boozer_coordinates_with_field(field_noncan)
+        end select
       case (MEISS)
         call get_meiss_coordinates
       case (ALBERT)
diff --git a/src/simple_main.f90 b/src/simple_main.f90
index cd7d8842..dd793ed9 100644
--- a/src/simple_main.f90
+++ b/src/simple_main.f90
@@ -126,6 +126,7 @@ end subroutine main
     subroutine init_field(self, vmec_file, ans_s, ans_tp, amultharm, aintegmode)
         use field_base, only: magnetic_field_t
         use field, only: field_from_file
+        use field_boozer_chartmap, only: boozer_chartmap_field_t, is_boozer_chartmap
         use timing, only: print_phase_time
         use magfie_sub, only: TEST, CANFLUX, VMEC, BOOZER, MEISS, ALBERT, &
                               REFCOORDS, set_magfie_refcoords_field
@@ -140,13 +141,43 @@ subroutine init_field(self, vmec_file, ans_s, ans_tp, amultharm, aintegmode)
         integer, intent(in) :: ans_s, ans_tp, amultharm, aintegmode
         class(magnetic_field_t), allocatable :: field_temp
         character(:), allocatable :: vmec_equilibrium_file
+        logical :: use_boozer_chartmap
 
         self%integmode = aintegmode
 
+        ! Check if field_input is a Boozer chartmap (no VMEC needed)
+        use_boozer_chartmap = .false.
+        if (len_trim(field_input) > 0) then
+            use_boozer_chartmap = is_boozer_chartmap(field_input)
+        end if
+
         ! TEST field is analytic - no VMEC or field files needed
         if (isw_field_type == TEST) then
             self%fper = twopi  ! Full torus for analytic tokamak
             call print_phase_time('TEST field mode - no input files required')
+        else if (use_boozer_chartmap) then
+            ! Boozer chartmap: file-based, no VMEC initialization needed
+            call init_reference_coordinates(coord_input)
+            call print_phase_time('Reference coordinate system '// &
+                                  'initialization completed')
+
+            block
+                use libneo_coordinates, only: chartmap_coordinate_system_t
+                select type (cs => ref_coords)
+                type is (chartmap_coordinate_system_t)
+                    self%fper = twopi / real(cs%num_field_periods, dp)
+                class default
+                    self%fper = twopi
+                end select
+            end block
+
+            call init_stl_wall_if_enabled(coord_input)
+            call print_phase_time('STL wall initialization completed')
+
+            if (self%integmode >= 0) then
+                call field_from_file(field_input, field_temp)
+                call print_phase_time('Boozer chartmap field loading completed')
+            end if
         else
             vmec_equilibrium_file = select_vmec_equilibrium_file(vmec_file, &
                                                                  field_input, &
diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt
index d249240c..dc617c98 100644
--- a/test/tests/CMakeLists.txt
+++ b/test/tests/CMakeLists.txt
@@ -331,6 +331,24 @@ add_test(NAME test_array_utils COMMAND test_array_utils.x)
     add_test(NAME test_chartmap_metadata COMMAND test_chartmap_metadata.x)
     set_tests_properties(test_chartmap_metadata PROPERTIES LABELS "unit")
 
+    # Boozer chartmap test: generate test file with Python, then run Fortran test
+    add_custom_command(
+        OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/test_boozer_chartmap.nc
+        COMMAND ${Python_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/generate_test_boozer_chartmap.py
+        WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+        DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/generate_test_boozer_chartmap.py
+        COMMENT "Generating test Boozer chartmap NetCDF"
+    )
+    add_custom_target(generate_test_boozer_chartmap
+        DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/test_boozer_chartmap.nc)
+    add_executable(test_boozer_chartmap.x test_boozer_chartmap.f90)
+    target_link_libraries(test_boozer_chartmap.x simple)
+    add_dependencies(test_boozer_chartmap.x generate_test_boozer_chartmap)
+    add_test(NAME test_boozer_chartmap COMMAND test_boozer_chartmap.x)
+    set_tests_properties(test_boozer_chartmap PROPERTIES
+        WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+        LABELS "unit")
+
     if(SIMPLE_ENABLE_CGAL)
         add_executable(test_stl_wall_intersection.x test_stl_wall_intersection.f90)
         target_link_libraries(test_stl_wall_intersection.x simple)
diff --git a/test/tests/generate_test_boozer_chartmap.py b/test/tests/generate_test_boozer_chartmap.py
new file mode 100644
index 00000000..3cce2455
--- /dev/null
+++ b/test/tests/generate_test_boozer_chartmap.py
@@ -0,0 +1,103 @@
+#!/usr/bin/env python3
+"""Generate a minimal Boozer chartmap NetCDF test file with known analytic values."""
+
+import numpy as np
+import netCDF4
+
+def generate_test_file(filename):
+    """Create a test Boozer chartmap with simple analytic field.
+
+    Uses a tokamak-like field: Bmod ~ B0 * (1 - epsilon * cos(theta))
+    with B_theta ~ constant, B_phi ~ constant, A_phi ~ -torflux * iota * s.
+    """
+    nfp = 2
+    n_rho = 17
+    n_theta = 33
+    n_phi = 33
+    B0 = 2.0e4  # 2 T in Gauss
+    R0 = 150.0  # 1.5 m in cm
+    a = 50.0    # 0.5 m in cm
+    epsilon = a / R0
+    iota = 0.5
+    torflux = B0 * np.pi * a**2 / (2.0 * np.pi)  # approximate
+
+    twopi = 2.0 * np.pi
+
+    rho_grid = np.linspace(0.0, 1.0, n_rho)
+    rho_grid[0] = 1e-6
+    theta_grid = np.linspace(0.0, twopi, n_theta, endpoint=False)
+    zeta_grid = np.linspace(0.0, twopi / nfp, n_phi, endpoint=False)
+
+    # Geometry: concentric circular cross-section tokamak
+    X = np.zeros((n_rho, n_theta, n_phi))
+    Y = np.zeros((n_rho, n_theta, n_phi))
+    Z = np.zeros((n_rho, n_theta, n_phi))
+    Bmod = np.zeros((n_rho, n_theta, n_phi))
+
+    for ir in range(n_rho):
+        rho = rho_grid[ir]
+        r = a * rho  # minor radius
+        for it in range(n_theta):
+            theta = theta_grid[it]
+            R = R0 + r * np.cos(theta)
+            Zval = r * np.sin(theta)
+            for ip in range(n_phi):
+                phi = zeta_grid[ip]
+                X[ir, it, ip] = R * np.cos(phi)
+                Y[ir, it, ip] = R * np.sin(phi)
+                Z[ir, it, ip] = Zval
+                Bmod[ir, it, ip] = B0 / (1.0 + rho * epsilon * np.cos(theta))
+
+    # Profiles
+    A_phi = -torflux * iota * rho_grid**2
+    B_theta = np.full(n_rho, B0 * epsilon * iota)
+    B_phi = np.full(n_rho, B0 * R0)
+
+    # Write
+    ds = netCDF4.Dataset(filename, "w", format="NETCDF4")
+
+    ds.createDimension("rho", n_rho)
+    ds.createDimension("theta", n_theta)
+    ds.createDimension("zeta", n_phi)
+
+    v = ds.createVariable("rho", "f8", ("rho",))
+    v[:] = rho_grid
+
+    v = ds.createVariable("theta", "f8", ("theta",))
+    v[:] = theta_grid
+
+    v = ds.createVariable("zeta", "f8", ("zeta",))
+    v[:] = zeta_grid
+
+    for name, arr in [("x", X), ("y", Y), ("z", Z)]:
+        v = ds.createVariable(name, "f8", ("zeta", "theta", "rho"))
+        v[:] = np.transpose(arr, (2, 1, 0))
+        v.units = "cm"
+
+    v = ds.createVariable("A_phi", "f8", ("rho",))
+    v[:] = A_phi
+
+    v = ds.createVariable("B_theta", "f8", ("rho",))
+    v[:] = B_theta
+
+    v = ds.createVariable("B_phi", "f8", ("rho",))
+    v[:] = B_phi
+
+    v = ds.createVariable("Bmod", "f8", ("zeta", "theta", "rho"))
+    v[:] = np.transpose(Bmod, (2, 1, 0))
+
+    ds.num_field_periods = np.int32(nfp)
+    ds.rho_convention = "rho_tor"
+    ds.zeta_convention = "cyl"
+    ds.rho_lcfs = 1.0
+    ds.boozer_field = np.int32(1)
+    ds.torflux = torflux
+
+    ds.close()
+    print(f"Generated {filename}")
+    print(f"  torflux={torflux:.6e}, B0={B0:.1f}, R0={R0:.1f}, a={a:.1f}")
+    print(f"  nfp={nfp}, nrho={n_rho}, ntheta={n_theta}, nphi={n_phi}")
+
+
+if __name__ == "__main__":
+    generate_test_file("test_boozer_chartmap.nc")
diff --git a/test/tests/test_boozer_chartmap.f90 b/test/tests/test_boozer_chartmap.f90
new file mode 100644
index 00000000..0e46688c
--- /dev/null
+++ b/test/tests/test_boozer_chartmap.f90
@@ -0,0 +1,105 @@
+program test_boozer_chartmap
+    use, intrinsic :: iso_fortran_env, only: dp => real64
+    use field_boozer_chartmap, only: boozer_chartmap_field_t, &
+                                     create_boozer_chartmap_field, &
+                                     is_boozer_chartmap
+
+    implicit none
+
+    type(boozer_chartmap_field_t), allocatable :: field
+    real(dp) :: x(3), Acov(3), hcov(3), Bmod
+    character(len=256) :: filename
+    logical :: is_bc
+    integer :: nfail
+
+    nfail = 0
+    filename = 'test_boozer_chartmap.nc'
+
+    ! Test detection
+    is_bc = is_boozer_chartmap(filename)
+    if (.not. is_bc) then
+        print *, 'FAIL: is_boozer_chartmap returned false'
+        nfail = nfail + 1
+    else
+        print *, 'PASS: is_boozer_chartmap detected file correctly'
+    end if
+
+    ! Test loading
+    call create_boozer_chartmap_field(filename, field)
+
+    if (.not. field%initialized) then
+        print *, 'FAIL: field not initialized'
+        error stop
+    end if
+
+    if (field%nfp /= 2) then
+        print *, 'FAIL: nfp =', field%nfp, ' expected 2'
+        nfail = nfail + 1
+    else
+        print *, 'PASS: nfp = 2'
+    end if
+
+    if (field%torflux <= 0.0_dp) then
+        print *, 'FAIL: torflux =', field%torflux, ' expected > 0'
+        nfail = nfail + 1
+    else
+        print *, 'PASS: torflux =', field%torflux
+    end if
+
+    ! Test field evaluation at mid-radius, theta=0, phi=0
+    x = [0.25_dp, 0.0_dp, 0.0_dp]  ! s=0.25, theta=0, phi=0
+    call field%evaluate(x, Acov, hcov, Bmod)
+
+    ! Bmod should be positive and reasonable (analytic tokamak ~ B0)
+    if (Bmod <= 0.0_dp) then
+        print *, 'FAIL: Bmod =', Bmod, ' expected > 0'
+        nfail = nfail + 1
+    else
+        print *, 'PASS: Bmod =', Bmod
+    end if
+
+    ! A_theta = torflux * s = torflux * 0.25
+    if (abs(Acov(2) - field%torflux * 0.25_dp) > 1.0e-8_dp * abs(Acov(2))) then
+        print *, 'FAIL: Acov(2) =', Acov(2), ' expected', field%torflux * 0.25_dp
+        nfail = nfail + 1
+    else
+        print *, 'PASS: Acov(2) = torflux * s'
+    end if
+
+    ! hcov(1) should be 0 (no radial component in Boozer)
+    if (abs(hcov(1)) > 1.0e-15_dp) then
+        print *, 'FAIL: hcov(1) =', hcov(1), ' expected 0'
+        nfail = nfail + 1
+    else
+        print *, 'PASS: hcov(1) = 0'
+    end if
+
+    ! hcov(2) and hcov(3) should be B_theta/Bmod and B_phi/Bmod
+    if (abs(hcov(2)) < 1.0e-15_dp .or. abs(hcov(3)) < 1.0e-15_dp) then
+        print *, 'FAIL: hcov(2) or hcov(3) is zero'
+        nfail = nfail + 1
+    else
+        print *, 'PASS: hcov(2) =', hcov(2), ' hcov(3) =', hcov(3)
+    end if
+
+    ! Test at another point
+    x = [0.5_dp, 1.0_dp, 0.5_dp]  ! s=0.5, theta=1, phi=0.5
+    call field%evaluate(x, Acov, hcov, Bmod)
+
+    if (Bmod <= 0.0_dp) then
+        print *, 'FAIL: Bmod at second point =', Bmod
+        nfail = nfail + 1
+    else
+        print *, 'PASS: Bmod at second point =', Bmod
+    end if
+
+    ! Summary
+    print *, ''
+    if (nfail == 0) then
+        print *, 'All tests passed.'
+    else
+        print *, nfail, ' tests failed.'
+        error stop 'Test failures'
+    end if
+
+end program test_boozer_chartmap

From 47c6abebcc3626349ad1711643b791d0b27ad3be Mon Sep 17 00:00:00 2001
From: Christopher Albert <albert@tugraz.at>
Date: Thu, 26 Mar 2026 23:53:57 +0100
Subject: [PATCH 2/9] fix: coordinate convention and chartmap validation for
 boozer chartmap

- evaluate() now takes x(1)=rho (matching chartmap coordinate system)
  instead of x(1)=s, with internal rho-to-s conversion for A_phi spline
- Add direct-evaluation fast path in field_splined.f90 for
  boozer_chartmap_field_t, avoiding unnecessary Cartesian roundtrip
- Store num_field_periods as a scalar NetCDF variable (not attribute)
  to match libneo chartmap validator expectations
- Fix test to use rho coordinates
---
 python/pysimple/boozer_chartmap.py          |  5 +++-
 src/boozer_converter.F90                    |  8 +++--
 src/field/field_boozer_chartmap.f90         | 33 ++++++++++++++-------
 src/field/field_splined.f90                 |  9 ++++++
 test/tests/generate_test_boozer_chartmap.py |  5 +++-
 test/tests/test_boozer_chartmap.f90         | 11 +++----
 6 files changed, 50 insertions(+), 21 deletions(-)

diff --git a/python/pysimple/boozer_chartmap.py b/python/pysimple/boozer_chartmap.py
index 858b3811..a340c1fa 100644
--- a/python/pysimple/boozer_chartmap.py
+++ b/python/pysimple/boozer_chartmap.py
@@ -299,8 +299,11 @@ def generate_boozer_chartmap(wout_file, output_file, n_rho=62, n_theta=63, n_phi
     v = ds.createVariable("Bmod", "f8", ("zeta", "theta", "rho"))
     v[:] = np.transpose(Bmod, (2, 1, 0))
 
+    # num_field_periods must be a scalar variable (not attribute) for libneo
+    v = ds.createVariable("num_field_periods", "i4")
+    v[:] = np.int32(nfp)
+
     # Global attributes
-    ds.num_field_periods = np.int32(nfp)
     ds.rho_convention = "rho_tor"
     ds.zeta_convention = "cyl"
     ds.rho_lcfs = rho_lcfs
diff --git a/src/boozer_converter.F90 b/src/boozer_converter.F90
index 447ecd79..e316b9b3 100644
--- a/src/boozer_converter.F90
+++ b/src/boozer_converter.F90
@@ -930,11 +930,13 @@ subroutine load_boozer_from_chartmap(filename)
         call nc_check(nf90_inq_varid(ncid, "zeta", varid), "var zeta")
         call nc_check(nf90_get_var(ncid, varid, zeta), "get zeta")
 
-        ! Read scalar attributes
+        ! Read scalar attributes and variables
         call nc_check(nf90_get_att(ncid, nf90_global, "torflux", torflux_val), &
                        "att torflux")
-        call nc_check(nf90_get_att(ncid, nf90_global, "num_field_periods", nfp_file), &
-                       "att num_field_periods")
+        call nc_check(nf90_inq_varid(ncid, "num_field_periods", varid), &
+                       "var num_field_periods")
+        call nc_check(nf90_get_var(ncid, varid, nfp_file), &
+                       "get num_field_periods")
 
         ! Read 1D profiles
         allocate (A_phi_arr(n_rho), B_theta_arr(n_rho), B_phi_arr(n_rho))
diff --git a/src/field/field_boozer_chartmap.f90 b/src/field/field_boozer_chartmap.f90
index e4c523d5..4be1e65a 100644
--- a/src/field/field_boozer_chartmap.f90
+++ b/src/field/field_boozer_chartmap.f90
@@ -110,11 +110,13 @@ subroutine create_boozer_chartmap_field(filename, field)
         call check_nc(nf90_inq_varid(ncid, "zeta", varid), "inq_var zeta")
         call check_nc(nf90_get_var(ncid, varid, zeta), "get_var zeta")
 
-        ! Read scalar attributes
+        ! Read scalar attributes and variables
         call check_nc(nf90_get_att(ncid, nf90_global, "torflux", torflux_val), &
                        "get_att torflux")
-        call check_nc(nf90_get_att(ncid, nf90_global, "num_field_periods", nfp_file), &
-                       "get_att num_field_periods")
+        call check_nc(nf90_inq_varid(ncid, "num_field_periods", varid), &
+                       "inq_var num_field_periods")
+        call check_nc(nf90_get_var(ncid, varid, nfp_file), &
+                       "get_var num_field_periods")
 
         ! Read 1D profiles
         allocate (A_phi_arr(n_rho), B_theta_arr(n_rho), B_phi_arr(n_rho))
@@ -196,26 +198,31 @@ subroutine create_boozer_chartmap_field(filename, field)
     end subroutine create_boozer_chartmap_field
 
     subroutine boozer_chartmap_evaluate(self, x, Acov, hcov, Bmod, sqgBctr)
-        !> Evaluate field at x = (s, theta_B, phi_B) in Boozer coordinates.
-        !> Returns covariant components for the non-canonical path.
+        !> Evaluate field at x = (rho, theta_B, phi_B) in chartmap coordinates.
+        !> Input x(1) = rho = sqrt(s), matching the chartmap coordinate system.
+        !> Returns covariant components in (rho, theta_B, phi_B) coordinates.
         class(boozer_chartmap_field_t), intent(in) :: self
         real(dp), intent(in) :: x(3)
         real(dp), intent(out) :: Acov(3), hcov(3), Bmod
         real(dp), intent(out), optional :: sqgBctr(3)
 
-        real(dp) :: s, rho_tor, theta_B, phi_B
+        real(dp) :: rho_tor, s, theta_B, phi_B
         real(dp) :: A_phi_val, A_theta_val
         real(dp) :: y1d_aphi(1), dy1d_aphi(1), d2y1d_aphi(1)
         real(dp) :: y1d_bc(2), dy1d_bc(2), d2y1d_bc(2)
         real(dp) :: x_3d(3), y_bmod(1)
         real(dp) :: B_theta_val, B_phi_val, Bmod_val
+        real(dp), parameter :: twopi = 8.0_dp * atan(1.0_dp)
 
-        s = x(1)
+        rho_tor = x(1)
         theta_B = x(2)
         phi_B = x(3)
-        rho_tor = sqrt(max(s, 0.0_dp))
+        s = rho_tor**2
 
         ! A_theta = torflux * s (linear, by definition)
+        ! In rho coordinates: A_theta(rho) = torflux * rho^2
+        ! Covariant component transforms: A_rho = A_s * ds/drho = A_s * 2*rho
+        ! But A_s = 0 (gauge), so A_rho = 0 regardless.
         A_theta_val = self%torflux * s
 
         ! A_phi from 1D spline at s
@@ -231,20 +238,24 @@ subroutine boozer_chartmap_evaluate(self, x, Acov, hcov, Bmod, sqgBctr)
 
         ! Bmod from 3D spline at (rho_tor, theta_B, phi_B)
         x_3d(1) = rho_tor
-        x_3d(2) = modulo(theta_B, 8.0_dp * atan(1.0_dp))
-        x_3d(3) = modulo(phi_B, 8.0_dp * atan(1.0_dp) / real(self%nfp, dp))
+        x_3d(2) = modulo(theta_B, twopi)
+        x_3d(3) = modulo(phi_B, twopi / real(self%nfp, dp))
 
         call evaluate_batch_splines_3d_no_der(self%bmod_spline, x_3d, y_bmod)
         Bmod_val = y_bmod(1)
 
         Bmod = Bmod_val
 
-        ! Covariant vector potential: A_s = 0 (gauge choice)
+        ! Covariant vector potential in (rho, theta, phi) coordinates.
+        ! A_theta and A_phi are the same in (s, theta, phi) and (rho, theta, phi)
+        ! since the angular coordinates are unchanged. A_rho = 0.
         Acov(1) = 0.0_dp
         Acov(2) = A_theta_val
         Acov(3) = A_phi_val
 
         ! Covariant unit field direction: h = B/|B|
+        ! B_theta and B_phi are covariant and independent of the radial
+        ! coordinate choice (they refer to the angular basis vectors).
         hcov(1) = 0.0_dp
         hcov(2) = B_theta_val / Bmod_val
         hcov(3) = B_phi_val / Bmod_val
diff --git a/src/field/field_splined.f90 b/src/field/field_splined.f90
index 4c5655ee..8909a25f 100644
--- a/src/field/field_splined.f90
+++ b/src/field/field_splined.f90
@@ -30,6 +30,7 @@ module field_splined
                            destroy_batch_splines_3d
     use field_base, only: magnetic_field_t
     use field_vmec, only: vmec_field_t
+    use field_boozer_chartmap, only: boozer_chartmap_field_t
     use libneo_coordinates, only: coordinate_system_t, vmec_coordinate_system_t, &
                                   chartmap_coordinate_system_t, RHO_TOR, RHO_POL, &
                                   PSI_TOR_NORM, PSI_POL_NORM
@@ -461,6 +462,14 @@ subroutine evaluate_at_ref_coords(source, ref_coords, x_scaled, x_ref, x_cart, &
             end select
             needs_scaling = .true.
             return
+        type is (boozer_chartmap_field_t)
+            ! Boozer chartmap evaluates directly in (rho, theta_B, phi_B),
+            ! which matches the chartmap reference coordinate system.
+            call source%evaluate(x_ref, Acov, hcov, Bmod)
+            x_cart = 0d0
+            e_cov = 0d0
+            needs_scaling = .true.
+            return
         class default
             continue
         end select
diff --git a/test/tests/generate_test_boozer_chartmap.py b/test/tests/generate_test_boozer_chartmap.py
index 3cce2455..b2b80fe3 100644
--- a/test/tests/generate_test_boozer_chartmap.py
+++ b/test/tests/generate_test_boozer_chartmap.py
@@ -86,7 +86,10 @@ def generate_test_file(filename):
     v = ds.createVariable("Bmod", "f8", ("zeta", "theta", "rho"))
     v[:] = np.transpose(Bmod, (2, 1, 0))
 
-    ds.num_field_periods = np.int32(nfp)
+    # num_field_periods must be a scalar variable (not attribute) for libneo
+    v = ds.createVariable("num_field_periods", "i4")
+    v[:] = np.int32(nfp)
+
     ds.rho_convention = "rho_tor"
     ds.zeta_convention = "cyl"
     ds.rho_lcfs = 1.0
diff --git a/test/tests/test_boozer_chartmap.f90 b/test/tests/test_boozer_chartmap.f90
index 0e46688c..3eebc808 100644
--- a/test/tests/test_boozer_chartmap.f90
+++ b/test/tests/test_boozer_chartmap.f90
@@ -47,7 +47,8 @@ program test_boozer_chartmap
     end if
 
     ! Test field evaluation at mid-radius, theta=0, phi=0
-    x = [0.25_dp, 0.0_dp, 0.0_dp]  ! s=0.25, theta=0, phi=0
+    ! x(1) = rho = sqrt(s), so rho=0.5 means s=0.25
+    x = [0.5_dp, 0.0_dp, 0.0_dp]  ! rho=0.5, theta=0, phi=0
     call field%evaluate(x, Acov, hcov, Bmod)
 
     ! Bmod should be positive and reasonable (analytic tokamak ~ B0)
@@ -58,12 +59,12 @@ program test_boozer_chartmap
         print *, 'PASS: Bmod =', Bmod
     end if
 
-    ! A_theta = torflux * s = torflux * 0.25
+    ! A_theta = torflux * s = torflux * rho^2 = torflux * 0.25
     if (abs(Acov(2) - field%torflux * 0.25_dp) > 1.0e-8_dp * abs(Acov(2))) then
         print *, 'FAIL: Acov(2) =', Acov(2), ' expected', field%torflux * 0.25_dp
         nfail = nfail + 1
     else
-        print *, 'PASS: Acov(2) = torflux * s'
+        print *, 'PASS: Acov(2) = torflux * rho^2'
     end if
 
     ! hcov(1) should be 0 (no radial component in Boozer)
@@ -82,8 +83,8 @@ program test_boozer_chartmap
         print *, 'PASS: hcov(2) =', hcov(2), ' hcov(3) =', hcov(3)
     end if
 
-    ! Test at another point
-    x = [0.5_dp, 1.0_dp, 0.5_dp]  ! s=0.5, theta=1, phi=0.5
+    ! Test at another point: rho=0.7, theta=1, phi=0.5
+    x = [0.7_dp, 1.0_dp, 0.5_dp]
     call field%evaluate(x, Acov, hcov, Bmod)
 
     if (Bmod <= 0.0_dp) then

From 0cc316d54244eb81b763dd23fd2fb6283bb7bb19 Mon Sep 17 00:00:00 2001
From: Christopher Albert <albert@tugraz.at>
Date: Fri, 27 Mar 2026 00:55:46 +0100
Subject: [PATCH 3/9] fix: proper rho/s coordinate transform for Boozer
 chartmap
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The Boozer integrator uses s as the radial coordinate, but chartmap
reference coordinates use rho = sqrt(s). The identity transform was
wrong — particles at rho=0.55 (s=0.3) were being integrated at s=0.55,
causing spurious losses.

Add integ_to_ref_boozer_chartmap (s -> rho) and
ref_to_integ_boozer_chartmap (rho -> s) transforms.

E2E test now shows exact match in total confined fraction between
VMEC-Boozer and chartmap-only paths.
---
 src/field/field_can_boozer.f90 | 24 ++++++++++++++++++++++++
 src/field_can.f90              |  9 +++++----
 2 files changed, 29 insertions(+), 4 deletions(-)

diff --git a/src/field/field_can_boozer.f90 b/src/field/field_can_boozer.f90
index bc2ea7c8..93b948ea 100644
--- a/src/field/field_can_boozer.f90
+++ b/src/field/field_can_boozer.f90
@@ -4,8 +4,32 @@ module field_can_boozer
 
     implicit none
 
+    public :: integ_to_ref_boozer_chartmap, ref_to_integ_boozer_chartmap
+
 contains
 
+    subroutine integ_to_ref_boozer_chartmap(xinteg, xref)
+        !> Convert integrator coords (s, theta_B, phi_B) to
+        !> chartmap reference coords (rho, theta_B, phi_B).
+        real(dp), intent(in) :: xinteg(3)
+        real(dp), intent(out) :: xref(3)
+
+        xref(1) = sqrt(max(xinteg(1), 0.0_dp))
+        xref(2) = xinteg(2)
+        xref(3) = xinteg(3)
+    end subroutine integ_to_ref_boozer_chartmap
+
+    subroutine ref_to_integ_boozer_chartmap(xref, xinteg)
+        !> Convert chartmap reference coords (rho, theta_B, phi_B) to
+        !> integrator coords (s, theta_B, phi_B).
+        real(dp), intent(in) :: xref(3)
+        real(dp), intent(out) :: xinteg(3)
+
+        xinteg(1) = xref(1)**2
+        xinteg(2) = xref(2)
+        xinteg(3) = xref(3)
+    end subroutine ref_to_integ_boozer_chartmap
+
     subroutine evaluate_boozer(f, r, th_c, ph_c, mode_secders)
         type(field_can_t), intent(inout) :: f
         real(dp), intent(in) :: r, th_c, ph_c
diff --git a/src/field_can.f90 b/src/field_can.f90
index c3e1fa01..d58607fd 100644
--- a/src/field_can.f90
+++ b/src/field_can.f90
@@ -7,7 +7,8 @@ module field_can_mod
   identity_transform, field_can_t
 use field_can_test, only : evaluate_test
 use field_can_flux, only : evaluate_flux, integ_to_ref_flux, ref_to_integ_flux
-use field_can_boozer, only : evaluate_boozer, integ_to_ref_boozer, ref_to_integ_boozer
+use field_can_boozer, only : evaluate_boozer, integ_to_ref_boozer, ref_to_integ_boozer, &
+  integ_to_ref_boozer_chartmap, ref_to_integ_boozer_chartmap
 use field_can_meiss, only : init_meiss, evaluate_meiss, &
   integ_to_ref_meiss, ref_to_integ_meiss
 use field_can_albert, only : evaluate_albert, init_albert, integ_to_ref_albert, &
@@ -145,9 +146,9 @@ subroutine init_field_can(field_id, field_noncan)
         select type (f => field_noncan)
         type is (boozer_chartmap_field_t)
           call load_boozer_from_chartmap(f%filename)
-          ! Integrator coords = reference coords (both Boozer), no conversion
-          integ_to_ref => identity_transform
-          ref_to_integ => identity_transform
+          ! Chartmap ref coords use rho=sqrt(s), integrator uses s
+          integ_to_ref => integ_to_ref_boozer_chartmap
+          ref_to_integ => ref_to_integ_boozer_chartmap
         class default
           call get_boozer_coordinates_with_field(field_noncan)
         end select

From 657235c2e129741db9f62193d7b1ebcf139b6d93 Mon Sep 17 00:00:00 2001
From: Christopher Albert <albert@tugraz.at>
Date: Fri, 27 Mar 2026 07:51:07 +0100
Subject: [PATCH 4/9] fix: initialize bmin/bmax in chartmap mode via
 init_starting_surf

The passing/trapped classification needs bmin and bmax from field line
tracing (init_starting_surf). In chartmap mode this was skipped, leaving
bmin=bmax=0 and producing garbage classification.

Fix: call init_magfie(BOOZER) first to set up magfie_boozer (which works
with the chartmap splines), then call init_starting_surf normally.

E2E test now shows exact match in total, passing, AND trapped fractions.
---
 src/simple_main.f90 | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/src/simple_main.f90 b/src/simple_main.f90
index dd793ed9..c430c675 100644
--- a/src/simple_main.f90
+++ b/src/simple_main.f90
@@ -91,6 +91,21 @@ subroutine main
             call print_phase_time('TEST field initialization completed')
             call sample_particles_test_field
             call print_phase_time('TEST field particle sampling completed')
+        else if (chartmap_mode) then
+            ! Boozer chartmap: no VMEC, use Boozer magfie for field line tracing.
+            if (startmode /= 2) then
+                error stop 'Boozer chartmap requires startmode=2 '// &
+                    '(load particles from start.dat)'
+            end if
+
+            call init_magfie(isw_field_type)
+            call print_phase_time('Boozer magfie initialization completed')
+
+            call init_starting_surf
+            call print_phase_time('Starting surface initialization completed')
+
+            call sample_particles
+            call print_phase_time('Particle sampling completed')
         else
             call init_magfie(VMEC)
             call print_phase_time('VMEC magnetic field initialization completed')

From 709df2af9f6e19bab07f5672e060371fb69e995d Mon Sep 17 00:00:00 2001
From: Christopher Albert <albert@tugraz.at>
Date: Sat, 28 Mar 2026 11:47:34 +0100
Subject: [PATCH 5/9] fix: address review findings in core Boozer chartmap path

---
 python/pysimple/boozer_chartmap.py  | 13 ++++++++++---
 src/field/field_boozer_chartmap.f90 | 13 ++-----------
 test/tests/test_boozer_chartmap.f90 | 27 +++++++++++++++++++++++++++
 3 files changed, 39 insertions(+), 14 deletions(-)

diff --git a/python/pysimple/boozer_chartmap.py b/python/pysimple/boozer_chartmap.py
index a340c1fa..000d47ad 100644
--- a/python/pysimple/boozer_chartmap.py
+++ b/python/pysimple/boozer_chartmap.py
@@ -1,9 +1,16 @@
 #!/usr/bin/env python3
 """Generate extended Boozer chartmap NetCDF from VMEC wout.nc.
 
-Uses SIMPLE's existing Boozer converter to compute the transformation,
-then writes the result as an extended chartmap file that can be used
-as both coord_input and field_input without VMEC or GVEC libraries.
+WARNING: This script uses an APPROXIMATE Boozer transformation based on
+VMEC lambda. For production use, prefer the Fortran export_boozer_chartmap
+subroutine (in boozer_converter.F90) which uses the exact Boozer coordinate
+transformation computed by get_boozer_coordinates(). The Fortran exporter
+also produces endpoint-included field grids matching the internal spline
+representation exactly.
+
+This Python script is a prototype for offline file generation and testing.
+The geometry (X, Y, Z) and Bmod are evaluated at VMEC angles, not true
+Boozer angles, leading to small but nonzero errors.
 
 Usage:
     python -m pysimple.boozer_chartmap --wout wout.nc --out boozer_chartmap.nc \\
diff --git a/src/field/field_boozer_chartmap.f90 b/src/field/field_boozer_chartmap.f90
index 4be1e65a..17cb27af 100644
--- a/src/field/field_boozer_chartmap.f90
+++ b/src/field/field_boozer_chartmap.f90
@@ -17,7 +17,7 @@ module field_boozer_chartmap
     use interpolate, only: BatchSplineData1D, BatchSplineData3D, &
                            construct_batch_splines_1d, construct_batch_splines_3d, &
                            evaluate_batch_splines_1d_der2, &
-                           evaluate_batch_splines_3d_der, &
+                           evaluate_batch_splines_3d, &
                            destroy_batch_splines_1d, destroy_batch_splines_3d
     use netcdf
 
@@ -241,7 +241,7 @@ subroutine boozer_chartmap_evaluate(self, x, Acov, hcov, Bmod, sqgBctr)
         x_3d(2) = modulo(theta_B, twopi)
         x_3d(3) = modulo(phi_B, twopi / real(self%nfp, dp))
 
-        call evaluate_batch_splines_3d_no_der(self%bmod_spline, x_3d, y_bmod)
+        call evaluate_batch_splines_3d(self%bmod_spline, x_3d, y_bmod)
         Bmod_val = y_bmod(1)
 
         Bmod = Bmod_val
@@ -276,15 +276,6 @@ subroutine boozer_chartmap_cleanup(self)
         end if
     end subroutine boozer_chartmap_cleanup
 
-    subroutine evaluate_batch_splines_3d_no_der(spl, x, y)
-        use interpolate, only: BatchSplineData3D, evaluate_batch_splines_3d
-        type(BatchSplineData3D), intent(in) :: spl
-        real(dp), intent(in) :: x(3)
-        real(dp), intent(out) :: y(:)
-
-        call evaluate_batch_splines_3d(spl, x, y)
-    end subroutine evaluate_batch_splines_3d_no_der
-
     subroutine read_3d_reordered(ncid, varid, n1, n2, n3, arr)
         !> Read a NetCDF variable with dims (zeta, theta, rho) into Fortran
         !> array (n1=rho, n2=theta, n3=zeta). NF90 reverses dimension order
diff --git a/test/tests/test_boozer_chartmap.f90 b/test/tests/test_boozer_chartmap.f90
index 3eebc808..7459d174 100644
--- a/test/tests/test_boozer_chartmap.f90
+++ b/test/tests/test_boozer_chartmap.f90
@@ -94,6 +94,33 @@ program test_boozer_chartmap
         print *, 'PASS: Bmod at second point =', Bmod
     end if
 
+    ! Test integ_to_ref / ref_to_integ coordinate transforms
+    block
+        use field_can_boozer, only: integ_to_ref_boozer_chartmap, &
+                                    ref_to_integ_boozer_chartmap
+        real(dp) :: x_ref(3), x_integ(3), x_roundtrip(3), err
+
+        x_ref = [0.7_dp, 2.5_dp, 1.3_dp]  ! rho, theta, phi
+        call ref_to_integ_boozer_chartmap(x_ref, x_integ)
+
+        ! x_integ(1) should be rho^2 = 0.49
+        if (abs(x_integ(1) - 0.49_dp) > 1.0e-14_dp) then
+            print *, 'FAIL: ref_to_integ s =', x_integ(1), ' expected 0.49'
+            nfail = nfail + 1
+        else
+            print *, 'PASS: ref_to_integ rho->s correct'
+        end if
+
+        call integ_to_ref_boozer_chartmap(x_integ, x_roundtrip)
+        err = maxval(abs(x_roundtrip - x_ref))
+        if (err > 1.0e-14_dp) then
+            print *, 'FAIL: coord transform roundtrip error =', err
+            nfail = nfail + 1
+        else
+            print *, 'PASS: coord transform roundtrip < 1e-14'
+        end if
+    end block
+
     ! Summary
     print *, ''
     if (nfail == 0) then

From bd69aceaf8d2b097c686d2da67fcf62f1f02d58f Mon Sep 17 00:00:00 2001
From: Christopher Albert <albert@tugraz.at>
Date: Sat, 28 Mar 2026 11:53:36 +0100
Subject: [PATCH 6/9] fix: restore chartmap-mode detection in core branch

---
 src/simple_main.f90 | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/src/simple_main.f90 b/src/simple_main.f90
index c430c675..9f23ac34 100644
--- a/src/simple_main.f90
+++ b/src/simple_main.f90
@@ -41,17 +41,19 @@ module simple_main
     subroutine main
         use params, only: read_config, netcdffile, ns_s, ns_tp, multharm, &
                           integmode, params_init, swcoll, generate_start_only, &
-                          isw_field_type, &
+                          isw_field_type, field_input, startmode, &
                           ntestpart, ntimstep, coord_input
         use timing, only: init_timer, print_phase_time
         use magfie_sub, only: TEST, VMEC, init_magfie
         use samplers, only: init_starting_surf
         use version, only: simple_version
+        use field_boozer_chartmap, only: is_boozer_chartmap
 
         implicit none
 
         character(256) :: config_file
         type(tracer_t) :: norb
+        logical :: chartmap_mode
 
         ! Print version on startup
         print '(A,A)', 'SIMPLE version ', simple_version
@@ -85,6 +87,11 @@ subroutine main
             call print_phase_time('Collision initialization completed')
         end if
 
+        chartmap_mode = .false.
+        if (len_trim(field_input) > 0) then
+            chartmap_mode = is_boozer_chartmap(field_input)
+        end if
+
         if (isw_field_type == TEST) then
             ! TEST field uses analytic tokamak - no VMEC needed for sampling
             call init_magfie(TEST)

From e23f1c04f6445e6a6419c244cfc7a6d57f8cb869 Mon Sep 17 00:00:00 2001
From: Christopher Albert <albert@tugraz.at>
Date: Sat, 28 Mar 2026 14:07:44 +0100
Subject: [PATCH 7/9] chore: remove dead approximate Boozer exporter

---
 python/pysimple/boozer_chartmap.py | 339 -----------------------------
 1 file changed, 339 deletions(-)
 delete mode 100644 python/pysimple/boozer_chartmap.py

diff --git a/python/pysimple/boozer_chartmap.py b/python/pysimple/boozer_chartmap.py
deleted file mode 100644
index 000d47ad..00000000
--- a/python/pysimple/boozer_chartmap.py
+++ /dev/null
@@ -1,339 +0,0 @@
-#!/usr/bin/env python3
-"""Generate extended Boozer chartmap NetCDF from VMEC wout.nc.
-
-WARNING: This script uses an APPROXIMATE Boozer transformation based on
-VMEC lambda. For production use, prefer the Fortran export_boozer_chartmap
-subroutine (in boozer_converter.F90) which uses the exact Boozer coordinate
-transformation computed by get_boozer_coordinates(). The Fortran exporter
-also produces endpoint-included field grids matching the internal spline
-representation exactly.
-
-This Python script is a prototype for offline file generation and testing.
-The geometry (X, Y, Z) and Bmod are evaluated at VMEC angles, not true
-Boozer angles, leading to small but nonzero errors.
-
-Usage:
-    python -m pysimple.boozer_chartmap --wout wout.nc --out boozer_chartmap.nc \\
-        [--nrho 62] [--ntheta 63] [--nphi 64]
-"""
-
-import argparse
-import numpy as np
-
-try:
-    import netCDF4
-except ImportError:
-    raise ImportError("netCDF4 is required: pip install netCDF4")
-
-
-def read_vmec_wout(filename):
-    """Read essential data from VMEC wout.nc."""
-    ds = netCDF4.Dataset(filename, "r")
-
-    nfp = int(ds.variables["nfp"][:])
-    ns = int(ds.dimensions["radius"].size)
-    xm = ds.variables["xm"][:]
-    xn = ds.variables["xn"][:]
-    rmnc = ds.variables["rmnc"][:]
-    zmns = ds.variables["zmns"][:]
-
-    # Check for non-stellarator-symmetric terms
-    has_rbc_sin = "rmns" in ds.variables
-    rmns = ds.variables["rmns"][:] if has_rbc_sin else None
-    zmnc = ds.variables["zmnc"][:] if has_rbc_sin else None
-
-    # Profiles
-    phi_edge = float(ds.variables["phi"][-1])  # total toroidal flux
-    iotaf = ds.variables["iotaf"][:] if "iotaf" in ds.variables else None
-    iotas = ds.variables["iotas"][:] if "iotas" in ds.variables else None
-
-    # Vector potential: bvco = B_phi (toroidal covariant), buco = B_theta (poloidal covariant)
-    bvco = ds.variables["bvco"][:] if "bvco" in ds.variables else None
-    buco = ds.variables["buco"][:] if "buco" in ds.variables else None
-
-    # Boozer spectrum if available
-    has_boozer = "bmnc_b" in ds.variables
-    bmnc_b = ds.variables["bmnc_b"][:] if has_boozer else None
-    xm_b = ds.variables["ixm_b"][:] if has_boozer else None
-    xn_b = ds.variables["ixn_b"][:] if has_boozer else None
-
-    # Lambda (straight field line angle correction)
-    lmns = ds.variables["lmns"][:] if "lmns" in ds.variables else None
-    lmnc = ds.variables["lmnc"][:] if has_rbc_sin and "lmnc" in ds.variables else None
-
-    ds.close()
-
-    return {
-        "nfp": nfp,
-        "ns": ns,
-        "xm": xm,
-        "xn": xn,
-        "rmnc": rmnc,
-        "zmns": zmns,
-        "rmns": rmns,
-        "zmnc": zmnc,
-        "phi_edge": phi_edge,
-        "iotaf": iotaf,
-        "iotas": iotas,
-        "bvco": bvco,
-        "buco": buco,
-        "lmns": lmns,
-        "lmnc": lmnc,
-        "has_rbc_sin": has_rbc_sin,
-        "has_boozer": has_boozer,
-        "bmnc_b": bmnc_b,
-        "xm_b": xm_b,
-        "xn_b": xn_b,
-    }
-
-
-def vmec_to_boozer_angles(vmec, s_idx, theta_v, phi_v):
-    """Compute Boozer angles from VMEC angles at a given flux surface.
-
-    Returns theta_B, phi_B and the angle corrections.
-    Uses the lambda transformation: theta_B = theta_V + lambda * iota,
-    phi_B = phi_V + lambda.
-    This is approximate -- for a proper transformation, use booz_xform.
-    """
-    xm = vmec["xm"]
-    xn = vmec["xn"]
-    lmns = vmec["lmns"]
-    lmnc = vmec["lmnc"]
-    iota = vmec["iotaf"][s_idx] if vmec["iotaf"] is not None else vmec["iotas"][s_idx]
-
-    n_theta = len(theta_v)
-    n_phi = len(phi_v)
-    lam = np.zeros((n_theta, n_phi))
-
-    for k in range(len(xm)):
-        m = xm[k]
-        n = xn[k]
-        for it in range(n_theta):
-            for ip in range(n_phi):
-                angle = m * theta_v[it] - n * phi_v[ip]
-                lam[it, ip] += lmns[s_idx, k] * np.sin(angle)
-                if lmnc is not None:
-                    lam[it, ip] += lmnc[s_idx, k] * np.cos(angle)
-
-    return lam
-
-
-def evaluate_vmec_position(vmec, s_idx, theta, phi):
-    """Evaluate R, Z at given (theta, phi) on flux surface s_idx."""
-    xm = vmec["xm"]
-    xn = vmec["xn"]
-    rmnc = vmec["rmnc"]
-    zmns = vmec["zmns"]
-    rmns = vmec["rmns"]
-    zmnc = vmec["zmnc"]
-
-    R = np.zeros_like(theta)
-    Z = np.zeros_like(theta)
-
-    for k in range(len(xm)):
-        m = xm[k]
-        n = xn[k]
-        angle = m * theta - n * phi
-        R += rmnc[s_idx, k] * np.cos(angle)
-        Z += zmns[s_idx, k] * np.sin(angle)
-        if rmns is not None:
-            R += rmns[s_idx, k] * np.sin(angle)
-        if zmnc is not None:
-            Z += zmnc[s_idx, k] * np.cos(angle)
-
-    return R, Z
-
-
-def evaluate_vmec_bmod(vmec, s_idx, theta, phi):
-    """Evaluate |B| at given (theta, phi) on flux surface s_idx.
-
-    Uses the Fourier representation of |B| from VMEC.
-    """
-    if vmec["has_boozer"]:
-        # Use Boozer spectrum
-        xm_b = vmec["xm_b"]
-        xn_b = vmec["xn_b"]
-        bmnc_b = vmec["bmnc_b"]
-        bmod = np.zeros_like(theta)
-        for k in range(len(xm_b)):
-            m = xm_b[k]
-            n = xn_b[k]
-            angle = m * theta - n * phi
-            bmod += bmnc_b[s_idx, k] * np.cos(angle)
-        return bmod
-    else:
-        raise ValueError("VMEC file does not contain Boozer spectrum (bmnc_b). "
-                         "Run booz_xform first or use a VMEC file with Boozer data.")
-
-
-def compute_aphi_profile(vmec, rho_grid):
-    """Compute A_phi profile from VMEC flux data.
-
-    A_phi = -torflux * integral(iota ds) from VMEC data.
-    In SIMPLE convention: A_phi is the toroidal vector potential.
-    """
-    ns = vmec["ns"]
-    s_half = np.linspace(0, 1, ns)
-
-    # VMEC stores phi (toroidal flux), from which we derive A_phi
-    # In the Boozer convention used by SIMPLE:
-    # A_theta = torflux * s (by definition)
-    # A_phi comes from the VMEC sA_phi array
-    # For now, approximate: A_phi = -torflux * integral(iota, ds)
-    torflux = abs(vmec["phi_edge"]) / (2.0 * np.pi)
-    iota = vmec["iotaf"] if vmec["iotaf"] is not None else vmec["iotas"]
-
-    # Integrate iota over s to get A_phi
-    A_phi = np.zeros(ns)
-    for i in range(1, ns):
-        ds = s_half[i] - s_half[i - 1]
-        A_phi[i] = A_phi[i - 1] - torflux * iota[i] * ds
-
-    # Interpolate onto rho grid
-    s_grid = rho_grid**2
-    A_phi_on_rho = np.interp(s_grid, s_half, A_phi)
-
-    return A_phi_on_rho, torflux
-
-
-def generate_boozer_chartmap(wout_file, output_file, n_rho=62, n_theta=63, n_phi=64):
-    """Generate extended Boozer chartmap from VMEC wout.nc."""
-    vmec = read_vmec_wout(wout_file)
-    nfp = vmec["nfp"]
-    ns_vmec = vmec["ns"]
-
-    print(f"VMEC file: {wout_file}")
-    print(f"  nfp={nfp}, ns={ns_vmec}")
-    print(f"Output grid: {n_rho} x {n_theta} x {n_phi}")
-
-    twopi = 2.0 * np.pi
-
-    # Build grids
-    rho_grid = np.linspace(0.0, 1.0, n_rho)
-    rho_grid[0] = 1e-6  # avoid axis singularity
-    theta_grid = np.linspace(0.0, twopi, n_theta, endpoint=False)
-    zeta_grid = np.linspace(0.0, twopi / nfp, n_phi, endpoint=False)
-
-    # Allocate output arrays
-    X = np.zeros((n_rho, n_theta, n_phi))
-    Y = np.zeros((n_rho, n_theta, n_phi))
-    Z = np.zeros((n_rho, n_theta, n_phi))
-    Bmod = np.zeros((n_rho, n_theta, n_phi))
-    B_theta_arr = np.zeros(n_rho)
-    B_phi_arr = np.zeros(n_rho)
-
-    # Interpolation from half-grid to rho grid
-    s_half = np.linspace(0, 1, ns_vmec)
-
-    # Get B_theta (buco) and B_phi (bvco) profiles
-    if vmec["buco"] is not None:
-        B_theta_arr = np.interp(rho_grid**2, s_half, vmec["buco"])
-    if vmec["bvco"] is not None:
-        B_phi_arr = np.interp(rho_grid**2, s_half, vmec["bvco"])
-
-    # Compute A_phi profile
-    A_phi_arr, torflux = compute_aphi_profile(vmec, rho_grid)
-
-    # For each flux surface, compute Boozer angle mapping and evaluate geometry
-    for ir in range(n_rho):
-        s = rho_grid[ir] ** 2
-        # Find nearest VMEC half-grid index
-        s_idx = min(int(s * (ns_vmec - 1) + 0.5), ns_vmec - 1)
-        s_idx = max(s_idx, 1)  # avoid axis
-
-        # For simplicity, use VMEC angles directly as approximate Boozer angles.
-        # A proper implementation would use booz_xform for the full transformation.
-        # The geometry (X,Y,Z) is evaluated at VMEC angles and the Bmod
-        # is evaluated using the Boozer spectrum if available.
-        for it in range(n_theta):
-            for ip in range(n_phi):
-                theta_v = theta_grid[it]
-                phi_v = zeta_grid[ip]
-
-                R, Zval = evaluate_vmec_position(vmec, s_idx, theta_v, phi_v)
-                X[ir, it, ip] = R * np.cos(phi_v)
-                Y[ir, it, ip] = R * np.sin(phi_v)
-                Z[ir, it, ip] = Zval
-
-        # Evaluate Bmod on the Boozer grid
-        if vmec["has_boozer"]:
-            theta_2d, zeta_2d = np.meshgrid(theta_grid, zeta_grid, indexing="ij")
-            Bmod[ir, :, :] = evaluate_vmec_bmod(vmec, s_idx, theta_2d, zeta_2d)
-        else:
-            # Fallback: use |B| from VMEC (not truly in Boozer coords)
-            Bmod[ir, :, :] = abs(B_phi_arr[ir]) * 0.01  # rough estimate
-
-    # Compute rho_lcfs (last closed flux surface in rho_tor)
-    rho_lcfs = rho_grid[-1]
-
-    # Write extended chartmap NetCDF
-    print(f"Writing {output_file}")
-    ds = netCDF4.Dataset(output_file, "w", format="NETCDF4")
-
-    # Dimensions
-    ds.createDimension("rho", n_rho)
-    ds.createDimension("theta", n_theta)
-    ds.createDimension("zeta", n_phi)
-
-    # Coordinate variables
-    v = ds.createVariable("rho", "f8", ("rho",))
-    v[:] = rho_grid
-
-    v = ds.createVariable("theta", "f8", ("theta",))
-    v[:] = theta_grid
-
-    v = ds.createVariable("zeta", "f8", ("zeta",))
-    v[:] = zeta_grid
-
-    # Geometry (stored as zeta, theta, rho for chartmap convention)
-    for name, arr in [("x", X), ("y", Y), ("z", Z)]:
-        v = ds.createVariable(name, "f8", ("zeta", "theta", "rho"))
-        # Transpose from (rho, theta, zeta) to (zeta, theta, rho)
-        v[:] = np.transpose(arr, (2, 1, 0))
-        v.units = "cm"
-
-    # Boozer field data - 1D profiles
-    v = ds.createVariable("A_phi", "f8", ("rho",))
-    v[:] = A_phi_arr
-
-    v = ds.createVariable("B_theta", "f8", ("rho",))
-    v[:] = B_theta_arr
-
-    v = ds.createVariable("B_phi", "f8", ("rho",))
-    v[:] = B_phi_arr
-
-    # Boozer field data - 3D
-    v = ds.createVariable("Bmod", "f8", ("zeta", "theta", "rho"))
-    v[:] = np.transpose(Bmod, (2, 1, 0))
-
-    # num_field_periods must be a scalar variable (not attribute) for libneo
-    v = ds.createVariable("num_field_periods", "i4")
-    v[:] = np.int32(nfp)
-
-    # Global attributes
-    ds.rho_convention = "rho_tor"
-    ds.zeta_convention = "cyl"
-    ds.rho_lcfs = rho_lcfs
-    ds.boozer_field = np.int32(1)
-    ds.torflux = torflux
-
-    ds.close()
-    print(f"Done. torflux={torflux:.6e}")
-
-
-def main():
-    parser = argparse.ArgumentParser(
-        description="Generate extended Boozer chartmap from VMEC wout.nc"
-    )
-    parser.add_argument("--wout", required=True, help="Input VMEC wout.nc file")
-    parser.add_argument("--out", required=True, help="Output Boozer chartmap .nc file")
-    parser.add_argument("--nrho", type=int, default=62, help="Number of radial points")
-    parser.add_argument("--ntheta", type=int, default=63, help="Number of poloidal points")
-    parser.add_argument("--nphi", type=int, default=64, help="Number of toroidal points")
-    args = parser.parse_args()
-
-    generate_boozer_chartmap(args.wout, args.out, args.nrho, args.ntheta, args.nphi)
-
-
-if __name__ == "__main__":
-    main()

From a5887c63ff170eb801270a49395309001ccab6cf Mon Sep 17 00:00:00 2001
From: Christopher Albert <albert@tugraz.at>
Date: Sat, 28 Mar 2026 14:22:05 +0100
Subject: [PATCH 8/9] fix: make Boozer chartmap generator executable

---
 test/tests/generate_test_boozer_chartmap.py | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 mode change 100644 => 100755 test/tests/generate_test_boozer_chartmap.py

diff --git a/test/tests/generate_test_boozer_chartmap.py b/test/tests/generate_test_boozer_chartmap.py
old mode 100644
new mode 100755

From f133009deac700d745f0cca75212d65fa9b40193 Mon Sep 17 00:00:00 2001
From: Christopher Albert <albert@tugraz.at>
Date: Sat, 28 Mar 2026 14:33:40 +0100
Subject: [PATCH 9/9] fix: remove Python netcdf dependency from chartmap
 fixture

---
 test/tests/CMakeLists.txt                   |   8 +-
 test/tests/generate_test_boozer_chartmap.c  | 164 ++++++++++++++++++++
 test/tests/generate_test_boozer_chartmap.py | 106 -------------
 3 files changed, 169 insertions(+), 109 deletions(-)
 create mode 100644 test/tests/generate_test_boozer_chartmap.c
 delete mode 100755 test/tests/generate_test_boozer_chartmap.py

diff --git a/test/tests/CMakeLists.txt b/test/tests/CMakeLists.txt
index dc617c98..46af006e 100644
--- a/test/tests/CMakeLists.txt
+++ b/test/tests/CMakeLists.txt
@@ -331,12 +331,14 @@ add_test(NAME test_array_utils COMMAND test_array_utils.x)
     add_test(NAME test_chartmap_metadata COMMAND test_chartmap_metadata.x)
     set_tests_properties(test_chartmap_metadata PROPERTIES LABELS "unit")
 
-    # Boozer chartmap test: generate test file with Python, then run Fortran test
+    # Boozer chartmap test: generate test file with a compiled NetCDF tool, then run Fortran test
+    add_executable(generate_test_boozer_chartmap.x generate_test_boozer_chartmap.c)
+    target_link_libraries(generate_test_boozer_chartmap.x netcdf m)
     add_custom_command(
         OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/test_boozer_chartmap.nc
-        COMMAND ${Python_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/generate_test_boozer_chartmap.py
+        COMMAND $<TARGET_FILE:generate_test_boozer_chartmap.x> ${CMAKE_CURRENT_BINARY_DIR}/test_boozer_chartmap.nc
         WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
-        DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/generate_test_boozer_chartmap.py
+        DEPENDS generate_test_boozer_chartmap.x
         COMMENT "Generating test Boozer chartmap NetCDF"
     )
     add_custom_target(generate_test_boozer_chartmap
diff --git a/test/tests/generate_test_boozer_chartmap.c b/test/tests/generate_test_boozer_chartmap.c
new file mode 100644
index 00000000..34f8e9ef
--- /dev/null
+++ b/test/tests/generate_test_boozer_chartmap.c
@@ -0,0 +1,164 @@
+#include <math.h>
+#include <netcdf.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#define CHECK_NC(call, context)                                                   \
+    do {                                                                          \
+        int status_ = (call);                                                     \
+        if (status_ != NC_NOERR) {                                                \
+            fprintf(stderr, "NetCDF error in %s: %s\n", context,                  \
+                    nc_strerror(status_));                                        \
+            return 1;                                                             \
+        }                                                                         \
+    } while (0)
+
+int main(int argc, char **argv)
+{
+    enum
+    {
+        nfp = 2,
+        n_rho = 17,
+        n_theta = 33,
+        n_zeta = 33
+    };
+
+    const double b0 = 2.0e4;
+    const double r0 = 150.0;
+    const double a = 50.0;
+    const double iota = 0.5;
+    const double twopi = 2.0 * acos(-1.0);
+    const double epsilon = a / r0;
+    const double torflux = 0.5 * b0 * a * a;
+    const char *filename = argc > 1 ? argv[1] : "test_boozer_chartmap.nc";
+
+    double rho[n_rho];
+    double theta[n_theta];
+    double zeta[n_zeta];
+    double a_phi[n_rho];
+    double b_theta[n_rho];
+    double b_phi[n_rho];
+    double *x = NULL;
+    double *y = NULL;
+    double *z = NULL;
+    double *bmod = NULL;
+    int ncid;
+    int dim_rho, dim_theta, dim_zeta;
+    int dims_3d[3];
+    int var_rho, var_theta, var_zeta, var_x, var_y, var_z, var_aphi, var_btheta;
+    int var_bphi, var_bmod, var_nfp;
+
+    x = malloc((size_t)n_rho * n_theta * n_zeta * sizeof(*x));
+    y = malloc((size_t)n_rho * n_theta * n_zeta * sizeof(*y));
+    z = malloc((size_t)n_rho * n_theta * n_zeta * sizeof(*z));
+    bmod = malloc((size_t)n_rho * n_theta * n_zeta * sizeof(*bmod));
+    if (x == NULL || y == NULL || z == NULL || bmod == NULL) {
+        fprintf(stderr, "Allocation failure\n");
+        free(x);
+        free(y);
+        free(z);
+        free(bmod);
+        return 1;
+    }
+
+    for (int ir = 0; ir < n_rho; ++ir) {
+        rho[ir] = (double)ir / (double)(n_rho - 1);
+        a_phi[ir] = -torflux * iota * rho[ir] * rho[ir];
+        b_theta[ir] = b0 * epsilon * iota;
+        b_phi[ir] = b0 * r0;
+    }
+    rho[0] = 1.0e-6;
+    a_phi[0] = -torflux * iota * rho[0] * rho[0];
+
+    for (int it = 0; it < n_theta; ++it) {
+        theta[it] = twopi * (double)it / (double)n_theta;
+    }
+
+    for (int iz = 0; iz < n_zeta; ++iz) {
+        zeta[iz] = (twopi / (double)nfp) * (double)iz / (double)n_zeta;
+    }
+
+    for (int iz = 0; iz < n_zeta; ++iz) {
+        const double zeta_val = zeta[iz];
+        for (int it = 0; it < n_theta; ++it) {
+            const double theta_val = theta[it];
+            const double cos_theta = cos(theta_val);
+            const double sin_theta = sin(theta_val);
+            for (int ir = 0; ir < n_rho; ++ir) {
+                const double rho_val = rho[ir];
+                const double r_minor = a * rho_val;
+                const double radius = r0 + r_minor * cos_theta;
+                const size_t idx = ((size_t)iz * n_theta + (size_t)it) * n_rho +
+                                   (size_t)ir;
+                x[idx] = radius * cos(zeta_val);
+                y[idx] = radius * sin(zeta_val);
+                z[idx] = r_minor * sin_theta;
+                bmod[idx] = b0 / (1.0 + rho_val * epsilon * cos_theta);
+            }
+        }
+    }
+
+    CHECK_NC(nc_create(filename, NC_NETCDF4, &ncid), "create file");
+    CHECK_NC(nc_def_dim(ncid, "rho", n_rho, &dim_rho), "def dim rho");
+    CHECK_NC(nc_def_dim(ncid, "theta", n_theta, &dim_theta), "def dim theta");
+    CHECK_NC(nc_def_dim(ncid, "zeta", n_zeta, &dim_zeta), "def dim zeta");
+
+    dims_3d[0] = dim_zeta;
+    dims_3d[1] = dim_theta;
+    dims_3d[2] = dim_rho;
+
+    CHECK_NC(nc_def_var(ncid, "rho", NC_DOUBLE, 1, &dim_rho, &var_rho), "def rho");
+    CHECK_NC(nc_def_var(ncid, "theta", NC_DOUBLE, 1, &dim_theta, &var_theta), "def theta");
+    CHECK_NC(nc_def_var(ncid, "zeta", NC_DOUBLE, 1, &dim_zeta, &var_zeta), "def zeta");
+    CHECK_NC(nc_def_var(ncid, "x", NC_DOUBLE, 3, dims_3d, &var_x), "def x");
+    CHECK_NC(nc_def_var(ncid, "y", NC_DOUBLE, 3, dims_3d, &var_y), "def y");
+    CHECK_NC(nc_def_var(ncid, "z", NC_DOUBLE, 3, dims_3d, &var_z), "def z");
+    CHECK_NC(nc_def_var(ncid, "A_phi", NC_DOUBLE, 1, &dim_rho, &var_aphi), "def A_phi");
+    CHECK_NC(nc_def_var(ncid, "B_theta", NC_DOUBLE, 1, &dim_rho, &var_btheta), "def B_theta");
+    CHECK_NC(nc_def_var(ncid, "B_phi", NC_DOUBLE, 1, &dim_rho, &var_bphi), "def B_phi");
+    CHECK_NC(nc_def_var(ncid, "Bmod", NC_DOUBLE, 3, dims_3d, &var_bmod), "def Bmod");
+    CHECK_NC(nc_def_var(ncid, "num_field_periods", NC_INT, 0, NULL, &var_nfp),
+             "def num_field_periods");
+
+    CHECK_NC(nc_put_att_text(ncid, var_x, "units", 2, "cm"), "put x units");
+    CHECK_NC(nc_put_att_text(ncid, var_y, "units", 2, "cm"), "put y units");
+    CHECK_NC(nc_put_att_text(ncid, var_z, "units", 2, "cm"), "put z units");
+    CHECK_NC(nc_put_att_text(ncid, NC_GLOBAL, "rho_convention", 7, "rho_tor"),
+             "put rho_convention");
+    CHECK_NC(nc_put_att_text(ncid, NC_GLOBAL, "zeta_convention", 6, "boozer"),
+             "put zeta_convention");
+    CHECK_NC(nc_put_att_double(ncid, NC_GLOBAL, "rho_lcfs", NC_DOUBLE, 1,
+                               (const double[]){1.0}),
+             "put rho_lcfs");
+    CHECK_NC(nc_put_att_int(ncid, NC_GLOBAL, "boozer_field", NC_INT, 1,
+                            (const int[]){1}),
+             "put boozer_field");
+    CHECK_NC(nc_put_att_double(ncid, NC_GLOBAL, "torflux", NC_DOUBLE, 1, &torflux),
+             "put torflux");
+
+    CHECK_NC(nc_enddef(ncid), "enddef");
+
+    CHECK_NC(nc_put_var_double(ncid, var_rho, rho), "put rho");
+    CHECK_NC(nc_put_var_double(ncid, var_theta, theta), "put theta");
+    CHECK_NC(nc_put_var_double(ncid, var_zeta, zeta), "put zeta");
+    CHECK_NC(nc_put_var_double(ncid, var_x, x), "put x");
+    CHECK_NC(nc_put_var_double(ncid, var_y, y), "put y");
+    CHECK_NC(nc_put_var_double(ncid, var_z, z), "put z");
+    CHECK_NC(nc_put_var_double(ncid, var_aphi, a_phi), "put A_phi");
+    CHECK_NC(nc_put_var_double(ncid, var_btheta, b_theta), "put B_theta");
+    CHECK_NC(nc_put_var_double(ncid, var_bphi, b_phi), "put B_phi");
+    CHECK_NC(nc_put_var_double(ncid, var_bmod, bmod), "put Bmod");
+    CHECK_NC(nc_put_var_int(ncid, var_nfp, (const int[]){nfp}), "put num_field_periods");
+    CHECK_NC(nc_close(ncid), "close");
+
+    printf("Generated %s\n", filename);
+    printf("  torflux=%14.6e, B0=%8.1f, R0=%8.1f\n", torflux, b0, r0);
+    printf("  nfp=%d, nrho=%d, ntheta=%d, nzeta=%d\n", nfp, n_rho, n_theta, n_zeta);
+
+    free(x);
+    free(y);
+    free(z);
+    free(bmod);
+    return 0;
+}
diff --git a/test/tests/generate_test_boozer_chartmap.py b/test/tests/generate_test_boozer_chartmap.py
deleted file mode 100755
index b2b80fe3..00000000
--- a/test/tests/generate_test_boozer_chartmap.py
+++ /dev/null
@@ -1,106 +0,0 @@
-#!/usr/bin/env python3
-"""Generate a minimal Boozer chartmap NetCDF test file with known analytic values."""
-
-import numpy as np
-import netCDF4
-
-def generate_test_file(filename):
-    """Create a test Boozer chartmap with simple analytic field.
-
-    Uses a tokamak-like field: Bmod ~ B0 * (1 - epsilon * cos(theta))
-    with B_theta ~ constant, B_phi ~ constant, A_phi ~ -torflux * iota * s.
-    """
-    nfp = 2
-    n_rho = 17
-    n_theta = 33
-    n_phi = 33
-    B0 = 2.0e4  # 2 T in Gauss
-    R0 = 150.0  # 1.5 m in cm
-    a = 50.0    # 0.5 m in cm
-    epsilon = a / R0
-    iota = 0.5
-    torflux = B0 * np.pi * a**2 / (2.0 * np.pi)  # approximate
-
-    twopi = 2.0 * np.pi
-
-    rho_grid = np.linspace(0.0, 1.0, n_rho)
-    rho_grid[0] = 1e-6
-    theta_grid = np.linspace(0.0, twopi, n_theta, endpoint=False)
-    zeta_grid = np.linspace(0.0, twopi / nfp, n_phi, endpoint=False)
-
-    # Geometry: concentric circular cross-section tokamak
-    X = np.zeros((n_rho, n_theta, n_phi))
-    Y = np.zeros((n_rho, n_theta, n_phi))
-    Z = np.zeros((n_rho, n_theta, n_phi))
-    Bmod = np.zeros((n_rho, n_theta, n_phi))
-
-    for ir in range(n_rho):
-        rho = rho_grid[ir]
-        r = a * rho  # minor radius
-        for it in range(n_theta):
-            theta = theta_grid[it]
-            R = R0 + r * np.cos(theta)
-            Zval = r * np.sin(theta)
-            for ip in range(n_phi):
-                phi = zeta_grid[ip]
-                X[ir, it, ip] = R * np.cos(phi)
-                Y[ir, it, ip] = R * np.sin(phi)
-                Z[ir, it, ip] = Zval
-                Bmod[ir, it, ip] = B0 / (1.0 + rho * epsilon * np.cos(theta))
-
-    # Profiles
-    A_phi = -torflux * iota * rho_grid**2
-    B_theta = np.full(n_rho, B0 * epsilon * iota)
-    B_phi = np.full(n_rho, B0 * R0)
-
-    # Write
-    ds = netCDF4.Dataset(filename, "w", format="NETCDF4")
-
-    ds.createDimension("rho", n_rho)
-    ds.createDimension("theta", n_theta)
-    ds.createDimension("zeta", n_phi)
-
-    v = ds.createVariable("rho", "f8", ("rho",))
-    v[:] = rho_grid
-
-    v = ds.createVariable("theta", "f8", ("theta",))
-    v[:] = theta_grid
-
-    v = ds.createVariable("zeta", "f8", ("zeta",))
-    v[:] = zeta_grid
-
-    for name, arr in [("x", X), ("y", Y), ("z", Z)]:
-        v = ds.createVariable(name, "f8", ("zeta", "theta", "rho"))
-        v[:] = np.transpose(arr, (2, 1, 0))
-        v.units = "cm"
-
-    v = ds.createVariable("A_phi", "f8", ("rho",))
-    v[:] = A_phi
-
-    v = ds.createVariable("B_theta", "f8", ("rho",))
-    v[:] = B_theta
-
-    v = ds.createVariable("B_phi", "f8", ("rho",))
-    v[:] = B_phi
-
-    v = ds.createVariable("Bmod", "f8", ("zeta", "theta", "rho"))
-    v[:] = np.transpose(Bmod, (2, 1, 0))
-
-    # num_field_periods must be a scalar variable (not attribute) for libneo
-    v = ds.createVariable("num_field_periods", "i4")
-    v[:] = np.int32(nfp)
-
-    ds.rho_convention = "rho_tor"
-    ds.zeta_convention = "cyl"
-    ds.rho_lcfs = 1.0
-    ds.boozer_field = np.int32(1)
-    ds.torflux = torflux
-
-    ds.close()
-    print(f"Generated {filename}")
-    print(f"  torflux={torflux:.6e}, B0={B0:.1f}, R0={R0:.1f}, a={a:.1f}")
-    print(f"  nfp={nfp}, nrho={n_rho}, ntheta={n_theta}, nphi={n_phi}")
-
-
-if __name__ == "__main__":
-    generate_test_file("test_boozer_chartmap.nc")