spectre.converters.vmec_helpers
===============================

.. py:module:: spectre.converters.vmec_helpers

.. autoapi-nested-parse::

   VMEC physics helpers and generic VMEC-to-SPECTRE builder.

   Contains everything that is shared across VMEC adapters:

   - Physical constants and polynomial utilities (:func:`power_series`)
   - Profile discretisation functions (pressure, iota, volume/sheet current)
   - Boundary Fourier-coefficient extraction from a simsopt surface object
   - :class:`ContinuousPhysics` – the common data structure that every adapter
     must produce
   - :func:`build_spectre_input` – turns a :class:`ContinuousPhysics` into a validated
     SPECTRE :class:`~spectre.file_io.input_parameters.InputParameters` object

   Adding a new VMEC interface
   ---------------------------
   Write a new adapter module that imports :class:`ContinuousPhysics` from here and
   implements its own extraction function::

       def extract_<name>_physics(source) -> ContinuousPhysics:
           ...
           return ContinuousPhysics(
               curtor=..., phiedge=..., nfp=...,
               rbc={...}, zbs={...},
               profile_type='iota',   # or 'current'
               pressure=...,          # callable pressure(s), mu0-scaled
               iota=...,              # callable iota(s)
               current_density=...,   # callable j(s) = dI/ds
           )

   Then call :func:`build_spectre_input` directly — no other file in this package
   needs to change.









Module Contents
---------------

.. py:data:: MU0

.. py:function:: power_series(coeffs)

   Return a callable that evaluates ``sum_i coeffs[i] * s**i``.

   VMEC stores pressure, iota, and current profiles as power-series
   coefficients in normalised toroidal flux s ∈ [0, 1].

   Args:
       coeffs: sequence of polynomial coefficients ordered from degree 0.

   Returns:
       Callable ``f(s)`` that accepts scalars or arrays.


.. py:function:: pressure_function(am, pres_scale)

   Return a callable ``pressure(s)`` for VMEC pressure coefficients.

   The returned callable evaluates the ``am`` power series scaled by
   ``pres_scale`` and mu0, so the result is in SPECTRE's internal pressure
   units.

   Args:
       am: VMEC pressure polynomial coefficients (array-like).
       pres_scale: VMEC pressure scale factor (multiplies the polynomial).

   Returns:
       Callable ``pressure(s)`` accepting scalars or arrays of normalised
       toroidal flux s ∈ [0, 1].


.. py:function:: compute_pressure_profile(pressure, tfluxes)

   Volume-averaged pressure in each sub-volume between consecutive interfaces.

   SPECTRE stores a single pressure value per volume.  We obtain it by
   integrating the continuous profile ``pressure(s)`` over the s-interval
   [s_{i-1}, s_i] of each volume and dividing by the interval width, giving
   the volume average.

   Args:
       pressure: callable ``pressure(s)`` returning the pressure in SPECTRE
           units (i.e. already mu0-scaled) at normalised toroidal flux s.
       tfluxes: 1-D array of interface toroidal-flux values in (0, 1].

   Returns:
       1-D array of length ``len(tfluxes)`` with <pressure> per volume.


.. py:function:: iota_function(vmec)

   Return a callable ``iota(s)`` for a simsopt ``Vmec`` object.

   Iota is the one profile whose *representation* differs by backend, so the
   dispatch lives here:

   - **indata source** (loaded from ``input.*``): evaluate the ``ai``
     power-series coefficients via :func:`power_series`, matching VMEC's own
     profile evaluation.
   - **wout source** (loaded from ``wout_*.nc``): cubic interpolation of the
     converged full-mesh ``iotaf`` grid.  The ``ai`` polynomial fit may not
     reproduce ``iotaf`` exactly, so the output grid is preferred when present.

   Args:
       vmec: simsopt ``Vmec`` instance with either an ``indata`` (input file)
           or ``wout`` (netCDF output) attribute.

   Returns:
       Callable ``iota(s)`` accepting scalars or arrays of normalised toroidal
       flux s ∈ [0, 1].

   Raises:
       RuntimeError: if ``vmec`` has neither ``indata`` nor ``wout``.


.. py:function:: normalized_enclosed_current(current_density, s)

   Normalised cumulative enclosed toroidal current Î(s) at flux value(s) ``s``.

   The current density j(s) = dI/ds is supplied as a callable.  The enclosed
   current is its integral I(s) = ∫₀ˢ j ds', and VMEC normalises so that I(1)
   corresponds to ``curtor``.  This function returns the dimensionless,
   normalised shape

       Î(s) = (∫₀ˢ j ds') / (∫₀¹ j ds')

   which runs from 0 at the axis to 1 at the plasma boundary.  Dividing by the
   total ∫₀¹ j is essential: the density is only a *shape*, and its overall
   scale is arbitrary (it is commonly — but not always — normalised so that
   ∫₀¹ j = 1).  Omitting this division is what made the old volume-current path
   wrong by the factor ∫₀¹ j.

   Args:
       current_density: callable ``j(s) = dI/ds`` (e.g. from ``power_series``
           of the VMEC ``ac`` coefficients).
       s: scalar or 1-D array of normalised toroidal-flux values.

   Returns:
       Î(s) with the same shape as ``s`` (scalar in → scalar out).  Returns all
       zeros when ∫₀¹ j ds == 0 (e.g. a zero current profile).


.. py:function:: distribute_current_to_interfaces(current_density, curtor, tfluxes, inface_current=False)

   Map a continuous VMEC current profile onto SPECTRE interfaces.

   .. warning::

       **Provisional / not final.**  How to map a *continuous* (VMEC) current
       profile onto the *discrete* interfaces of an MRxMHD (SPEC/SPECTRE)
       equilibrium is a fundamental open question.  The control-volume scheme
       below is one physically-motivated choice; expect it to change.

   inface_current = False (default) -> all current carried in the volumes;
   there is no *net* toroidal current on the interfaces.
   inface_current = True -> all current carried on the interfaces (``isurf``).

   - **Volume current** (``inface_current=False``) ``ivolume[i] = (−curtor·μ₀) · Î(sᵢ)``
     — the *cumulative* enclosed current at interface i (Î from
     :func:`normalized_enclosed_current`); ``isurf`` is then zero.
   - **Sheet current** (``inface_current=True``) ``isurf[i] = (−curtor·μ₀) · (Î(eᵢ) − Î(eᵢ₋₁))``
     — the current density integrated over the *control volume* owned by
     interface i, whose edges sit at the midpoints to its neighbours::

         e₀ = 0                     # magnetic axis
         eᵢ = (sᵢ + sᵢ₊₁)/2         # midpoint to the next interface
         e_N = s_N                  # outer edge of the last cell = last interface

     So each interface collects current from halfway to its inner neighbour to
     halfway to its outer neighbour — an *average of inside and outside*, not
     just the current of the volume just inside it.  For equidistant interfaces
     this reduces to the spot-evaluation j(sᵢ)·Δ; for non-equidistant
     interfaces a sparsely-spaced interface automatically owns a wider band and
     carries more of the surrounding current.

   Current is conserved: the control-volume increments telescope to Î(s_N), so
   ``ivolume[-1] + sum(isurf) = (−curtor·μ₀)·Î(s_N)`` in either mode.  When the last
   interface is the true plasma boundary (s_N = 1) this is the full −curtor·μ₀;
   when tfluxes[-1] < 1 the total is reduced by the appropriate amount, because
   the current flowing beyond the (interior) SPECTRE boundary is not enclosed.

   The last interface carries ``K_N = (−curtor·μ₀)(Î(s_N) − Î(e_{N-1}))`` — the
   current from the last midpoint out to the boundary.  In a fixed-boundary run
   SPECTRE ignores ``Isurf(Mvol)``, but this value is preserved so that
   :func:`~spectre.converters.fix2free.fix2free` carries it over: when the
   plasma boundary is promoted to an interior interface in the free-boundary
   equilibrium, the slot becomes physical.

   Args:
       current_density: callable ``j(s) = dI/ds`` (e.g. from ``power_series``
           of the VMEC ``ac`` coefficients).
       curtor: total toroidal current from VMEC [A].
       tfluxes: 1-D array of interface toroidal-flux values in (0, 1].
       inface_current: if ``False`` (default), all current is distributed as
           volume current (``ivolume``); if ``True``, all current is placed on
           the interfaces as sheet current (``isurf``).

   Returns:
       Tuple ``(ivolume, isurf)`` of 1-D arrays of length ``len(tfluxes)``,
       both in μ₀·A.


.. py:function:: vmec_boundary_to_fourier_coeffs(boundary, angle_fix=True)

   Extract RBC and ZBS Fourier coefficients from a simsopt boundary surface.

   A simsopt ``SurfaceRZFourier`` stores the full (m, n) mode spectrum
   including negative-m modes.  We iterate over the positive-m half
   (num_modes = len(boundary.m) // 2 + 1) and read each cosine (R) and
   sine (Z) coefficient via the surface's accessor methods.

   Args:
       boundary: simsopt ``SurfaceRZFourier`` object (e.g. ``vmec.boundary``).
       angle_fix: if ``True`` (default), apply the θ→-θ handedness correction
           for m>0 modes: the toroidal index is negated (n→-n) and the ZBS
           coefficient is negated.  This pre-applies what ``lchangeangle=True``
           would do at Fortran runtime, so the output TOML needs no
           ``lchangeangle`` flag.  Set to ``False`` to keep the raw VMEC
           convention.

   Returns:
       Tuple ``(rbc, zbs)`` where each is a dict keyed by ``'(m, n)'`` strings
       mapping to the corresponding Fourier coefficient.


.. py:function:: compute_interface_guesses(wout, tfluxes, mpol, ntor)

   Interpolate VMEC wout surfaces at each interface position.

   Reads ``rmnc``/``zmns`` from the simsopt wout object, interpolates
   cubically on the VMEC radial grid, and evaluates at each value in
   ``tfluxes``.  The θ→-θ angle fix is applied (n→-n for m>0, ZBS negated)
   so the returned coefficients are in SPECTRE's right-handed coordinate
   convention — identical to what ``lchangeangle=True`` would produce at
   Fortran runtime.

   Only handles stellarator-symmetric equilibria (``lasym=False``).

   Args:
       wout: simsopt wout object (``vmec_obj.wout``).
       tfluxes: 1-D array of interface toroidal-flux values in (0, 1].
       mpol: maximum poloidal mode number to include in the output.
       ntor: maximum toroidal mode number (reduced, without nfp) to include.

   Returns:
       dict suitable for ``PhysicsParameters.allrzrz``:
       ``{'interface_1': {'rbc': {...}, 'zbs': {...}, 'rbs': {}, 'zbc': {}},
       ...}``


.. py:class:: ContinuousPhysics

   SPECTRE-relevant physics extracted from a VMEC equilibrium.

   This is the contract between a VMEC adapter and :func:`build_spectre_input`.
   Every adapter must populate all fields; the builder has no knowledge of
   where the data came from.

   The profiles are stored as **callables** ``f(s)`` on normalised toroidal
   flux s ∈ [0, 1], not as backend-specific coefficient arrays.  This is what
   keeps the builder backend-agnostic: each adapter decides how to evaluate a
   profile (power series, grid interpolation, spline, …) and hands over a plain
   function.  Being a mutable dataclass, a caller can override a profile in
   place, e.g. ``physics.iota = lambda s: 0.9 + 0.1 * s``.

   Attributes:
       curtor: total toroidal current [A].
       phiedge: toroidal magnetic flux at the plasma boundary [Wb].
       nfp: number of toroidal field periods.
       rbc: cosine boundary coefficients as ``{'(m, n)': value}`` dict.
       zbs: sine boundary coefficients as ``{'(m, n)': value}`` dict.
       profile_type: ``'iota'`` for a rotational-transform constraint,
           ``'current'`` for an integrated-current constraint.
       pressure: callable ``pressure(s)`` returning the pressure in SPECTRE
           units (mu0-scaled); see :func:`pressure_function`.
       iota: callable ``iota(s)`` returning the rotational transform; see
           :func:`iota_function`.  Only used when ``profile_type == 'iota'``.
       current_density: callable ``j(s) = dI/ds`` (current-profile shape; its
           overall scale is irrelevant, it is normalised internally).  Only
           used when ``profile_type == 'current'``.


   .. py:attribute:: curtor
      :type:  float


   .. py:attribute:: phiedge
      :type:  float


   .. py:attribute:: nfp
      :type:  int


   .. py:attribute:: rbc
      :type:  dict


   .. py:attribute:: zbs
      :type:  dict


   .. py:attribute:: profile_type
      :type:  str


   .. py:attribute:: pressure
      :type:  collections.abc.Callable


   .. py:attribute:: iota
      :type:  collections.abc.Callable


   .. py:attribute:: current_density
      :type:  collections.abc.Callable


.. py:function:: build_spectre_input(physics, tfluxes, mpol, ntor, inface_current=False, interface_guesses=None)

   Assemble a SPECTRE :class:`InputParameters` object from VMEC physics.

   This function has no dependency on any VMEC interface library; it operates
   purely on the :class:`ContinuousPhysics` contract (profile callables + geometry).

   Args:
       physics: :class:`ContinuousPhysics` dataclass instance.
       tfluxes: 1-D array of interface toroidal-flux values in (0, 1].
       mpol: poloidal mode resolution for SPECTRE.
       ntor: toroidal mode resolution for SPECTRE.
       inface_current: if ``False`` (default), the toroidal current is
           distributed as volume current (``ivolume``); if ``True``, it is
           placed entirely on the interfaces as sheet current (``isurf``).
           Only used when ``physics.profile_type == 'current'``.  See
           :func:`distribute_current_to_interfaces` for the (provisional)
           discretisation scheme.
       interface_guesses: optional dict in ``allrzrz`` format returned by
           :func:`compute_interface_guesses`.  When provided, it is stored as
           the initial geometry guess and ``linitialize`` is set to 0 (use
           supplied surfaces).  When ``None``, ``linitialize=1`` causes SPECTRE
           to linearly interpolate an initial guess from the boundary.

   Returns:
       Tuple ``(dict_pydantic, profiles)`` where ``profiles`` is a dict of
       intermediate arrays (``'pressure'``, ``'iota'``, ``'ivolume'``,
       ``'isurf'``) useful for plotting or inspection.

   Raises:
       RuntimeError: if ``profile_type`` is unrecognised.