Add eq.poloidalBeta, docs, version 0.1.5

Made some routines more robust to different inputs, particularly
`eq.q()`. This now handles case of no inputs, generating and returning
an array of normalised psi values.

Documetation on diagnostics

Added poloidal beta calculation

01-freeboundary.py

@@ -18,8 +18,8 @@
 #########################################
 # Plasma profiles
 
-profiles = freegs.jtor.ConstrainPaxisIp(1e4, # Plasma pressure on axis [Pascals]
-                                        1e6, # Plasma current [Amps]
+profiles = freegs.jtor.ConstrainPaxisIp(1e3, # Plasma pressure on axis [Pascals]
+                                        2e5, # Plasma current [Amps]
                                         2.0) # Vacuum f=R*Bt
 
 #########################################
@@ -55,6 +55,10 @@
 # Forces on the coils
 eq.printForces()
 
+print("\nSafety factor:\n\tpsi \t q")
+for psi in [0.01, 0.9, 0.95]:
+    print("\t{:.2f}\t{:.2f}".format(psi, eq.q(psi)))
+
 ##############################################
 # Save to G-EQDSK file
 
@@ -64,8 +68,15 @@
     geqdsk.write(eq, f)
 
 ##############################################
-# Final plot
+# Final plot of equilibrium
 
 axis = eq.plot(show=False)
 constrain.plot(axis=axis, show=True)
 
+# Safety factor
+
+import matplotlib.pyplot as plt
+plt.plot(*eq.q())
+plt.xlabel(r"Normalised $\psi$")
+plt.ylabel("Safety factor")
+plt.show()

05-fixed-boundary.py

@@ -10,8 +10,8 @@
 # Boundary conditions
 import freegs.boundary as boundary
 
-profiles = freegs.jtor.ConstrainPaxisIp(1e4, # Plasma pressure on axis [Pascals]
-                                        1e6, # Plasma current [Amps]
+profiles = freegs.jtor.ConstrainPaxisIp(1e3, # Plasma pressure on axis [Pascals]
+                                        1e5, # Plasma current [Amps]
                                         1.0) # fvac = R*Bt
 
 eq = freegs.Equilibrium(Rmin=0.1, Rmax=2.0,
@@ -25,9 +25,12 @@
 
 print("Done!")
 
+# Some diagnostics
+print("Poloidal beta: {}".format(eq.poloidalBeta()))
+print("Pressure on axis: {} Pa".format(eq.pressure(0.0)))
+
 # Plot equilibrium
 from freegs.plotting import plotEquilibrium
-
 plotEquilibrium(eq)
 
 

docs/creating_equilibria.rst

@@ -1,3 +1,5 @@
+.. _creating_equilibria
+
 Creating equilibria
 ===================
 
@@ -327,7 +329,9 @@ The total toroidal plasma current is calculated by integrating the toroidal curr
 
 The integrals in these two constraints are done numerically,
 and then rearranged to get :math:`L` and :math:`\beta_0`. 
-
+
+.. _constrain_betap_ip
+
 Constrain poloidal beta and current
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

docs/diagnostics.rst

@@ -1,20 +1,69 @@
 Diagnostics
 ===========
 
-Once an equilibrium has been generated (see `Creating equilibria`_)
+Once an equilibrium has been generated (see creating_equilibria_)
 there are routines for diagnosing and calculating derived quantities.
 Here the ``Equilibrium`` object is assumed to be called ``eq`` and the
 ``Tokamak`` object called ``tokamak``.
 
 Safety factor, q
 ----------------
 
-The safety factor at a given normalised poloidal flux (0 on the
-magnetic axis, 1 at the separatrix) can be calculated using the
-``q(psinorm)`` function::
+The safety factor :math:`q` at a given normalised poloidal flux
+:math:`\psi_{norm}` (0 on the magnetic axis, 1 at the separatrix) can
+be calculated using the ``q(psinorm)`` function::
 
   eq.q(0.9)  # safety factor at psi_norm = 0.9
 
+Note that calculating :math:`q` on either the magnetic axis or separatrix
+is problematic, so values calculated at :math:`\psi_{norm}=0` and :math:`\psi_{norm}=1`
+are likely to be inaccurate.
+
+This function can be used to print the safety factor on a set of flux
+surfaces::
+  
+  print("\nSafety factor:\n\tpsi \t q")
+  for psi in [0.01, 0.9, 0.95]:
+    print("\t{:.2f}\t{:.2f}".format(psi, eq.q(psi)))
+
+If no value is given for the normalised psi, then a uniform array of
+values between 0 and 1 is generated (not including the end points). In
+this case both the values of normalised psi and the values of q are returned::
+
+  psinorm, q = eq.q()
+
+which can be used to make a plot of the safety factor::
+
+  import matplotlib.pyplot as plt
+  
+  plt.plot(*eq.q())
+  plt.xlabel(r"Normalised $\psi$")
+  plt.ylabel(r"Safety factor $q$")
+  plt.show()
+
+Poloidal beta
+-------------
+
+The poloidal beta :math:`\beta_p` is given by::
+
+  betap = eq.poloidalBeta()
+
+This is calculated using the expression
+
+.. math::
+
+   \beta_p = \frac{8\pi}{\mu_0} \frac{1}{I_p^2}\iint p\left(\psi\right) dRdZ 
+
+i.e. the same calculation as is done in the poloidal beta constraint constrain_betap_ip_.
+
+Plasma pressure
+---------------
+
+The pressure at a specified normalised psi is::
+
+  p = eq.pressure(0.0)  # Pressure on axis
+
+
 Separatrix location
 -------------------
 
@@ -33,10 +82,19 @@ A set of points on the separatrix, measured in meters::
 Currents in the coils
 ---------------------
 
+The coil objects can be accessed and their currents queried. The
+current in a coil named "P1L" is given by::
+
+  eq.tokamak["P1L"].current
+
 The currents in all coils can be printed using::
 
-  tokamak.printCurrent()
+  tokamak.printCurrents()
+
+which is the same as::
 
+  for label, coil in eq.tokamak.coils:
+    print(label + " : " + str(coil))
 
 Forces on the coils
 -------------------

freegs/__init__.py

@@ -28,7 +28,7 @@
 
 """
 
-__version__ = "0.1.4"
+__version__ = "0.1.5"
 
 from .equilibrium import Equilibrium
 

freegs/critical.py

@@ -448,7 +448,11 @@ def find_safety(eq, npsi=1, psinorm=None, ntheta=128, psi=None, opoint=None, xpo
     if (opoint is None) or (xpoint is None):
         opoint, xpoint = find_critical(eq.R, eq.Z, psi)
 
-    psinormal = (psi - opoint[0][2])/(xpoint[0][2] - opoint[0][2])
+    if (xpoint is None) or (len(xpoint) == 0):
+        # No X-point
+        raise ValueError("No X-point so no separatrix")
+    else:
+        psinormal = (psi - opoint[0][2])/(xpoint[0][2] - opoint[0][2])
 
     psifunc = interpolate.RectBivariateSpline(eq.R[:,0], eq.Z[0,:], psinormal)
 
@@ -469,10 +473,14 @@ def find_safety(eq, npsi=1, psinorm=None, ntheta=128, psi=None, opoint=None, xpo
 
     if psinorm is None:
         npsi = 100
-        psirange = linspace(0.,1.0,npsi)
+        psirange = linspace(1./(npsi+1), 1.0, npsi, endpoint=False)
     else:
-        psirange = psinorm
-        npsi = len(psinorm)
+        try:
+            psirange = psinorm
+            npsi = len(psinorm)
+        except TypeError:
+            npsi = 1
+            psirange = [psinorm]
 
     psisurf = zeros([npsi,ntheta,2])
 

freegs/equilibrium.py

@@ -3,12 +3,13 @@
 state, including plasma and coil currents
 """
 
-from numpy import meshgrid, linspace, exp, zeros, nditer, array
+from numpy import pi, meshgrid, linspace, exp, zeros, nditer, array
+import numpy as np
 from scipy import interpolate
 from scipy.integrate import romb, quad # Romberg integration
 
 from .boundary import fixedBoundary, freeBoundary
-from .critical import find_separatrix, find_safety
+from . import critical
 
 # Operators which define the G-S equation
 from .gradshafranov import mu0, GSsparse
@@ -165,7 +166,48 @@ def plasmaCurrent(self):
         Plasma current [Amps]
         """
         return self._current
-
+
+    def poloidalBeta(self):
+        """ 
+        Return the poloidal beta 
+        betap = (8pi/mu0) * int(p)dRdZ / Ip^2
+        """
+
+        # Total poloidal flux (plasma + coils)
+        psi = self.psi()
+
+        # Analyse the equilibrium, finding O- and X-points
+        opt, xpt = critical.find_critical(self.R, self.Z, psi)
+        if not opt:
+            raise ValueError("No O-points found!")
+        psi_axis = opt[0][2]
+
+        if xpt:
+            psi_bndry = xpt[0][2]
+            mask = critical.core_mask(R, Z, psi, opt, xpt)
+        else:
+            # No X-points
+            if psi[0,0] > psi_axis:
+                psi_bndry = np.amax(psi)
+            else:
+                psi_bndry = np.amin(psi)
+            mask = None
+
+        dR = self.R[1,0] - self.R[0,0]
+        dZ = self.Z[0,1] - self.Z[0,0]
+
+        # Normalised psi
+        psi_norm = (psi - psi_axis)  / (psi_bndry - psi_axis)
+
+        # Plasma pressure
+        pressure = self.pressure(psi_norm)
+        if mask is not None:
+            # If there is a masking function (X-points, limiters)
+            pressure *= mask
+
+        # Integrate pressure in 2D
+        return ((8.*pi)/mu0)*romb(romb(pressure))*dR*dZ / (self.plasmaCurrent()**2)
+
     def plasmaBr(self, R,Z):
         """
         Radial magnetic field due to plasma
@@ -214,11 +256,33 @@ def fvac(self):
         """
         return self._profiles.fvac()
 
-    def q(self, psinorm):
+    def q(self, psinorm = None, npsi=100):
         """
         Returns safety factor q at specified values of normalised psi
-        """
-        return find_safety(self,psinorm=psinorm)
+        
+        psinorm is a scalar, list or array of floats betweem 0 and 1.
+        
+        >>> safety_factor = eq.q([0.2, 0.5, 0.9])
+        
+        If psinorm is None, then q on a uniform psi grid will be returned,
+        along with the psi values
+
+        >>> psinorm, q = eq.q()
+        
+        Note: psinorm = 0 is the magnetic axis, and psinorm = 1 is the separatrix.
+              Calculating q on either of these flux surfaces is problematic,
+              and the results will probably not be accurate.
+        """
+        if psinorm is None:
+            # An array which doesn't include psinorm = 0 or 1
+            psinorm = linspace(1./(npsi+1), 1.0, npsi, endpoint=False)
+            return psinorm, critical.find_safety(self, psinorm=psinorm)
+
+        result = critical.find_safety(self, psinorm=psinorm)
+        # Convert to a scalar if only one result
+        if len(result) == 1:
+            return np.asscalar(result)
+        return result
 
     def pprime(self, psinorm):
         """
@@ -243,7 +307,7 @@ def separatrix(self, ntheta=20):
         Returns an array of ntheta (R, Z) coordinates of the separatrix,
         equally spaced in geometric poloidal angle.
         """
-        return array(find_separatrix(self, ntheta=ntheta, psi=self.psi()))[:, 0:2]
+        return array(critical.find_separatrix(self, ntheta=ntheta, psi=self.psi()))[:, 0:2]
 
     def solve(self, profiles, Jtor=None, psi=None, psi_bndry=None):
         """
@@ -341,7 +405,7 @@ def print_forces(forces, prefix=""):
                     print(prefix + label + " (circuit)")
                     print_forces(force, prefix=prefix + "  ")
                 else:
-                    print(prefix + label+ " : R = {0:.2f} MN , Z = {1:.2f} MN".format(force[0]*1e-6, force[1]*1e-6))
+                    print(prefix + label+ " : R = {0:.2f} kN , Z = {1:.2f} kN".format(force[0]*1e-3, force[1]*1e-3))
 
         print_forces(self.getForces())
 

freegs/jtor.py

@@ -40,7 +40,7 @@
 from .gradshafranov import mu0
 
 from numpy import clip, zeros, reshape, sqrt, pi
-
+import numpy as np
 
 class Profile(object):
     """
@@ -90,8 +90,8 @@ def fpol(self, psinorm, out=None):
         
         """
 
-        if not hasattr(psinorm, 'shape'):
-            # Assume  a single value
+        if not hasattr(psinorm, '__len__'):
+            # Assume a single value
 
             val, _ = quad(self.ffprime, psinorm, 1.0)
             # Convert from integral in normalised psi to integral in psi
@@ -101,7 +101,8 @@ def fpol(self, psinorm, out=None):
             # Apply boundary condition at psinorm=1 val = fvac**2
             return sqrt(2.*val + self.fvac()**2)
 
-        # Assume it's a NumPy array
+        # Assume it's a NumPy array, or can be converted to one
+        psinorm = np.array(psinorm)
 
         if out is None:
             out = zeros(psinorm.shape)

freegs/machine.py

@@ -203,7 +203,7 @@ def getForces(self, equilibrium):
                           -total_current * Br * Ltor]) # Jphi x Br = - Fz
 
     def __repr__(self):
-        return ("Coil(R={0}, Z={1}, current={2}, turns={3}, control={4})"
+        return ("Coil(R={0}, Z={1}, current={2:.1f}, turns={3}, control={4})"
                 .format(self.R, self.Z, self.current, self.turns, self.control))
 
     def __eq__(self, other):