From 2b7bf0404134fcc40334a7109f3b6eec1f5dc589 Mon Sep 17 00:00:00 2001
From: Christopher Albert <albert@tugraz.at>
Date: Fri, 23 Oct 2026 23:46:54 +0200
Subject: [PATCH] Equilibrium - BUGFIX - Preserve projected radii in direct
 flux coordinates

(cherry picked from commit ff6edda2a6415000c27e6d88338b7fbc737284a0)
---
 equil/direct.f                                |  3 ++
 regression/test/direct_projected_radius.f90   | 35 +++++++++++++++
 .../test/test_direct_projected_radius.py      | 45 +++++++++++++++++++
 3 files changed, 83 insertions(+)
 create mode 100644 regression/test/direct_projected_radius.f90
 create mode 100644 regression/test/test_direct_projected_radius.py

diff --git a/equil/direct.f b/equil/direct.f
index c81c4aef..299d2dab 100644
--- a/equil/direct.f
+++ b/equil/direct.f
@@ -159,6 +159,8 @@ SUBROUTINE direct_run
          CALL direct_fl_int(ipsi,y_out,bf,y4tot,ys)
          DO itheta=0,mtheta
             theta=rzphi%ys(itheta)
+            CALL direct_refine(ys(itheta,2),ys(itheta,0),
+     $           psio*(1-sq%xs(ipsi)))
             CALL direct_fl_der(4,ys(itheta,0),ys(itheta,1:4),dy)
             rzphi%fs(ipsi,itheta,1)=ys(itheta,2)**2
             rzphi%fs(ipsi,itheta,2)=ys(itheta,0)/twopi-theta
@@ -513,6 +515,7 @@ SUBROUTINE direct_fl_int(ipsi,y_out,bf,y4tot,ys)
          z=zo+rfac*SIN(eta)
          CALL direct_get_bfield(r,z,bf,2)
          y_out(istep,:)=(/eta,y/)
+         y_out(istep,2)=rfac
          err=(bf%psi-psi0)/bf%psi
 c-----------------------------------------------------------------------
 c     compute and print output for each step.
diff --git a/regression/test/direct_projected_radius.f90 b/regression/test/direct_projected_radius.f90
new file mode 100644
index 00000000..6ef143a2
--- /dev/null
+++ b/regression/test/direct_projected_radius.f90
@@ -0,0 +1,35 @@
+program direct_projected_radius
+  use read_eq_mod
+  implicit none
+  integer :: j, idx
+  real(r8) :: radius, angle, r, z, expected, residual
+
+  call get_command_argument(1,eq_filename)
+  mpsi=32
+  mtheta=64
+  psilow=.05_r8
+  psihigh=.95_r8
+  etol=1e-5_r8
+  grid_type='ldp'
+  jac_type='hamada'
+  power_bp=0
+  power_b=0
+  power_r=0
+  verbose=.false.
+  open(out_unit,file='gpec.out',status='new')
+  call read_eq_efit()
+
+  idx=mpsi/2
+  expected=psio*(1-sq%xs(idx))
+  residual=0
+  do j=0,mtheta-1
+    call bicube_eval(rzphi,sq%xs(idx),real(j,r8)/mtheta,0)
+    radius=sqrt(rzphi%f(1))
+    angle=twopi*(real(j,r8)/mtheta+rzphi%f(2))
+    r=ro+radius*cos(angle)
+    z=zo+radius*sin(angle)
+    call bicube_eval(psi_in,r,z,0)
+    residual=max(residual,abs(psi_in%f(1)-expected)/psio)
+  end do
+  write(*,'(es25.16)') residual
+end program
diff --git a/regression/test/test_direct_projected_radius.py b/regression/test/test_direct_projected_radius.py
new file mode 100644
index 00000000..7549fec8
--- /dev/null
+++ b/regression/test/test_direct_projected_radius.py
@@ -0,0 +1,45 @@
+"""Check that the delivered direct chart stays on its input flux surface."""
+import os
+from pathlib import Path
+import shutil
+import subprocess
+import tempfile
+import unittest
+
+
+class DirectProjectedRadius(unittest.TestCase):
+    @unittest.skipUnless(shutil.which('gfortran'), 'gfortran is required')
+    def test_chart_radius_follows_input_flux_surface(self):
+        root = Path(__file__).resolve().parents[2]
+        scratch = Path(os.environ.get('GPEC_TEST_SCRATCH', root / 'build-test'))
+        scratch.mkdir(exist_ok=True)
+        names = ('io', 'local', 'utils', 'spline', 'bicube', 'cspline',
+                 'global', 'grid', 'inverse', 'direct', 'read_eq')
+        sources = [root / 'equil' / (name + ('.f90' if name == 'grid' else '.f'))
+                   for name in names]
+        direct_override = os.environ.get('GPEC_DIRECT_SOURCE')
+        if direct_override:
+            sources[names.index('direct')] = Path(direct_override)
+        sources += sorted((root / 'lsode').glob('*.f'))
+        fixture = root / 'docs/examples/a5_tearing_example/gfile_a5_cir'
+        driver = root / 'regression/test/direct_projected_radius.f90'
+        with tempfile.TemporaryDirectory(dir=scratch) as directory:
+            work = Path(directory)
+            for source in sources:
+                subprocess.run(['gfortran', '-O2', '-fallow-argument-mismatch',
+                                '-frecursive', '-ffixed-line-length-none', '-c',
+                                str(source)], cwd=work, check=True,
+                               capture_output=True)
+            subprocess.run(['gfortran', '-O2', str(driver),
+                            *map(str, sorted(work.glob('*.o'))),
+                            '-llapack', '-lblas', '-o', str(work / 'probe')],
+                           cwd=work, check=True, capture_output=True)
+            result = subprocess.run([str(work / 'probe'), str(fixture)],
+                                    cwd=work, check=True, capture_output=True,
+                                    text=True)
+            residual = float(result.stdout.splitlines()[-1])
+            self.assertLess(residual, 1e-5)
+
+
+if __name__ == '__main__':
+    unittest.main()