Correct q sign-convention docs and enforce positive q (replaces #268)

CLAUDE.md

@@ -257,7 +257,7 @@ decomposition exp(im*theta - in*phi). A comprehensive reference is in
 - **helicity** = ipd * btd (+1 = RH, -1 = LH), computed in `gpec_main` (`gpec/gpec.f`)
 - **nn** (toroidal mode number) is always positive; resonant **m** is always positive
 - **F = R*Bt** is forced positive via ABS() in `read_eq_efit` (`equil/read_eq.f`)
-- **q** is NOT forced positive (read directly from EFIT)
+- **q** is recomputed by field-line integration (`direct_run` in `equil/direct.f`; `inverse_run` in `equil/inverse.f` for inverse formats); the input file's q profile is unused, and the recomputed q is always positive
 - **omega_E** is positive for rotation in the direction of the toroidal coordinate zeta
 - The code does NOT use the COCOS standard
 - For SURFMN interface: `m_surfmn = helicity * m_gpec`

docs/sign_conventions.rst

@@ -96,9 +96,9 @@ Helicity is computed in the ``gpec_main`` program (``gpec/gpec.f``):
 
 - ``ip_direction`` and ``bt_direction`` are set in ``coil.in``.
   "positive" means CCW viewed from above; "negative" means CW.
-- **helicity = +1**: right-handed (RH) --- :math:`B_t` and :math:`I_p` in the
+- **helicity = +1**: right-handed (RH), :math:`B_t` and :math:`I_p` in the
   same direction.
-- **helicity = -1**: left-handed (LH) --- :math:`B_t` and :math:`I_p` opposed.
+- **helicity = -1**: left-handed (LH), :math:`B_t` and :math:`I_p` opposed.
 - The helicity value is stored in the ``gpec_control_output`` netcdf file.
 
 
@@ -121,12 +121,22 @@ Safety Factor :math:`q`
 ========================
 
 - Defined as :math:`q = (\mathbf{B} \cdot \nabla\zeta) / (\mathbf{B} \cdot \nabla\theta)`.
-- **Not forced positive.** Read directly from the EFIT g-file without sign
-  manipulation (unlike :math:`F`).
-- In standard tokamak operation, EFIT provides :math:`q > 0`.
-- Can be adjusted via ``newq0`` in ``equil.in``, which modifies :math:`F`
-  to match while preserving the Grad-Shafranov solution
-  (``direct_run`` in ``equil/direct.f``).
+- **Recomputed, not read.** For direct equilibria (EFIT g-files), the q
+  profile in the file is read but never used. ``direct_run`` in
+  ``equil/direct.f`` rebuilds q on each flux surface from a field-line
+  integration.
+- **Always positive for g-file input.** The integration uses the
+  sign-normalized :math:`\psi` map and :math:`|F|`, so q comes out positive
+  for any combination of signs in the g-file. Reversing all signed
+  quantities in a g-file (:math:`\psi`, ``simag``, ``sibry``, ``cpasma``,
+  q) reproduces the original DCON results bit for bit.
+- **Inverse equilibria recompute q too.** The main inverse path
+  (``inverse_run`` in ``equil/inverse.f``, used by the CHEASE, JSOLVER,
+  TRANSP, and analytic formats) rebuilds q on each surface by a flux-surface
+  integration, exactly as the direct path does; the input file's q profile
+  is read into ``sq_in`` but not used, and the recomputed q is positive. The
+  one exception is the ``chease4`` format (``inverse_chease4_run``), which
+  takes q directly from the input file.
 
 
 Mode Numbers :math:`m` and :math:`n`
@@ -157,7 +167,8 @@ Why Positive :math:`m` Is Always Resonant
 
 Resonant surfaces are found by ``sing_find`` in ``dcon/sing.f``.
 It performs a binary search for flux surfaces where :math:`m = n \cdot q`.
-Since :math:`n > 0` (by convention) and :math:`q > 0` (standard tokamak):
+Since :math:`n > 0` (by convention) and :math:`q > 0` (the recomputed q is
+always positive, see the safety factor section above):
 
 .. math::
 
@@ -284,6 +295,19 @@ system. The conventions described in this document are GPEC-native and predate
 COCOS. Users interfacing with COCOS-aware codes must manually translate between
 conventions.
 
+For g-file input the sign part of the COCOS choice cannot matter:
+``read_eq_efit`` sign-normalizes the :math:`\psi` map, forces
+:math:`F > 0`, and discards the file's q profile, so g-files that differ
+only in sign conventions produce identical results. Field and current
+directions enter solely through ``ip_direction`` and ``bt_direction`` in
+``coil.in``. The main inverse path likewise recomputes q rather than
+inheriting it from the input file. What does matter is the flux
+normalization: GPEC assumes the
+g-eqdsk standard of poloidal flux per radian (Wb/rad). A file carrying the
+full flux (the COCOS 11-18 family) yields a wrong q magnitude. Reports
+that g-files must be supplied in one specific COCOS trace to this unit
+requirement, since the sign choices are normalized away on read.
+
 
 Quick Reference
 ===============
@@ -308,8 +332,8 @@ Quick Reference
      - Always positive
      - Yes, ABS()
    * - :math:`q` (safety factor)
-     - From EFIT, typically positive
-     - No
+     - Recomputed by field-line integration
+     - Always positive
    * - :math:`n` (toroidal mode)
      - Always positive
      - By convention
@@ -330,6 +354,8 @@ Source Code References
 - **Poloidal flux sign**: ``read_eq_efit`` in ``equil/read_eq.f``
 - **F = R*Bt forced positive**: ``read_eq_efit`` in ``equil/read_eq.f``
 - **q definition**: ``dcon/README``
+- **q from field-line integration**: ``direct_run`` in ``equil/direct.f``
+- **q recomputed (inverse path)**: ``inverse_run`` in ``equil/inverse.f``
 - **Helicity computation**: ``gpec_main`` program in ``gpec/gpec.f``
 - **ip/bt direction**: ``input/coil.in``
 - **Poloidal mode range**: ``dcon`` program in ``dcon/dcon.F``

equil/direct.f

@@ -195,7 +195,6 @@ SUBROUTINE direct_run
       sq%name="  sq  "
       sq%title=(/" psi  ","  f   ","  p   ","  q   "/)
       q0=sq%fs(0,4)-sq%fs1(0,4)*sq%xs(0)
-      IF(newq0 == -1)newq0=-q0
 c-----------------------------------------------------------------------
 c     revise q profile.
 c-----------------------------------------------------------------------
@@ -204,7 +203,7 @@ SUBROUTINE direct_run
          f0fac=f0**2*((newq0/q0)**2-1)
          q0=newq0
          DO ipsi=0,mpsi
-            ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2)*SIGN(one,newq0)
+            ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2)
             sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac
             sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac
             rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac

equil/equil.f

@@ -55,6 +55,11 @@ SUBROUTINE equil_read(unit, op_psihigh, op_psilow)
          READ(UNIT=in_unit,NML=equil_output)
       ENDIF
       CALL ascii_close(in_unit)
+c-----------------------------------------------------------------------
+c     guard against negative newq0 (negative q is not supported).
+c-----------------------------------------------------------------------
+      IF(newq0 < 0)CALL program_stop("newq0 < 0 is not supported: "//
+     $     "q must be positive. Use newq0 = 0 (default) or newq0 > 0.")
       IF(PRESENT(op_psihigh))THEN
          psihigh = op_psihigh
          IF(verbose) WRITE(*,*) "Reforming equilibrium with new psihigh"

equil/inverse.f

@@ -237,7 +237,6 @@ SUBROUTINE inverse_run
       ENDDO
       CALL spline_fit(sq,"extrap")
       q0=sq%fs(0,4)-sq%fs1(0,4)*sq%xs(0)
-      IF(newq0 == -1)newq0=-q0
 c-----------------------------------------------------------------------
 c     revise q profile.
 c-----------------------------------------------------------------------
@@ -246,7 +245,7 @@ SUBROUTINE inverse_run
          f0fac=f0**2*((newq0/q0)**2-one)
          q0=newq0
          DO ipsi=0,mpsi
-            ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2)*SIGN(one,newq0)
+            ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2)
             sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac
             sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac
             rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac
@@ -579,7 +578,11 @@ SUBROUTINE inverse_chease4_run(jaceq)
       ENDDO
       CALL spline_fit(sq,"extrap")
       q0=sq%fs(0,4)-sq%fs1(0,4)*sq%xs(0)
-      IF(newq0 == -1)newq0=-q0
+c-----------------------------------------------------------------------
+c     enforce positive q (chease4 takes q directly from the input file).
+c-----------------------------------------------------------------------
+      IF(MINVAL(sq%fs(:,4)) <= 0)CALL program_stop
+     $     ("chease4 input q must be positive (GPEC requires q > 0).")
 c-----------------------------------------------------------------------
 c     revise q profile.
 c-----------------------------------------------------------------------
@@ -588,7 +591,7 @@ SUBROUTINE inverse_chease4_run(jaceq)
          f0fac=f0**2*((newq0/q0)**2-one)
          q0=newq0
          DO ipsi=0,mpsi
-            ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2)*SIGN(one,newq0)
+            ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2)
             sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac
             sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac
             rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac

equil/read_eq.f

@@ -1498,13 +1498,12 @@ SUBROUTINE read_eq_dump
 c-----------------------------------------------------------------------
 c     revise q profile.
 c-----------------------------------------------------------------------
-      IF(newq0 == -1)newq0=-q0
       IF(newq0 /= 0)THEN
          f0=sq%fs(0,1)-sq%fs1(0,1)*sq%xs(0)
          f0fac=f0**2*((newq0/q0)**2-one)
          q0=newq0
          DO ipsi=0,mpsi
-            ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2)*SIGN(one,newq0)
+            ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2)
             sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac
             sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac
             rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac