Merge pull request #121 from freegs-plasma/R0-geqdsk

Add R0 machine property for GEQDSK output

Author: bendudson

freegs/geqdsk.py

@@ -4,7 +4,7 @@
 Writing is relatively straightforward, but reading requires inferring
 the currents in the PF coils
 
-Copyright 2016 Ben Dudson, University of York. Email: [email protected]
+Copyright 2016-2025 FreeGS contributors
 
 This file is part of FreeGS.
 
@@ -85,7 +85,7 @@ def write(eq, fh, label=None, oxpoints=None, fileformat=geqdsk.write):
     zmax = eq.Zmax
 
     fvac = eq.fvac()  # Vacuum f = R*Bt
-    R0 = 1.0  # Reference location
+    R0 = eq.tokamak.R0  # Reference location
 
     data = {
         "nx": nx,
@@ -337,6 +337,9 @@ def q_func(psinorm):
     # 1 = plasma boundary
     psi_norm = clip((psi - psi_axis) / (psi_bndry - psi_axis), 0.0, 1.1)
 
+    # rcentr reference location
+    machine.R0 = data["rcentr"]
+
     # Create an Equilibrium object
     eq = Equilibrium(
         tokamak=machine,

freegs/machine.py

@@ -4,7 +4,7 @@
 License
 -------
 
-Copyright 2016-2019 Ben Dudson, University of York. Email: [email protected]
+Copyright 2016-2025 FreeGS contributors.
 
 This file is part of FreeGS.
 
@@ -637,15 +637,18 @@ class Machine:
 
     """
 
-    def __init__(self, coils, wall=None, sensors=None, nlimit=500):
+    def __init__(self, coils, wall=None, sensors=None, nlimit=500, R0=1.0):
         """
         coils - A list of coils [(label, Coil|Circuit|Solenoid)]
         sensors - A list of sensors
+        R0 - The major radius to be written to GEQDSK [meters].
+             By convention EFIT has a value of R0 for each tokamak.
         """
 
         self.coils = coils
         self.wall = wall
         self.sensors = sensors
+        self.R0 = R0
 
         self.limit_points_R = None
         self.limit_points_Z = None
@@ -1018,7 +1021,8 @@ def DIIID():
         -1.00000e-03,0.00000e+00])
     wall = Wall(R,Z)
 
-    return Machine(coils,wall=wall)
+    # Note: Using same R0 convention as EFIT on DIII-D.
+    return Machine(coils, wall=wall, R0=1.6995)
 
 
 def MAST():