From 79223666fde28bed663b4729d0e9aa48ae60f6c8 Mon Sep 17 00:00:00 2001
From: Aaron Bader <aaron.bader@typeoneenergy.com>
Date: Mon, 26 Feb 2024 14:34:10 -0600
Subject: [PATCH 1/8] add magnetic potential to mgrid

---
 src/simsopt/field/magneticfield.py |  9 ++++++++-
 src/simsopt/field/mgrid.py         | 21 ++++++++++++++++++++-
 2 files changed, 28 insertions(+), 2 deletions(-)

diff --git a/src/simsopt/field/magneticfield.py b/src/simsopt/field/magneticfield.py
index cfe03ba68..62aa7d89c 100644
--- a/src/simsopt/field/magneticfield.py
+++ b/src/simsopt/field/magneticfield.py
@@ -127,17 +127,24 @@ def to_mgrid(self, filename, nr=10, nphi=4, nz=12, rmin=1.0, rmax=2.0, zmin=-0.5
         # get field on the grid
         self.set_points_cyl(RPhiZ)
         B = self.B_cyl()
+        A = self.A_cyl()
 
         # shape the components
         br, bp, bz = B.T
         br_3 = br.reshape((nphi, nz, nr))
         bp_3 = bp.reshape((nphi, nz, nr))
         bz_3 = bz.reshape((nphi, nz, nr))
+        
+        # shape the Potential components
+        ar, ap, az = A.T
+        ar_3 = ar.reshape((nphi, nz, nr))
+        ap_3 = ap.reshape((nphi, nz, nr))
+        az_3 = az.reshape((nphi, nz, nr))
 
         mgrid = MGrid(nfp=nfp,
                       nr=nr, nz=nz, nphi=nphi,
                       rmin=rmin, rmax=rmax, zmin=zmin, zmax=zmax)
-        mgrid.add_field_cylindrical(br_3, bp_3, bz_3, name='simsopt_coils')  
+        mgrid.add_field_cylindrical(br_3, bp_3, bz_3, ar_3, ap_3, az_3, name='simsopt_coils')  
 
         mgrid.write(filename)
 
diff --git a/src/simsopt/field/mgrid.py b/src/simsopt/field/mgrid.py
index bb946c18b..6ca39e333 100644
--- a/src/simsopt/field/mgrid.py
+++ b/src/simsopt/field/mgrid.py
@@ -69,7 +69,11 @@ def __init__(self,  # fname='temp', #binary=False,
         self.bz_arr = [] 
         self.bp_arr = [] 
 
-    def add_field_cylindrical(self, br, bp, bz, name=None):
+        self.ar_arr = []
+        self.az_arr = []
+        self.ap_arr = []
+
+    def add_field_cylindrical(self, br, bp, bz, ar, ap, az, name=None):
         '''
         This function saves the vector field B.
         B is defined by cylindrical components.
@@ -88,6 +92,9 @@ def add_field_cylindrical(self, br, bp, bz, name=None):
             br: the radial component of B field
             bp: the azimuthal component of B field
             bz: the axial component of B field
+            ar: the radial component of the magnetic potential
+            ap: the azimuthal component of the magnetic potential
+            az: the axial component of the magnetic potential
             name: Name of the coil group
         '''
 
@@ -95,6 +102,11 @@ def add_field_cylindrical(self, br, bp, bz, name=None):
         self.br_arr.append(br)
         self.bz_arr.append(bz)
         self.bp_arr.append(bp)
+        
+        # add potential
+        self.ar_arr.append(ar)
+        self.az_arr.append(az)
+        self.ap_arr.append(ap)
 
         # add coil label
         if (name is None):
@@ -104,6 +116,7 @@ def add_field_cylindrical(self, br, bp, bz, name=None):
         self.coil_names.append(label)
         self.n_ext_cur = self.n_ext_cur + 1
 
+
         # TO-DO: this function could check for size consistency, between different fields, and for the (nr,nphi,nz) settings of a given instance
 
     def write(self, filename):
@@ -172,10 +185,16 @@ def write(self, filename):
                 var_br_001 = ds.createVariable('br'+tag, 'f8', ('phi', 'zee', 'rad'))
                 var_bp_001 = ds.createVariable('bp'+tag, 'f8', ('phi', 'zee', 'rad'))
                 var_bz_001 = ds.createVariable('bz'+tag, 'f8', ('phi', 'zee', 'rad'))
+                var_ar_001 = ds.createVariable('ar'+tag, 'f8', ('phi', 'zee', 'rad'))
+                var_ap_001 = ds.createVariable('ap'+tag, 'f8', ('phi', 'zee', 'rad'))
+                var_az_001 = ds.createVariable('az'+tag, 'f8', ('phi', 'zee', 'rad'))
 
                 var_br_001[:, :, :] = self.br_arr[j]
                 var_bz_001[:, :, :] = self.bz_arr[j]
                 var_bp_001[:, :, :] = self.bp_arr[j]
+                var_ar_001[:, :, :] = self.ar_arr[j]
+                var_az_001[:, :, :] = self.az_arr[j]
+                var_ap_001[:, :, :] = self.ap_arr[j]
 
     @classmethod 
     def from_file(cls, filename):

From cb9e3ef2c95cdd8964d082e8dc805872375bdba6 Mon Sep 17 00:00:00 2001
From: Aaron Bader <aaron.bader@typeoneenergy.com>
Date: Mon, 26 Feb 2024 16:24:10 -0600
Subject: [PATCH 2/8] update tests

---
 tests/field/test_magneticfields.py | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/tests/field/test_magneticfields.py b/tests/field/test_magneticfields.py
index 3985708c7..3b2eb469c 100644
--- a/tests/field/test_magneticfields.py
+++ b/tests/field/test_magneticfields.py
@@ -1062,12 +1062,18 @@ def test_to_mgrid(self):
                 Br = f.variables["br_001"][()]
                 Bphi = f.variables["bp_001"][()]
                 Bz = f.variables["bz_001"][()]
+                Ar = f.variables["ar_001"][()]
+                Aphi = f.variables["ap_001"][()]
+                Az = f.variables["az_001"][()]
                 assert nr == f.dimensions["rad"]
                 assert nphi == f.dimensions["phi"]
                 assert nz == f.dimensions["zee"]
                 assert Br.shape == (nphi, nz, nr)
                 assert Bphi.shape == (nphi, nz, nr)
                 assert Bz.shape == (nphi, nz, nr)
+                assert Ar.shape == (nphi, nz, nr)
+                assert Aphi.shape == (nphi, nz, nr)
+                assert Az.shape == (nphi, nz, nr)
                 r = np.linspace(rmin, rmax, nr)
                 phi = np.linspace(0, 2 * np.pi / nfp, nphi, endpoint=False)
                 z = np.linspace(zmin, zmax, nz)
@@ -1078,6 +1084,9 @@ def test_to_mgrid(self):
                             np.testing.assert_allclose(Br[jphi, jz, jr], bs.B_cyl()[0, 0])
                             np.testing.assert_allclose(Bphi[jphi, jz, jr], bs.B_cyl()[0, 1])
                             np.testing.assert_allclose(Bz[jphi, jz, jr], bs.B_cyl()[0, 2])
+                            np.testing.assert_allclose(Ar[jphi, jz, jr], bs.A_cyl()[0, 0])
+                            np.testing.assert_allclose(Aphi[jphi, jz, jr], bs.A_cyl()[0, 1])
+                            np.testing.assert_allclose(Az[jphi, jz, jr], bs.A_cyl()[0, 2])
 
     def test_poloidal_field(self):
         B0 = 1.1

From de04ed9a15fa924b6cbbac01a6c1f51d69fa6eee Mon Sep 17 00:00:00 2001
From: Aaron Bader <aaron.bader@typeoneenergy.com>
Date: Thu, 4 Apr 2024 10:46:07 -0500
Subject: [PATCH 3/8] add potential as an option and update tests

---
 src/simsopt/field/magneticfield.py | 23 +++++++++++++++--------
 src/simsopt/field/mgrid.py         | 30 +++++++++++++++++++-----------
 tests/field/test_magneticfields.py |  4 ++--
 tests/field/test_mgrid.py          |  2 +-
 4 files changed, 37 insertions(+), 22 deletions(-)

diff --git a/src/simsopt/field/magneticfield.py b/src/simsopt/field/magneticfield.py
index 62aa7d89c..9897244a2 100644
--- a/src/simsopt/field/magneticfield.py
+++ b/src/simsopt/field/magneticfield.py
@@ -92,7 +92,8 @@ def to_vtk(self, filename, nr=10, nphi=10, nz=10, rmin=1.0, rmax=2.0, zmin=-0.5,
         contig = np.ascontiguousarray
         gridToVTK(filename, X, Y, Z, pointData={"B": (contig(vals[..., 0]), contig(vals[..., 1]), contig(vals[..., 2]))})
 
-    def to_mgrid(self, filename, nr=10, nphi=4, nz=12, rmin=1.0, rmax=2.0, zmin=-0.5, zmax=0.5, nfp=1):
+    def to_mgrid(self, filename, nr=10, nphi=4, nz=12, rmin=1.0, rmax=2.0, zmin=-0.5, zmax=0.5, nfp=1, 
+                 include_potential=False):
         """Export the field to the mgrid format for free boundary calculations.
 
         The field data is represented as a single "current group". For
@@ -111,6 +112,7 @@ def to_mgrid(self, filename, nr=10, nphi=4, nz=12, rmin=1.0, rmax=2.0, zmin=-0.5
             zmin: Minimum value of z for the grid.
             zmax: Maximum value of z for the grid.
             nfp: Number of field periods.
+            include_potential: Boolean to include the magnetic potential, defaults to false.
         """
 
         rs = np.linspace(rmin, rmax, nr, endpoint=True)
@@ -127,7 +129,6 @@ def to_mgrid(self, filename, nr=10, nphi=4, nz=12, rmin=1.0, rmax=2.0, zmin=-0.5
         # get field on the grid
         self.set_points_cyl(RPhiZ)
         B = self.B_cyl()
-        A = self.A_cyl()
 
         # shape the components
         br, bp, bz = B.T
@@ -135,16 +136,22 @@ def to_mgrid(self, filename, nr=10, nphi=4, nz=12, rmin=1.0, rmax=2.0, zmin=-0.5
         bp_3 = bp.reshape((nphi, nz, nr))
         bz_3 = bz.reshape((nphi, nz, nr))
         
-        # shape the Potential components
-        ar, ap, az = A.T
-        ar_3 = ar.reshape((nphi, nz, nr))
-        ap_3 = ap.reshape((nphi, nz, nr))
-        az_3 = az.reshape((nphi, nz, nr))
+        if include_potential:
+            A = self.A_cyl()
+            # shape the Potential components
+            ar, ap, az = A.T
+            ar_3 = ar.reshape((nphi, nz, nr))
+            ap_3 = ap.reshape((nphi, nz, nr))
+            az_3 = az.reshape((nphi, nz, nr))
 
         mgrid = MGrid(nfp=nfp,
                       nr=nr, nz=nz, nphi=nphi,
                       rmin=rmin, rmax=rmax, zmin=zmin, zmax=zmax)
-        mgrid.add_field_cylindrical(br_3, bp_3, bz_3, ar_3, ap_3, az_3, name='simsopt_coils')  
+        if include_potential:
+            mgrid.add_field_cylindrical(br_3, bp_3, bz_3, ar=ar_3, ap=ap_3, az=az_3, name='simsopt_coils')  
+        else:
+            mgrid.add_field_cylindrical(br_3, bp_3, bz_3, name='simsopt_coils')  
+
 
         mgrid.write(filename)
 
diff --git a/src/simsopt/field/mgrid.py b/src/simsopt/field/mgrid.py
index 6ca39e333..982e60c16 100644
--- a/src/simsopt/field/mgrid.py
+++ b/src/simsopt/field/mgrid.py
@@ -73,7 +73,7 @@ def __init__(self,  # fname='temp', #binary=False,
         self.az_arr = []
         self.ap_arr = []
 
-    def add_field_cylindrical(self, br, bp, bz, ar, ap, az, name=None):
+    def add_field_cylindrical(self, br, bp, bz, ar=None, ap=None, az=None, name=None):
         '''
         This function saves the vector field B.
         B is defined by cylindrical components.
@@ -104,9 +104,12 @@ def add_field_cylindrical(self, br, bp, bz, ar, ap, az, name=None):
         self.bp_arr.append(bp)
         
         # add potential
-        self.ar_arr.append(ar)
-        self.az_arr.append(az)
-        self.ap_arr.append(ap)
+        if ar is not None:
+            self.ar_arr.append(ar)
+        if ap is not None:
+            self.az_arr.append(az)
+        if az is not None:
+            self.ap_arr.append(ap)
 
         # add coil label
         if (name is None):
@@ -119,7 +122,7 @@ def add_field_cylindrical(self, br, bp, bz, ar, ap, az, name=None):
 
         # TO-DO: this function could check for size consistency, between different fields, and for the (nr,nphi,nz) settings of a given instance
 
-    def write(self, filename):
+    def write(self, filename): 
         '''
         Export class data as a netCDF binary.
 
@@ -185,16 +188,21 @@ def write(self, filename):
                 var_br_001 = ds.createVariable('br'+tag, 'f8', ('phi', 'zee', 'rad'))
                 var_bp_001 = ds.createVariable('bp'+tag, 'f8', ('phi', 'zee', 'rad'))
                 var_bz_001 = ds.createVariable('bz'+tag, 'f8', ('phi', 'zee', 'rad'))
-                var_ar_001 = ds.createVariable('ar'+tag, 'f8', ('phi', 'zee', 'rad'))
-                var_ap_001 = ds.createVariable('ap'+tag, 'f8', ('phi', 'zee', 'rad'))
-                var_az_001 = ds.createVariable('az'+tag, 'f8', ('phi', 'zee', 'rad'))
 
                 var_br_001[:, :, :] = self.br_arr[j]
                 var_bz_001[:, :, :] = self.bz_arr[j]
                 var_bp_001[:, :, :] = self.bp_arr[j]
-                var_ar_001[:, :, :] = self.ar_arr[j]
-                var_az_001[:, :, :] = self.az_arr[j]
-                var_ap_001[:, :, :] = self.ap_arr[j]
+    
+                #If the potential value is not an empty cell, then include it
+                if len(self.ar_arr) > 0 :
+                    var_ar_001 = ds.createVariable('ar'+tag, 'f8', ('phi', 'zee', 'rad'))
+                    var_ar_001[:, :, :] = self.ar_arr[j]
+                if len(self.ap_arr) > 0 :
+                    var_ap_001 = ds.createVariable('ap'+tag, 'f8', ('phi', 'zee', 'rad'))
+                    var_ap_001[:, :, :] = self.ap_arr[j]
+                if len(self.az_arr) > 0 :
+                    var_az_001 = ds.createVariable('az'+tag, 'f8', ('phi', 'zee', 'rad'))
+                    var_az_001[:, :, :] = self.az_arr[j]
 
     @classmethod 
     def from_file(cls, filename):
diff --git a/tests/field/test_magneticfields.py b/tests/field/test_magneticfields.py
index 3b2eb469c..965ab01bb 100644
--- a/tests/field/test_magneticfields.py
+++ b/tests/field/test_magneticfields.py
@@ -1041,14 +1041,14 @@ def test_to_vtk(self):
         bs = BiotSavart(coils)
         bs.to_vtk('/tmp/bfield')
 
-    def test_to_mgrid(self):
+    def test_to_mgrid(self, include_potential=True):
         curves, currents, ma = get_ncsx_data()
         nfp = 3
         coils = coils_via_symmetries(curves, currents, nfp, True)
         bs = BiotSavart(coils)
         with tempfile.TemporaryDirectory() as tmpdir:
             filename = Path(tmpdir) / "mgrid.bfield.nc"
-            bs.to_mgrid(filename, nfp=nfp)
+            bs.to_mgrid(filename, nfp=nfp, include_potential=True)
 
             # Compare the B data in the file to a separate evaluation here
             with netcdf_file(filename, mmap=False) as f:
diff --git a/tests/field/test_mgrid.py b/tests/field/test_mgrid.py
index 8a913f07d..6b9bf63f8 100644
--- a/tests/field/test_mgrid.py
+++ b/tests/field/test_mgrid.py
@@ -49,7 +49,7 @@ def test_add_field_cylinder(self):
         assert mgrid.coil_names[0] == '__________test_coil___________'
         assert np.allclose(mgrid.br_arr[0], br)
 
-    def test_write(self):
+    def test_write(self): 
         mgrid = MGrid.from_file(test_file)
         with tempfile.TemporaryDirectory() as tmpdir:
             filename = Path(tmpdir) / 'mgrid.test.nc'

From 157839aec7105f99a12b5413beb7714831291d2b Mon Sep 17 00:00:00 2001
From: Aaron Bader <aaron.bader@typeoneenergy.com>
Date: Thu, 4 Apr 2024 11:01:30 -0500
Subject: [PATCH 4/8] adjust the mgrid reading file to handle mgrids with
 potential

---
 src/simsopt/field/mgrid.py | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/src/simsopt/field/mgrid.py b/src/simsopt/field/mgrid.py
index 982e60c16..bed84eb34 100644
--- a/src/simsopt/field/mgrid.py
+++ b/src/simsopt/field/mgrid.py
@@ -239,6 +239,10 @@ def from_file(cls, filename):
             bp_arr = []
             bz_arr = []
 
+            ar_arr = []
+            ap_arr = []
+            az_arr = []
+
             nextcur = mgrid.n_ext_cur
             for j in range(nextcur):
                 idx = '{:03d}'.format(j+1)
@@ -255,6 +259,27 @@ def from_file(cls, filename):
                     bp_arr.append(bp)
                     bz_arr.append(bz)
 
+                #check for potential in mgrid file
+                if 'ar_'+idx in f.variables:
+                    ar = f.variables['ar_'+idx][:]
+                    if mgrid.mode == 'S':
+                        ar_arr.append(ar * mgrid.raw_coil_current[j])
+                    else:
+                        ar_arr.append(ar)
+                        
+                if 'ap_'+idx in f.variables:
+                    ap = f.variables['ap_'+idx][:]
+                    if mgrid.mode == 'S':
+                        ap_arr.append(ap * mgrid.raw_coil_current[j])
+                    else:
+                        ap_arr.append(ap)
+                if 'az_'+idx in f.variables:
+                    az = f.variables['az_'+idx][:]
+                    if mgrid.mode == 'S':
+                        az_arr.append(az * mgrid.raw_coil_current[j])
+                    else:
+                        az_arr.append(az)
+
             mgrid.br_arr = np.array(br_arr)
             mgrid.bp_arr = np.array(bp_arr)
             mgrid.bz_arr = np.array(bz_arr)

From 054319cc8c006a80c640d59099319de690b43eff Mon Sep 17 00:00:00 2001
From: Aaron Bader <aaron.bader@typeoneenergy.com>
Date: Fri, 5 Apr 2024 15:06:18 -0500
Subject: [PATCH 5/8] allow loading mgrids when coil_groups has two dimensions
 or 1

---
 src/simsopt/field/mgrid.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/simsopt/field/mgrid.py b/src/simsopt/field/mgrid.py
index bed84eb34..a6ec67dbd 100644
--- a/src/simsopt/field/mgrid.py
+++ b/src/simsopt/field/mgrid.py
@@ -230,7 +230,10 @@ def from_file(cls, filename):
             mgrid.filename = filename
             mgrid.n_ext_cur = int(f.variables['nextcur'].getValue())
             coil_data = f.variables['coil_group'][:]
-            mgrid.coil_names = [_unpack(coil_data[j]) for j in range(mgrid.n_ext_cur)] 
+            if len(f.variables['coil_group'].dimensions) == 2:
+                mgrid.coil_names = [_unpack(coil_data[j]) for j in range(mgrid.n_ext_cur)] 
+            else:
+                mgrid.coil_names = [_unpack(coil_data)]
 
             mgrid.mode = f.variables['mgrid_mode'][:][0].decode()
             mgrid.raw_coil_current = np.array(f.variables['raw_coil_cur'][:])

From 1ef6e4ffe85ddaa2e407d62e1d01d3752c454273 Mon Sep 17 00:00:00 2001
From: Aaron Bader <aaron.bader@typeoneenergy.com>
Date: Fri, 5 Apr 2024 15:19:47 -0500
Subject: [PATCH 6/8] fix loading of potential values in mgrids

---
 src/simsopt/field/mgrid.py | 26 ++++++++++++++++++++++++++
 1 file changed, 26 insertions(+)

diff --git a/src/simsopt/field/mgrid.py b/src/simsopt/field/mgrid.py
index a6ec67dbd..01c6ce34a 100644
--- a/src/simsopt/field/mgrid.py
+++ b/src/simsopt/field/mgrid.py
@@ -246,6 +246,10 @@ def from_file(cls, filename):
             ap_arr = []
             az_arr = []
 
+            ar = []
+            ap = []
+            az = []
+
             nextcur = mgrid.n_ext_cur
             for j in range(nextcur):
                 idx = '{:03d}'.format(j+1)
@@ -276,6 +280,7 @@ def from_file(cls, filename):
                         ap_arr.append(ap * mgrid.raw_coil_current[j])
                     else:
                         ap_arr.append(ap)
+                
                 if 'az_'+idx in f.variables:
                     az = f.variables['az_'+idx][:]
                     if mgrid.mode == 'S':
@@ -287,20 +292,41 @@ def from_file(cls, filename):
             mgrid.bp_arr = np.array(bp_arr)
             mgrid.bz_arr = np.array(bz_arr)
 
+            mgrid.ar_arr = np.array(ar_arr)
+            mgrid.ap_arr = np.array(ap_arr)
+            mgrid.az_arr = np.array(az_arr)
+            
             # sum over coil groups
             if nextcur > 1:
                 br = np.sum(br_arr, axis=0)
                 bp = np.sum(bp_arr, axis=0)
                 bz = np.sum(bz_arr, axis=0)
+                
+                if len(mgrid.ar_arr) > 0:
+                    ar = np.sum(ar_arr, axis=0)
+                if len(mgrid.ap_arr) > 0:
+                    ap = np.sum(ap_arr, axis=0)
+                if len(mgrid.az_arr) > 0:
+                    az = np.sum(az_arr, axis=0)
             else:
                 br = br_arr[0]
                 bp = bp_arr[0]
                 bz = bz_arr[0]
+                if len(mgrid.ar_arr) > 0:
+                    ar = ar_arr[0]
+                if len(mgrid.ap_arr) > 0:
+                    ap = ap_arr[0]
+                if len(mgrid.az_arr) > 0:
+                    az = az_arr[0]
 
             mgrid.br = br
             mgrid.bp = bp
             mgrid.bz = bz
 
+            mgrid.ar = ar
+            mgrid.ap = ap
+            mgrid.az = az
+
             mgrid.bvec = np.transpose([br, bp, bz])
 
         return mgrid

From 5b53d27e0b2ff81e5c4c00a97999c02140e56a4b Mon Sep 17 00:00:00 2001
From: Aaron Bader <aaron.bader@typeoneenergy.com>
Date: Fri, 5 Apr 2024 15:20:15 -0500
Subject: [PATCH 7/8] add test of mgrid with potential

---
 tests/field/test_mgrid.py                  |   9 ++++++++-
 tests/test_files/mgrid_ncsx_lowres_test.nc | Bin 0 -> 24008 bytes
 2 files changed, 8 insertions(+), 1 deletion(-)
 create mode 100644 tests/test_files/mgrid_ncsx_lowres_test.nc

diff --git a/tests/field/test_mgrid.py b/tests/field/test_mgrid.py
index 6b9bf63f8..0799f051f 100644
--- a/tests/field/test_mgrid.py
+++ b/tests/field/test_mgrid.py
@@ -18,7 +18,7 @@
 
 TEST_DIR = (Path(__file__).parent / ".." / "test_files").resolve()
 test_file = TEST_DIR / 'mgrid.pnas-qa-test-lowres-standard.nc'
-
+test_file2 = TEST_DIR / 'mgrid_ncsx_lowres_test.nc'
 
 class Testing(unittest.TestCase):
 
@@ -36,6 +36,13 @@ def test_from_file(self):
         assert mgrid.br_arr.shape == (1, 6, 11, 11)
         assert mgrid.br[0, 0, 0] == -1.0633399551863771  # -0.9946816978184079
 
+        mgrid = MGrid.from_file(test_file2)
+        assert mgrid.rmin == 1.0
+        assert mgrid.bvec.shape == (10, 12, 4, 3)
+        assert mgrid.bz[1, 1, 1] == -1.012339153040808
+        assert mgrid.ap[1, 1, 1] == -0.3719177477496187
+
+
     def test_add_field_cylinder(self):
         N_points = 5
         br = np.ones(N_points)
diff --git a/tests/test_files/mgrid_ncsx_lowres_test.nc b/tests/test_files/mgrid_ncsx_lowres_test.nc
new file mode 100644
index 0000000000000000000000000000000000000000..ec40923c31e9b6f8f8242d926531bd788e2318c1
GIT binary patch
literal 24008
zcmbWf1z42r*X~UifE}2KsGz8Tfe4DSFfjl{#YQ9y5D^tn6a_^I#RLTrkPfB07v0?{
zjdZ7|>~XH^|9<;nf6w>6`<vq+oO@W#nsu)8{4v}!`}QAXWg!<X#+T`g53akqcQmxG
z-PL$-g?WIN@tIRy<EHXv#$^lhDwe-Lv+&A&-79yrRkW1VZfR&KU%PWlPv<V<`oG6!
z|Kr%-YsvNOcU08>e*gFXVt;Uj`5(;u`@PO}jlb7V`TL1pW-z|aXMAw|KjoQ6{&Ron
z%~;qN-|#U?TxWcprh12QH_4r3|NC#o$5h70|6Z4I-5kbM?B;)utLpquZ$F>h{=xtB
z_N$n;SNR{`pUvd<|Ks}u<o5sL`y>DN_FOlw-O*52zIjWX%=YhN`M*D&e;-T6W2N->
zwlnXj+#{n^)4M}&_xCaV&*kr9%fYy<CF3ib#@}1~^YZsG|69wN@s&;U!9Nd@W6a~W
ze~;hL`Pbw2j3eyYSO1S2{HMMX<M?!KhTzEE{70St+{TS@?UXw=HMIY6l<4Dc4_<$d
zKlu0K0p#aZ?*GfrN06WY*W>ZZ^m13@=G|L5y3EJ*?tiXSW?_M$&-s529DfPDYah+!
zzPSim_1qVXshUE=*&1OH&TUY2AZm9}wE+}A;yIjibOiD~3-Ua0IRF{6dIg=2b3jt;
z_tl<9*&wQEkD+bmQuyuIxxC`oR2clW<AVNzVdyZwYcTO(8#G?KBALXh2UUlU{>l7p
z4MjF4Au<&YAP36#$k)z;)Z}ZT(>H}d;_NrkY<*P_^{7QqeyuC~G1RD8x@H7^9eyGN
z3!g&Q)~Pu+In1Fk=rh;XO*v4}bigcQ-#IAw^i+SJr2%9pROxd*cmPQUU}|;5FvMLJ
z(I4C#4`KBl($UhvP?lGmJVU1j2m4M~3~s4|)^iG1j}^H>!`=ANut#g5^wNcp4GGGS
zZLSkH#YF;=5@er>2b_Y~YX@)dQxSx)?sXlXL#9HYO6Sb^yZ6KW+Y>?-hbpkoAx?0+
zOeRDZXmA~4Pl9UKX0c`Qc97?F$mXHsLr84kNliO01d%cim8Vq+L(u*0nd|OHfRAaz
ziZvJ4gS*|V7NJdY_-sZ%_$SX%OtF_RnzwfhGu_wTUj1bsq*`rmDflS|;hNWvdX#p+
z*H7MWDu<_ovyM`OMSdK75|_1|5SIe0i!ox^hBx5V{tbH8dg(ywbpfe245XeEl6vw;
z>UAQiHv>t%c%=1^)YC*#FBzHr9^Mb6-*hDX#v|$XGm?G{ko3C|Nxwo!`i(@=?>Z#?
z9z{A2AoFNIG7mST^FTTeB=gvTbRJ0Nu?Fcpkj?{%UVb9ciwP3FSRmC461|*3suv`B
zF+`%5F(i5kN1~TNBzk5?qUS**dKN>fXC!*ILZatKNc21(sh*MQ8Ht`<knS6jeOruV
z-<~7cH!&pp_5jJgjUw5%y-4>B$-c28*|#N-AX~gT>b*ZC9brvi^$vqh7D>53o(rJ;
zB9G_q%eGMc<k02!C$2-q{2L(Yz7Z-|4;@*3p&LrR<*GFqtc2_>EWv9z;z47z<ozk-
z3|%ZR&bpg?8{8iywR|4cm!9d-r;wf+3f_%Yj}HjdK;)?pl06OekfNdDu)TEwWQxv<
zIubk&(gqBpf9$k@s6?|*E0$GZXdinwx1t%=@2Y#;`O*loxEB~^O19x2Tkj2g_YAO=
zwcv*_j$-zy<Vv3reY87um$zzaCIkpY37c~SK<I~K;akh*fmi0CyvOg4<JZ!fuN9%e
z*w?sV>4a8ZJeT{;YsNn=;Ljbqmujo-!!B|5D&wkGSh;YcKL4CPOujnH`Oe~_=-iS!
zZ^Tm--j6Qp>$}nhpQ2iL`kMB_?Gf|xk#$ly9zT9UmQ^&KJ>iD_>KDIpl<#B3icEg&
z+{_`;tTc$<ch4HnTGoR(ChHfy${NP-2G)7=r}Ltfki)u}gKy!3A)id5)^G3^xh^JT
zkbr&0N@1&(e#fClJ72df+kkzRCk<qeFy=h(bN5d{b*$NO=KEZqW0*hV`Mzk-#5g&7
z*}AA{=)Cip<-=WP!BMfkX7le+h{+vTm1bFj4N}{x9-O_3?NSw|_!qy$*3s_S&5I6V
z<LZRPd!C)cio5f2wnY8F+znTD`+QfzxD9sN67D<j%htKqJjNJvJ23wk$CAsC{L^l|
z*uGpWyL3DHc-md8H}qLxsTqj%F}jNewpd_|KxlW&VlOQIcI?KGggs_BaLRuOy@Zi_
zmRbi&KSZ;l>6$NcUqMLaMD0C$d&qW|VHxn9jY(BgwzY2Cgr&AerkzQ@hUFy+rO7Mz
zV96_Iet&TV%&Urbd?6r#NkyX<)_EPqfFJYj?aDd_AN?-Io5#dMyu-TMJ16Hs(Zi&q
z?;iUwKul<<rY|4nNKcm(RC$2~>VpFZ4G&`818I(r3uj`Q%?Y9M-ANdI=ln%!-fpzH
zE1X`(ehGryn-UK>azVzx!pc7_*-#}f?$@<@HTsDi-G0VD3sa=a7FCz~W5!!H?Y$i9
zG0phqk|g;iOxUhwt;uPL{#l`l+sfR)G5ffbL2e+!uSM3MkKlyj=g+R)!ZBz%nK|v<
zVsi{%VR@Wm>3vMJv%GO&tQ(Vuf(Fj-QN*Op4<d{oKf<WK(}G7@zvJhd0y|%B7l&|1
zCtc2Q3&`*96x^_ADb(MeNdK6k482D_-rMut6T9ZwhsL5h_Rjy(Rkq3#J9obC&MBIQ
zt*#e>m(TLX#?LH0MQc;BTvq1$N8ii%F;PceRbBy_rssrJY8`-HzX7dcYgx=Son3mJ
z<v6yPwd{Pj(+=AYcl2*q^$>q}dT4~i|HP)9S?eF_Nn+idlGd8Hu~>0EWX9spNm!g)
zy8dj5F4oP;x^pn$Rs77ncW1s<J;T;jGPkxQSYp%N@dXMIrPy?o!+GH2V{A}OKNEdG
z1#3UV@7yr`D^^-O?^C*d49iZjh5N*>!AcIvs$;1~@kd#u+KK#`IN~>dpM=Le9I6_Z
zXO~)xHLF5;geDBJ+VK~+@Xda#TyQvIXzvd!bq-I|dewzRZ=HuEE48ukO2-VTo%gV$
z`igk$UQevQ$E(10rVhI`SwXeseLVYR$@fL)Vz5B<V)5{Y#aO@=bwM#E1M?QQXrAel
z$1I6cT3g#oFl|&jd6``mrU_QPl6X~%Su@w)y*IZF3ngn+7O&L6!XlNFXOA^u)Ghml
z`8mXJF)h~_iSIGYRx#}||4$6@9$&i0KMI33ijH@v9l?Or!8?}mJVpPjA4D9pMKLgM
zd_&RR3JfuDb?kN8g;w8jMnN(o_jgU1C+KAdq@FgAdT~JNNdc+17DzoQAoXSgsn-gm
z-XmbvOZ7xj?+%cDrvT|U2uQy&K>GCu(r*%weyf4>`v*wBra<})M$&IRka=7HGLLc~
z^T-D>k7^+ExC~?-_kqmA3&=d$fy|>on#>~=h+Z;)=;a0wy_^D~7fT>|*$zZ6d_eRf
z1w=0kfaqlp5WTDcqGuBzdcFcg&&@#eEDA)=K|u5z0YuN2f#_LEnyKe%Eg*W<1F~=O
zK=#cRnEO^Z3S{4Gfb3f^kbRQ@vTtvZxo<c6k?tG(-M4EK?LhqCG!Q>%0OAK?K>R=z
zh#&X?(+|Y9BhwFl%|YS^Qb_zj8;Ku00peF?K>W%Ln0_^V0*GJb1M#ctK>X@7GW{y=
z3WK-)`jrqe{VIMI62H;~;-{;C_~~sRe!2~qe!8F&h@bia@zX|N`YBrk5<fkR#82Ch
z_$e<Czvlwt_osmP{dORJF9XEy*8uT*aUgy#j!eIgwL;?e+mZPFF(7%Q1SF5#faH-8
zkUV+^B#+(z$)giM@<<h!d6biZB#%On<k5B@d3gm$UYY^PO9>!(c^pVyRs+e)OF;5+
z1(3W{MUt0(Nb+(RNS@aK$#Vf9c|He7o}UMj=lg)<xf(F@JaQOFp35OK&;P_CGtXZg
z0>T>t(u6naf$+vUAiPlzgg2Ui$s5@-fyo=wMu5p113t*)4eiCq<c$aIKzOK3n($B%
z5FU~T!b287c&HNy52*p+p~cALp{ruZ<e_F3Wb)9~IY4;L6q&rX%NYo-Z3Mz=6F_*)
z9tf}1AmOz#B)oPL39tPGCa(p*1;Ud?$mGefHStWIOwvatPj0x1geN~E;mMgucycuo
zp3FfePcByjCQqKMM#8%eNO<=b65hRrgm>+b@UA)%-W5Q?yS+$wR}z`L%gKv>d3WmO
zBQP>$lU_ub9*mBewXmsq!l<V4rag)8VdTc61sfV;VOYe}tFu=G23y(U^|`jdKm_df
zPF8|`k02*Q9b@S85S^9j=K{T!zV5sIwV^M4lb5w|K8#03k6v2P3_p#381K-#0b`u^
z`FBR_hY{sH^WQivfuYz>2`4Q0VDQL5!ZZmM7~q^Dps6SaeN7p=IX|{SZ^%T@ueDm3
zoO#|tgkMXV&9v{gzjhG(mfwHYZth3;*%jD+^nMtOdZyjfEbxWl*W958!;)bzE7#ye
z5(f;NG`(tm-3I!nw*5+(K9x~#YpeeC7VI%O{x##pcj<)z#}xc#OqFJt{<%f5Y&HD)
z>E`_Fz&#in=U*uhc?d?1H^n~_ya+=(>MSlTwSWQ3msfhUt)YMMGb_y*-=Wvlpm^iK
zBiI(YgX8_QF&y9=)3}*&41V(M@Ge|A8~&J?2B@#z591qzQ?7o|hmoG=i(dvZ>b+V}
zzNe!N20~sh@xLMl{d@tRtSk7T_nqmYQeH=F=hpflTz>+4eW$9Fx4nQ8)q8Th*KWce
z?_UPmzA7-@7?9?vBLJgEN=26Ual=s5Tdf1teK3$=6!lfS3Htf{O79s2K(DE{W$--f
zN%daRdRnxeDy_#)>n)m8uX$2EE86cb+VAoIrr!&b`rY<_^h@SZGie^H=sfOEn#XxM
z5A%PS$6=}$pnCED7rppUy(Ce+Oq-;aY5&m+)pHTmb0gKWCe`!7NqSyM^ejF>^(;&E
z-1Z+m(|yyW`{qIS&FcSh-(2awX<?eyZ4Rvu4F7p@WwS<OB6JPKC%xOd0LJzml3I99
z5{8YIMv1L$gu%7Kc?%RfVSro1ntPiF49xv(nDt5*29@uxc_$kPKaDgC&p*z<IN_GF
zf}y7HDbs!5T7$&+`Q9JB%?hKXS$y(6ElcEJY)RLWg_q{RP+IfxM_YEo;Ko-M6xbfY
zfJ^aV)5Z2Mcu1Ii?duOP5?<kacYY6MaP{snJ2V^fr+*QYSThb|K^tqH)YvihPJShF
zErzil4XtvooM70aT<oL`7Yx>ZFdiPC4})AOtSiNJU~t{3YIlnhFgzngXL{}&EHS*K
zkT~Rl<EwT_TG*_E-);N9wLUuzzZ9b1RqR;`W9KacRBqgX;rbYk8x{LuNGiNgJ531&
z<#XSp4-UZKqK*yKNylMm<x8<M+6q{Z7aQa$R02JxZ_YYXa36k7-=E*dIRInDPacEi
z0~oz+b1kLiEDXnc^)ILqhM~*Tl&V+VgF&4G7a!?Agu$r+@Aia`!r&7QrQQ|((C@nZ
zbWHU<7|KaJ<bUx5j8ti>Tx>rBqbG*%Th}n={B89Ny<Q<0HqJR`^VJE4E?PUJDtp7A
z-umrUui0Rr?Uk}XaWxFG3s@ZP?}Q=ulk@x5nZxjkBTkj23_e`rviZ!yei(Vv5+P{u
z1V+{hM5yQ9gW>&Bvem}^Fm!}drg4)G4Bp`SQQk8O1DU;3w;phU0k>V+%Erzxw0q!*
zafApApEj9@iVcO4d5@yfUV6X?u*KHyy$-|W6&mL(cf;_KjR7t18(~Odz=YG8;s0tq
z#-BnL!9bkYy5fll7`Xef`_w9c;j-1wm4yaiWZLk>-Fh*MI@2<f9w@>HzeA^l+$tFM
zA1V5&k_tl;g1se6cwtB^=X_+BBn;k-R63iK0RuHT#t-zIVc?lomfrp!Fgjb#ZdX-0
zj5u!how>aQMzk+@)w-y`$ch*Rt*|dJ{5h|8v1%<0jlRDWWzTs2RyHSKLggC_hKel^
z+4vp?Pn_JUk(~wu?Q^1s-?qW1t$xex@P{xm^mF1$wmXcxu0Ol-a4C%N%vn6wTNH-h
z^(fottbk$mr=NE}W8BYoQcp}y48zbRKi6RYgD|A~){Zrv@pxLze=E0?7e;?Ot0#=N
z!|0BqO}f58FyiaAZdS`Mqn^*|Ic(cu*mR}z;ZwI^c=36*HK{*fn15&Q{XoX!J<&8;
z_kA~vf(Ta2>cK$#FWCXl7-+AriAiwtgaMCkDGtbok=V=!l8NQ;^NMkKMM5Hs=RUfj
zBP9#Jc7E48SRw=C1=*MC+CyM??-aJ9X^vR={)MX8Y7uZrbXuDn0MJ`8A)j0)1cQZq
zgPV*iU?g}{0b2LM*y+Ih6-G)hR`s*$;W5U)-?{3R^gV<zp`|l6Y{-Lg<FKEamFiF<
zYO+oBL?tBdizq7#-VXyMo%=MzvSDb<b9SWqCKz#^U*%WA=;MCl)&ZVlFluNpJ+2;M
zq)grA!8V3YYeyv`s}{ge%xUREje^isQJQAAQ69?vblo)*-v&cXS6JigI$=2BtanQW
zD~wpo?O9yT0wcxy4}~pz1S8KRi)BnU!mwMYqqs*h49Q4nRh{gIf!WV$3T~N0hvU$>
z{VQCd#4PEZLi`38?)935&x~P2a4=!SkTK`mCZFyM)WgVi70y3h?_s!Q%HY|@i(#m-
zBWUf51u%I2+4f`Cx}cvUNyBi<3+UvS{^iKl4yc>mynD^^l`wKz^mXBBOBhjD_v>gW
zLx-x3-8<$o>WT`g{izj%;ZKKnSvT^*klFJjGdv9gz2~&e!&9Lz;GX%8-s{kH{8N|y
z_jvemE+U0Jw{23rY+CO=t;e|k%zCS6y)(343axj8)=Q=J25G;KX}^p*%zhhbzwc<j
z%O>^PNBezC`W0oT{jR6;=%(|SLgrESi_U}b|1jsVp3WnQ&SQ|y<1L*>9G%Bus+R()
zm(5f!QdBRpR4)vlVCv<?B)y1Gy=YRsY@m7;rFyQKr00v1^sG<y%-9>Io@I!h6;2X8
z2PM&cOQ8E!ME9+d?pq4ow|#Wqw$ptxqx*K0>>FPc**AYijx+tBi~2z@^@CCB2NKi|
zbg3UObkFnyhVL-_U<dVsbJPzWQNJ>ze#P+RfA|%{UzmQyNBt_3`qdEiD@NWi{mO;-
z)shkFr)iV?RGs>%GxgJY>Zh#KPZ@s4^wT%gPZO!17Ebc}IO_MOsNb)pexFADekJw$
zUDWRxzRdJ{d+PUPG>^8@Jc^=u<W2KPlje~W&7&roN1JFKann2sqItA~=A|jk%lVV?
zlHn7~yws<8=|}UDk-yBm{7LhYasQZkX-o55iRQT+&2vSX=higO8UD`9a~GQDUNq0y
zX`WZmJa4AFv1}4=^ibY-NO|J|<qZ+a8}})1FuwkSH@qni$q^p<C_#DX+$0{lOL^$P
zBp&LcJS0nbNRjZGOeNtp{(Y3!zEWPBPkD`@8z!&Kq`bC}@*0CD|G{hXC{L<Vo;*)^
zay#Wo`AIw}MS1ce<w-lrlW!<b_EO$0q`d1xcz4ei%DcNL?>1B3?V!9XIEi<^Qr?Zm
z=HVoc@aXkecx;rj(`FAwZF+jdB62Ez^m6J=smTVtWX0z;N;g2yKzhnAkrnu8+6><J
zTIOi;_~^;y+`i~?ZvX3j#+~TOZsGUlb0_w^SlRK6<1A)NKAYGaq6TN;cGoloJcfV*
z=Z6Ajt03s1{i+l1D&VVT{c}5^!(by={jhT3PJEEEjJE*l(Mjd`x85b&8PCH=MJgU>
zz|s0RuZ5L&AWGZ$Y0-zdkY+!{dX}mvq-pzhv_+>wlDRwU;p{~a^|5iJwuudVtT$%g
zigyFcy{sRH+q>}N5;-p|hr8&#X@d?&^#|;TIscqx!CUFMVQanKw)#We%@cPdIhCR0
z373@ts|e&RIcC4naSbHTXpOP7p8;W7DyimXJ>a~^*Yb#iIwNYW3?6S%MGqf|w69s3
zn5=gsQqFfOG!$$&nEvrN^zJ{>;_&4u)Rl7ozU^fIB~tOf^3)w6qo7IC&Py9&+H&5U
z19tFgw_s~JR)qU6o>G&T;*PH4DYKQuFM(U$F3A%tU66Y>ap5+#htMsjb^lQULW?%*
zwhECoP$9e`RETv0<m_EPjZ1ze#2<V%-B2bPe79+txTT1}o$!ZOr|AWv%Z1C|B+YGU
zy&hW6@L%e!|Cf5(X+6LH)T8}Qq5U44)Gy=t;s4Km3u(U^v|lAUk2n7^k4<zQ{B#~%
zljfmH=W&_NLyPJ~hv>zc@%$!JFCtVg$y6`gR4*)4FLqQfom4N3=ckx@v8H-HN%gEr
z^z7SB^}L+wxs>Wzg6erI(X;e2s%KHE=PJ5yf|K@b8r`>4x^K7WzVXw2dqMZjlkQsw
z-M21m7MOFx<NQwS?LVVuztIcQ{htVJ|HXK|rh}c&RHzy{cs9M&pOOUiqRDTTZ+s2~
zI=L+FuY17h<n-?Rw_TXh?#VW+CxaEomR+4O?tnk^4O-#i8|>Zs;)T#D&3K;GiR;>{
z1K_6`mmBY|Q_#`hto|l^Bh+uTzxqQ$7>WuOI-X(tevrU*%e*{p6Bz2G3p|&zMDIt3
zje>QWuu~mQPyaIo2i{%u*|qo%RP*@^R?iTEzTd~^?o@pVKXy)e6Ze1@zAFoDIydnP
zij-DGEUV0i)JEHI>9ms&<t-HGCL+!_za;0Z+2$44e@^9YEL@DA<$kO6iuMFF7FVR1
z*dKu=j#>rd1Kv<4R$6qfSr^K_CiPrd83*~k;x}ua+=H~kvWIJTvp}NoJ|1WGF-Tb;
ztHUm6f#XxSwqNnvk3%dm3MWq^>YFMZ&i3|#{EyMIK8UtKVc4tVdqi6xfBG{uwoBU~
z<9Hrt$$KS8=6KU1Yj7447n*usTw(>O++Ss;d@YJ+JsqHwdX5Ks+BTimvE6`0!Y<tk
zUo7D6bibRc8=gU|OsoHx)+R`}$<Zm~vjyU1SKqTOdJD0>26HymR70G9*oLS3#2|UK
zuKWpA5$sowEYI52jBRVnV)-qDv7%@fub`JICdqrV-QaqTR%vTaD!zRTt}HpPr^J?k
zZ(sJ7xtrcY(7b8FM?3r>Tz~B0aw}&@NN9V<D!mhXPNy36>>tGzO|j(@iv+QHg->05
zI}hdtPwxzulET=N8V+UyrRb5hkY_;N25(PXR$S4h51-fb*L!cOhkyW!l!m}55N~wk
zy8ll#Y?oE#be|G|O&qsgJ8pc5)glJ`>o4EILbE+bLt4LJnvSi@oVA@8Q&2BHtu6@t
z-!|qy)P0Mu*yqeMVLVUnd8A5sx@Zt2IF#hDh^k=I)M-=Y?4z)L-@dK92Q0BlaH_@4
zrZFtOXQsQO&<=C;@4SxB@WWK)Q{KBC$77se;X)3VHTbpdLh%cyF>u|sL{{;hBcz=j
zWSw~3i{B>>CFl2l$C^goDSoqpu~K051x-64EL~HbDiAvd3;j!r=he@`yty{#-Sm!P
znjWjj&Z1I`m~rK)kyAK)J=4T@utXAy8h#y~s&xn}S&r(7-eKIYSzc3l66~;|MnWp<
z?piEAz&U&H&@h(vKTF(ixeiP8vLiJ=yu?Cb<0+pn$YHV^Ys40ly?_hf^Lo@aL-V`l
z8;|s5@kjqz$%Wf_vHn&_CA(!Smf7iaG;#W2mU4AQ%svkERXirMo_88vUAwCO*++RS
z@Y?X4_gf40#fVwUoaca2!?~q54;Di0uS=CYf-2bN?YDA41uHfT;>)7jAFy2a7hAD0
z3uf4cT&?-4jo;+1oZZM73hzw(y&Qhe!l2ybyAv<WvEj(hC>KjrY4#h+>JNV&hE^fT
zgou^p*q&I$oo{K0jdwb>v^f@FdHur0(^l|eM&OqkHIW|tx_{0&t7=K`GL&tb*PjKU
zUWx7-5~LuyXu41BHbrRqU90v@FaetPtmo$axB%M|+&@^n8p1~1Wht!o_ppNR@CL89
zDwrvGWtsAW3G}{Tcs=*{CJ1)xHql8}f(+I$$vo}%P;pyYdHF6G=vtgnf6q7@`Z|}c
z-Crbt?e{9X`o4a{#!Kxjd|_@_A?u#*tUZL83C@E?-)^8+|3L4J`IQi{cANj0Of%$}
zHC%bIHyggMmKBn8`wd-6oHyrGTftb(8@IR8f3RJpvffBl0vpfrO3G<6>Rn$a+?XeX
zS@U*mk!8P)Ufr_-BR6k{$hwb1%V%za{P#1j+&bX^)iL3XC&$-A(~{D8!+z;d<}U5w
z-1MH-+wm{;?$LUOXuTD*9^<?nX1!IkUIwlAiT3MB`*orH_R)SB_wR4N3s%y8Q)s{b
zwBG^RuOsPqj}o0nC!L2boyP?-kB77AJg(DuET{8eoKM7@hcKN-2Au~x)r&FJOV%X4
z{HA(|r+QgQ^b+q$^dh>E>LrZoWh2pZ`5@J^CDF6d1l6+%(KFX&qUWEzM9(jeQ$0Hp
zJ&&)XdLE<uHjnOGBi*;Jbl>>MzO9@~_ATf-**BlXWZ%w(kbOJg^<Vo&{J@m@ferP8
zHPjCVsUNJReqc`hK$`f$(|+Oy_0y;yY@~iwP5mm6`qc*NS2L+!aZK{7H0oCx#IH)Z
zsb9%bKaHe*`jYzT5cN}A>ZjkRpGH$ZolgBUg!t)w4&tY881rTNy&v^^vq^sMNd3Np
z`h6(%dsFK7d#K+xQ@_`sc{EJ(XbsJyWi*d=(LAcBdE`R#NSx+T1I?ouG><ORy!=7)
zGKl7-Ce6ziG%urRUdGbAd`<H*h~}j$&CAU+&pBzHuc3MVndEt$7R~cnG|xq8p1-Df
zUP<$OBhB;GlsC3f-uOazBiMlQMk(cuYRVg;lsC>&-e8<(z~qfrlsDog@sJ(mp*@s`
z>?jZEP##LAJanJ(5Ig0e0m?(SD6dVYy!MXr+6~HUPyfMd=ULb(uT@fBQ=+`Knef`C
zSjv+xDNkxsp46l~d7kp5I1ru`r92r*c~Xk-q^tnp$vaCa?{1>JYejf>T%Pi77v<ez
z%DYyScYhMzeRrPn?s3YyIruUzT1JIy7`^j%e~X_!jA7h;CqgGaV;r0NxWKM9Oq7_J
zBW@{*$w6utQr>uCs^$G_vK!xE+7FZMDP0FJy;|%T&w5kLxV=&C)K~?0cpHhT@&FiW
zr#T7jWI<QACthNj`53H~kmxAKgHfp;++2hLF`i3b_^9eROlak-HCf+*NunuV)vDAm
zd43Vw#rd6>Y;BNJbYcL)RqCC%AB#cI;K*wmzd?AKveadfp%l6bH$G8v55d6pD9%9^
zJ`7ji^vFgYF=|L+QCCtX#vWa7k-w!G<1#iX7A-PIwA;=5b<-h;VEM4-WGgqM`ia-~
zEbWJH^KE~dvFM<F?lv<G+iB>Q_3%{remV4C^YZPUxJnH2IP}L!SO-HaeWYGz7-4AW
ztc2<%moQ9TR&%FlCj``G%zMDi1F@?t2OoP1L-4(oBQp%Izze~>BX=LKL91=`*Uz>j
zptJFf?1Ee?bkDvmFrIh-y^d!es(kwmy|biwO`=%vt7&-2#G`7k^%UjoG&u~eeC1nt
zOH;sgW8urapZ(zJp3vXjsx_#q(LAshq|sE}=J(s7wP@+QPF+rt72k3XJPpp6g722q
z)%y7j;k!@b3xdP?kks3Yq+T+TdOk?%RRgJa1ZX`V^;m$^V+UFf{!@?iy9`Oc8<F;l
zq~B#g`i%$D?|UHq?g!d0(0+mTi*z1H<}rYD9zf<X8|XZM&I9N?fX)NxJb=uj3yEHY
zk?5rxiC(q<(MvcGy$Ar&OA-*hWB}1i3lP0n0@2H1BzhJ>rk-~^0HWtAAbMU2M9-l>
z^n3w`p1XnQSsjR;yMX9<5Xinw0NJ-8p!)`7-&%p}+g2d^<^g2iUIW>;{Xq8ZIgouT
z$BdCMdGQ&aF?YpW|IvO`%;&KdnHAxJ+10MEB=+sYO!i0nzLjue`c>ZpDzmvU?Fq*V
z{%35MwyER(I{6)#=JKs~g=+(5B()YAZK%Vnjbb%-bhlzo#GVI(^@p%P>*}Sx^n;jX
zF8uq=(|k;e-<2c&h7D63X6{>>W`N13m;2^_)yJfuRX6xwhhtK_`mt6A7fc?M&_wkQ
zm_4&jDnfG<a}4Fa=C#RU_C*!bX5M;C=X!s<cR%Aiqvw8StcMLTaf<`^=GkDp?}{s(
z2k#=jw^_ya<38ea`I>iTcQJl`(LldEKj!#e635R+F>mI8mLaPmrq`J%Z>w8}$&0Lb
z=5&9?c%Omoi3%+kV;r>Vs-G`L3dAV22+Y8+U%hX&#eZUGm)F6nV|y`dH4E>oE>Fxi
z9&$<#dWhL7qt}%LLNK}8dd^4Z?TC>R={z3}VniYy_~Ph>!7-2Zwr#$H0RoQ{OI}*w
zw<~Plcc!YN*L}By+rr1t^J#nV#5)PhpBreC%_oMb@yDuEJ~tzlOmq!+9mB94HYy5j
zyy*ARQ=she6m*k2{c5Z4HvHr^_1Q5QD||iwH?PabXf&RgY8$R_9d*nCm!!+uV@mdm
z+Q8u!#QA}<-Qv<QWShO^%4liyT$pF(k?(_6s&_8E9rplv@iVpkIve2aGUZRYA`jv7
zN}n&kzg!3B8<vOeSe^wp7A^DNLN_sLUrI#!xl<VM+~l;X#s_qq+d4zS;}9J6<tfiz
zH5)#SnVvkfp&dN(S?oeAZh_z76y?*mc0pj9SR9Ya7YJ_O$g=*)I|x->>9YFJc67B9
z&G{+fh!0;2s44Nk1)ISV_3ESw@QxefurrN<psyEmOt<<$_{jznp@9mB>R;-hDHsT`
zUyq4<&YuPtz+TDQmH_ci*X)CKorBLRpNu{p>V>bp8~mmwwnGT-3gN;<sSpMBvDx3h
z1HPLkT_Kna3FE$PdG8o|7S;c4>zV{enO}G7`;r%s%43|SaH0=Vdy>8gbbCO!kmO#8
z>qQV#>?5)A{RAW^KG$oS!2`*QA9!_-voZE7<8ye_6G->)ioG(g88X$3vyAHILRQX#
zmWNOHA$u$<{ie7N<dlmpefwb<B)u|PZ+uG!Qu$^yw_Xi_^rSGgkOnTuI@A#}XqOE+
zAHU}WS09BuZ*z$dqYsdu$}u*4eHaS7Ocg70K0x84*!pqVR47bGj?u10%um>rl&g9c
z3m>ICkdENPqK-_q-5GPR_)nQZq0KPkJaU#_wnN9U#A%c4CH<pVa$auX*}#`r?6vpJ
zrFEKEaI(nEpYH;u2Pa&aSKGmO{l)Q)NT=7BCzTU*#z+_QHTi<w>)tTVw^TprbLtot
z2yE}w{>^~}QzDC<Gb1tICq-;~1UKg2zWk(HqXqM&K0PnZI*X|nTc3P6(TQojU6}&`
z{FvV0vP8Oi6f<0GdF~5(VCDk;tv}CaVkXZkZj-Ugn1TF8DRLQ@vDm&Vk8z%Ty1bK-
z8A}4D+Qsa<7N3rZcgrP?%xA>}k2ZD@Z@@&3o^G#OCYX3;%I$)IQ<!-FmY9FAJ|<qO
z-0vm$0~5K<c%L<?#RT){=P$hrF}^v#>$q$LV$qd78kR>e=1}N!Zk1aYYO;iNpZOCE
z%@H`m{$>k?u?vldnk>UGt_udgF&0Bh=Ij@H)rO(!;&D!gtuZ7e{Kx!SR}A4=u%-C9
zE(VJlUTO5YfI(LuznZ#>5A7?q{P?2u4IS2gFZuHA4mzA&la#vi0XoPvn}nuwp#7L&
z?5(*M(f;DrEA=Z>@Uwa5%QtyLXlJ`XdSz2Ce)_<B&6}eZZKC$fRpmSeq@FU6ddGp(
zTMf*5X+c2h?Eq454KVBFF9lLB9GLYu1%cE%3Z&l~K>A$<q~EDP`V9lpZzGU?U4ZoK
z0Hoi0K>EE3q~Faz=CK{fJamD~;{=d-+y*j_7NGM0GLH-(^GE?Qk8?og@fe6+76Q>r
zHxRwJ0MSbl5WOe^(MvKAy(9zG3lP2J0MW}_AbNfdM9(*X=-Hg{dV;@tmIb2csX+7`
z0YuM>fatjqh@Ni((en_HeN$lU;op5*0c7810ok`9#(BJd_e}=KzA?Ts_w5moeY*%`
z-%bG8w>BhxV2H#I;*t2l5+r_5fW!|@BJl%9Bz`c)_`UkCAG~Mm{a-&wK;j2wNc^f2
znSNz8ABkV_A@Qpz4F3M>S0+gOsuzi0Ek)v2>yh|X0205NiA+B|yakD${z2lW`AGbf
zkxxuN%|YU)ok;vN35lOtBJtB+B!21$#P6R2)9=?@K&Ia(^CR*56r_HS#P5F~@%tzw
zem?_=-^U~Id&cie%skoyB#$hCnMdy5k>t^9WaiOERwQ|3jwFvxAju;KBzg1-NgnM7
zl9#PO@-h%eUP=JT%NQVenSji^)Y3+hmv4~dWgwEgR6~-N_kra3S|E8|2_(<ufaLjR
zVCMOwCLnpf2}qtF2a@M;K=OPYF!Nka7zl6N2f`a(KzL&-5Z*WglsACzh7J(kcmjks
zE&<_<qd<6LKM)=|284%Xf$)$B5FT;_!b1Z<cqkSK4{-pKha?!!KQVbo0tgSe0h8Bm
zOasDe3P5;m2neqs5MDC^!fQ@Ic+CI^ubBhkH4PxVwgCuF?gPS;jD2PDqyiA0ECtGw
zKzQ;M5T1MjgeMb#@Fe4LV)CRV5Z?6w!n+rN@U90C-aQS3cUJ@9-DN;{mlFu@3IXBW
zIv~9J9pg5qhIx-|#RS%0<;GhsV0_Dmfd|Dm5Y47^vB;ZYj5%*X>gnwmX&bzrZ{$6O
zFDX`I`|g6FBIx$}<}M6w9hfy7T#muVg_6{;15+1@ep&KG3DZ9M#Js)~i7Dfnl3Kb!
zm}GymQFYlzOjuJg-CbT0W83<tyj-^&qcj{%^YV9Nc>cY_`@CZqw(ykTLEk70-LTg=
zzSRVC4tpJ1qBDf~yv9d!)9o?qR{!3mdQUO!;q~r&u0JrTq1ogiO9{sJSZuc8;>Xx?
zW9gn7f-!2L<mopR%Q3tus?C&bJBE4JDdcDrWA2+{%lk%6vGn?^mPCW?n6Kviz$CW_
zvv@)>f^;J=^|-WeS-c`9e%4=hb}$w(!T!!0X;F-B_<5_^J|82GO%WIE*TwMc5f9O`
zwwT?NvXWcP8jBW`3l;5Mh=r=1i$mrpVD17jz07D2Oc(MtdxG+qR3g5)L%j#%V@?|P
zo|nPcwGQV`&k)3@Lwtix@7H3)+xCI-t45f;tgZJ_>j^9f`8`Kws~Q%L3C0^P6~w%p
zKrs>X0L&1zD+#U{#AHK>tz1)uF+t?VGw$<k7@KpQU$aILqte#vMwY(8h=C{by5`QJ
z^%&=QGV7`S59_tldZ#AUn?w6$yx#Tys^5dO-}SU#f&c3F|9u{X|IIwa=sfcNk9qu4
zFLqQfM*odoLa1JNsa`~>UP`E*&rv;xQ$71n((@^*=OvT$JVNyx{a<<}`?h(~zG?o;
zzFE<IbN-io>%_cA1Nxz&1z3_<preq<f)y7x=N3t7W0Ax@3!(JYm_871x<FD2QN?(4
zLE(7}_TiDgA>)ERHx#2d-5;aR&Ldm-YCd65#-c?A#*Z++!1rdiC=Zt97%wT0ZNi2+
z7eApX4p@?>4&rShnAPCfs;11&;Mu#ce4mS8#E9m>Hk;QNB)Dgx`@EItcfaS8m-`3w
zP1@N~=_Sh8+qiwAyAES{f52*kIqFzdVx(7h_7N8EnZtADn=NMH$oQ?7cQHw2AdFk;
z9!7u6T>bT89EO}fzjB?EBl`P4`_c0{0KLT}RbQ@(!%QFfi@8>du|jcb46pMOEM?u@
zz3!YD7OZ+DJ$vQ}%y@V3mZhv0CXT-PYEj;YF=trsM15L_p?bk0&#!l2fJ$`*7wcv8
z{I*isN!uCYm4=;)gMBbpUF2ef!bi*-Ol-Nb{tsrw_iniF{RvZzW6$1VJdc-<5v%(<
zVL3(*F7??J-G`w=AB|KVPhcRITNy`UCwkq}`<1o11-(aRe6)J-023oma1`DW$JG3S
zU4rkeF?H_APo`Clm~>e}v2kN5#-E7UY&-J|Msroi3cYZ`u&oXQuBSF*U~PW%DFbf&
zT9Gho(*bF;*F1Q8#Z(!LNICY}dp}@&oYMWrg&Q!zlSj=l#RcQP3iadYbBwL*cr~xa
zfT3fk9X681Fdg4C-Zv{TsO`k;$X{&exAUv~^a@0CBN6kZiGCRPWYLO*;xvp5RdW0C
z)&^tQ9?YFHEWy~*oppGw1fyMAn*Aggv8>T6HNHRs!#a6w8UrmbMC!WJuto<4Zs}oL
zJlu^Y)Ao3(ENR1Ux=szsAq>76d$q;jIOF~nhiH^YFwQ4BpxA43w+y4i?;V+;dJiMz
zN0ZK9bYQ%W`a(eG*5w#Bv))N`+HwqztG;=*V*uaUZkJp0jTgQ2whR1BXvCl+KiZat
zjAPjTrL#UP>%)k5Uj+j=voSIt?2035VdT&4xiho`F!Ix-h`BPd81XqdQ1^m8<N1#K
zbJHJkpzZO1#SYto(bJ@3dc3Iy2IS{#lsRLAA$ex1;h9MoUL3Z;#a#lU7S#UmxvGrO
zuG)|LV(wv#j#+SLVFSio4DPG&bHV7ttX5Uc{b(=r>fz@H0_c%DPkH=iJNifS=Nw(d
zgTcS6rpw&Z$MB3p@44!vF#4>CS<$N>7-wBna<D21<1M6imKYtwgk{Qp>d5$g?q!Je
z1(hF&!<*ua=d@rln0ZPr(#MprRV?i(DVQQWO~~|&4<_dESLQ@Ap2t<U?VaAd5F;))
z@q1<0V|3m+LH{bo-)kz_9=*rY=q>W_yqNrahF|zN_|NJ?49@T^D^kFOw4r(><p50h
z#w~O)sU5LDc6HD-NsOv!zPipv2E(dwp3L#r7^TK@BqDJ$rhY6qr=)lVoldQP>c`=W
zzGuSUCSUc!@Up#6S8U#mF|W*R3RiE&SlM!4*8?pW&F&NJd@2gV8eP?=H1%MRp@X5C
zn;!bJ$;!`M+k>vkE221DEx_@%>9+eBndn%s!?7yX7y}Jo&MH&4WxRibiu#X)HjFSj
zaHUYJ0>j4Q*dISUi@{qf0z@sM@tcYF9EZ-Gj5)3nttyxX&dYc0IB#nJiM}Pv2Hq5)
zsbFDftvrM87drjwky6LtsR`elWJfVnaG?uNj~a$7OG+`Fau5S<zk%`M-RS0~{?10;
z5jCuMXUNUm3jQZDAEcfVg7nZGT+hw&@R0hEYnCY*=x)q0B9K*yfm4sj4wa=~@NK(}
zIjJ`>$RY6JZz(VIV;TM|F1Q1ok2a>Q${hwB`^OgLt$V;%#cP!-`wobHsS+mK+)V1_
zOsDngXuVik?>epbmDUrW^*U)iR$6Z+tyfR`<)i(cr2Y2Oe$UZ<8Ly*f_A5vG%_jZM
zzD4_;NBga#^EgT8VNK`JLg&#!=W&A0BZ1CCn$BY-orfZwM>5q5V-J{m*+%ukL-kTi
z^@3C{Ra7sAL@(<Ye=qr~7dJMl7ZIxGCaPzkdUmIJE~k3lI!VvgRL^o$&;3NtoLqF@
zGU&c((tV4d`{qLTt%mN~O1f`7bl-gFzMZ7|wx9ZeKJ|kk>IaO!3o!lQF7<=G)DL<l
z`N2u*2l><w>ZxB9Qok~we&tF1DvJ6QP`|oN{puC<t9I&FZq%<jsGlyPe%eC(G?b0{
zX%F?&T<WI}sGqK-etMkxsoj71De?O;;`gERsNWx;et(eqJ>zwoOuv6j{oaxK{R-;$
zYpCCk&^#KXd6Yr(h>>Ui$fF9HM@MNMF?gPtM}{<y`e+`tlf0b5IDdzkmwPAWr6bMD
zXEZMv@7KW0OQ3nVjO69`1kKA-n&)?Eo^PdjZb9?>1I_asn&*l%&p*>VUq|zN5zTWs
z${S9UH!LV`bWz@TO?YF+DCG@a${XJ)Z;VsksHeP9OL=HM<sog#Lr*9V9i=?vPI+j6
z@(?HGp%%(RGL(nxDX+;;Ub{hgO^@>02g++&l-KT4Uh|>6_JQ)6JLNUT^F2(STsDa(
zZ75G}AUt_#g7PHe_W~wQKBYX#M|tuF<;f+KcZ&$`X8feQt511%Gv(d=ly{XV?=s$Z
zh{?OUly@2RnY^nGZQ)`I1=5<J#r#gr#Nlbs$az}%>fF8Xy|en0zoRvj^>DT!{|FSE
zxb&k;k@0tKmH0WvMXC_CtHIDcS{fdFT(oUn`zH+O3bt@6mWR%1G9wuRN1%21i@pKB
zAT&?gGBU<^ol`yQw*#J2^Pp1lxl{Q(UnpWbVsQG@3dZw<S}X&uQ4lrT`qZ!CW_VxB
zIdao!1^Ps)EnTDj8oDI}zTOjRg--a{mJ=%pKQ5gX_|3-;4H2Hpg3Lvs+F_vc&|MD3
z^K)|5RyQ>uYe(hzu8i*xYyNAcvauyN<gzws=8L1(tUWiAc4R=elf;+14*Q@t&u>Aj
zTNJb}iC#Y?at$<D+o>K7G>7j)Q)f=wqza|$wv`@VnFKjf$5sWcxeW0V9J8WrZh-6c
zy(fn~&!d}!_1qnT^PwBsZi%0BgPxwH2ZiT|KxZk3XZyBdXt~kWvNmEC)EQO=xGqp&
zyx+s+%oRK_kc%%qS1R;F!u~ADknuF|u=k%U|0@-ppI!?v%r1rQMfP`^8|FdJ3d47)
z{)eD@hvJ8$@qX|l;``X|jX$A&*>e*q)><h0ZN6N?^*ZFenc;68w+j*qPst@I+y^h2
zhl%T4h0$qQKmUCTqyL}vBxyZ^N%bDle%=42-<h=E<Nr;+iL~D`I*+|{9?_HLVL|7i
zLg$e@X&#z%9vA35xF+dkH`U9UNqU(!NiQ{%^s<BMC7<e<aXuhZ&ugikqp6;A|3%Lc
zRL@OR&r7MEpH14g>2%-rP1?87N&BWYY2O^^zR5zn*g20cHX_h|-H+#9*am1jv?J0y
z?+!GrJ@9+seNOnkXXD{e?uAgwYyI6vVhQ95J6ib6<AQ|wJJ$2pj=+~g+%*;twJ_|7
zyr8X8CUneley~R32Xt(IZU4F?9=a6rYM&ZeLsRh%q5RS?`0lHpR^@+-F^^~goAuiv
zC-HnwR?!)Ve>`~jmC01FuiW!!nF?ZvE>EwpmnL*rF1P1&J_wzMWv?qv_(8jpky1Cy
z2sB-_i&gAC4BzE?>;k^sgOXh7us_GvLiTt_v3=wpz(;%;?3RM?X~S9DDHoG4*l5em
z&mV6?hwY`qdy0ag%k}pO*R}{~buO5vVJ!`f7mmE-E<OY`(lx8kMW#XtcaGGH@jS@N
zTIVG<>l(xwSPj^XKL#s{!S}+UyD`{O>Uwv#6?Ck6e(hY?IcPsy@$^+e6~hMy+m*%a
zpkYOJ&C!o%penh1!)LAoP&`k<MEckWWJ;LdJ*yYb_&cV)PTfZY^Sx`dY*Pj><goCB
ztw0R44VPrTSyl=^nm2p=+Qtt}u4mE$n%+ZQ<%{VJPi&!bZsFH$?kAwIY?QrLt`^cy
z?`%ytKM3I?8#s1>2WVJbIZ)hTkD;dbZ*0EJ0j+9mC%b=qhn6KPPRa6jL*t^^`3ldB
zp!S5zrL1SFP<~`<;OQx|pumXLVX6CTNR{|-$n1M5gxE&q#ZLD|L(kr|i#1&^d}m%&
ztWGZc5R0vi+qe*#qL-u=-#Z5N8F4Sf_vk^*N`;GBvBFULT&wol;vmT9MwU#|*O1h=
zb#0Q>LkMV!d&{@j5Utp24++=AW8|Gye+At|&>|z*&|}MZf00l70tIe=hB{C6x4x&`
zpo+UK$uwvl6laO1N7)>NoTL2}t~+)@!mF#|mu6RjcgE+YmaDbsB>ifQ+C%`x91si^
z3LJyx+i3&8!$P3pW8wj)pjS{UJ?p@GJuk-VK6LNB3QmBc6;>w~o)>{k&2>BZ)_Fpl
zkZOHiTO&B{cxPD1BZi*CPp_@(;zh=Zg1+JwjPv51cDAH<G0y*YU$)>=%q6HfS*=)n
zz!AzW#x>{hwL`%}cE{HvYas1~xAFEj`Vh&tDTEXDz}p!!JvP1jh5;rYcWK(qz(hgr
zs=3~_(3q|6y4|)4>Y|%<6%U?+s)vOlH4izV#4ZG86);|}5kEuu_-@8|AAT~zHRV$w
zc<8De`}%WuTUYyNlKmA7JK`}@T89mj8-^~HXDowe_ZGV`{v*(+-fZM_{WjEHw^{0-
zv;!&+ZJXaP*8z$woAp{*jzLD5_U2xzHxRX-J8-RlA{bxa-NN#c6P*)u+QZwH<FiTy
zDR>+T&DGyyT(+q|<EeMDT95BS-Q_6<tE<hR^5z9I{Vv8lxD?V^P2WM5uf2ZRqSt_*
zE5A1u3WCGaF+OX{Objest32SrkJc$lA?m`((EQW#Ug|9_Xq4Kozs!vT>W<9NpHN;5
zl{#XgA+aV<tY8q;cU&B@%e$oBOU{9W*lL@a`RU-B!`hd#(*jS<jembodm36Sz3o`(
zU<EBxtLJ{*V*!nV(RX}r?Sr~)3)k&08Gy>GPgjmyUIWDjc)!Zk267zr5{%rqA!&v6
z){Lp65FA%q&6@oMK0OY;-f^-8%nQtun)vvjWqO0duZy#xaY<exw@DJzt?Ac&vZ@v;
zFS33tm>vVgZ+-+{<G&BNOV8=ea#;yUzi*WeIX{3f%N3nUrW4@$jO)phjXhwSq_tS&
zf;yv~USrKeKWLn5UFP~>Kh!PvaJaTP1S%El^Fml{p!ie6)GAXK$UPNdGB@WKBx`b7
zOP$Vu@RD<p2O=`Tz5bTJnnfI#DKGk#@jGl%y#usf7_Dbb>$%Z-hiSdcN%g+bdWy8)
zS+w5-+HVc*H;VRqkM_Hi_Nzzx?Wg@Fl77$E(|MTEdF-Y0=%e%Cqw}z$^O#BJ(MRX;
zlFkFjJa|Q@UO1^<rc%8upn5T)df7_#!cO&)K=mR=^zvvC(aRS$s%IstXC12NXH?H~
zsGcpUo}W-X^H4p%pnCQudKNiL_w69vx0`g|#OS`+(tVSm`({e_?FiX7%gtopTKmbq
z86Bp6P)+@Skrzxqm`eRXnEHV-^#gY52MW{=+=(CD>L-4XK0L{<l&D{Qp?>v)`qgUc
zSGv@%B&lE7Q@>)Pe)Wv{X)yKE*VIp4sh>tuKNX~Y`jz@=4fRt|>ZgmTpNdew-|;Vg
zKa=`B<L`Y;zu!&$o^d|VKm0z7`hDA^Jd&e%q)YQ?Bh4d!nn!PF9^IpPWK8p@i008;
znwQlyFFR;nKB9T4MDvnkQeGC)ymTjdd1D>T%XFIOMKsU(Xr8~Od2UVfT%P9ncbez@
zG|xZNJlCap&O&)Zk@7|%<qaFk8>c95?4i8zlk!H)B;LrOyb(2thx{oI0p+3nl!uZh
z@zB*tJj6$N=o95N8Om#cl-GtP@fzcGolIU^M|tfh<u&0+ympK7WGCfG#{1JSdGa;o
z$?!=$xr6egBIQXA%99@`Pijx%T_MW5cPa03P~PRHyvs{@myt_M-u**)m+|~KlXt`A
t8E+sBy^K%f^7sAG|GqE!-}kTn`<~mZ|2c+rjDImM(#+3{pQ$6?{4Xf<v&aAd

literal 0
HcmV?d00001


From 377a793581bf8acfd0c54d78b5e3a39348d82f42 Mon Sep 17 00:00:00 2001
From: Matt Landreman <mattland@umd.edu>
Date: Sun, 7 Apr 2024 08:58:44 -0400
Subject: [PATCH 8/8] minor tweaks to mgrid test and docstrings

---
 src/simsopt/field/magneticfield.py |  4 ++--
 src/simsopt/field/mgrid.py         | 14 +++++++-------
 tests/field/test_magneticfields.py | 29 ++++++++++++++++++-----------
 tests/field/test_mgrid.py          |  2 +-
 4 files changed, 28 insertions(+), 21 deletions(-)

diff --git a/src/simsopt/field/magneticfield.py b/src/simsopt/field/magneticfield.py
index 9897244a2..8a647633c 100644
--- a/src/simsopt/field/magneticfield.py
+++ b/src/simsopt/field/magneticfield.py
@@ -112,7 +112,7 @@ def to_mgrid(self, filename, nr=10, nphi=4, nz=12, rmin=1.0, rmax=2.0, zmin=-0.5
             zmin: Minimum value of z for the grid.
             zmax: Maximum value of z for the grid.
             nfp: Number of field periods.
-            include_potential: Boolean to include the magnetic potential, defaults to false.
+            include_potential: Boolean to include the vector potential A. Defaults to false.
         """
 
         rs = np.linspace(rmin, rmax, nr, endpoint=True)
@@ -138,7 +138,7 @@ def to_mgrid(self, filename, nr=10, nphi=4, nz=12, rmin=1.0, rmax=2.0, zmin=-0.5
         
         if include_potential:
             A = self.A_cyl()
-            # shape the Potential components
+            # shape the potential components
             ar, ap, az = A.T
             ar_3 = ar.reshape((nphi, nz, nr))
             ap_3 = ap.reshape((nphi, nz, nr))
diff --git a/src/simsopt/field/mgrid.py b/src/simsopt/field/mgrid.py
index 01c6ce34a..8e7528cd2 100644
--- a/src/simsopt/field/mgrid.py
+++ b/src/simsopt/field/mgrid.py
@@ -75,15 +75,15 @@ def __init__(self,  # fname='temp', #binary=False,
 
     def add_field_cylindrical(self, br, bp, bz, ar=None, ap=None, az=None, name=None):
         '''
-        This function saves the vector field B.
-        B is defined by cylindrical components.
+        This function saves the vector field B, and (optionally) the vector potential A, to the Mgrid object.
+        B and A are provided on a tensor product grid in cylindrical components.
 
         The Mgrid array assumes B is sampled linearly first in r, then z, and last phi.
         Python arrays use the opposite convention such that B[0] gives a (r,z) square at const phi
         and B[0,0] gives a radial line and const phi and z.
 
-        It is assumed that the (br,bp,bz) inputs for this function is already in a
-        (nphi, nz, nr) shaped array.
+        It is assumed that each of the inputs ``br``, ``bp``, and ``bz`` for this function are already
+        arrays of shape ``(nphi, nz, nr)``. The same is true for ``ar``, ``ap``, and ``az`` if they are provided.
 
         This function may be called once for each coil group, 
         to save sets of fields that can be scaled using EXTCUR in VMEC.
@@ -92,9 +92,9 @@ def add_field_cylindrical(self, br, bp, bz, ar=None, ap=None, az=None, name=None
             br: the radial component of B field
             bp: the azimuthal component of B field
             bz: the axial component of B field
-            ar: the radial component of the magnetic potential
-            ap: the azimuthal component of the magnetic potential
-            az: the axial component of the magnetic potential
+            ar: the radial component of the vector potential A
+            ap: the azimuthal component of the vector potential A
+            az: the axial component of the vector potential A
             name: Name of the coil group
         '''
 
diff --git a/tests/field/test_magneticfields.py b/tests/field/test_magneticfields.py
index 965ab01bb..ed13a4d60 100644
--- a/tests/field/test_magneticfields.py
+++ b/tests/field/test_magneticfields.py
@@ -1041,14 +1041,14 @@ def test_to_vtk(self):
         bs = BiotSavart(coils)
         bs.to_vtk('/tmp/bfield')
 
-    def test_to_mgrid(self, include_potential=True):
+    def subtest_to_mgrid(self, include_potential):
         curves, currents, ma = get_ncsx_data()
         nfp = 3
         coils = coils_via_symmetries(curves, currents, nfp, True)
         bs = BiotSavart(coils)
         with tempfile.TemporaryDirectory() as tmpdir:
             filename = Path(tmpdir) / "mgrid.bfield.nc"
-            bs.to_mgrid(filename, nfp=nfp, include_potential=True)
+            bs.to_mgrid(filename, nfp=nfp, include_potential=include_potential)
 
             # Compare the B data in the file to a separate evaluation here
             with netcdf_file(filename, mmap=False) as f:
@@ -1062,18 +1062,19 @@ def test_to_mgrid(self, include_potential=True):
                 Br = f.variables["br_001"][()]
                 Bphi = f.variables["bp_001"][()]
                 Bz = f.variables["bz_001"][()]
-                Ar = f.variables["ar_001"][()]
-                Aphi = f.variables["ap_001"][()]
-                Az = f.variables["az_001"][()]
                 assert nr == f.dimensions["rad"]
                 assert nphi == f.dimensions["phi"]
                 assert nz == f.dimensions["zee"]
                 assert Br.shape == (nphi, nz, nr)
                 assert Bphi.shape == (nphi, nz, nr)
                 assert Bz.shape == (nphi, nz, nr)
-                assert Ar.shape == (nphi, nz, nr)
-                assert Aphi.shape == (nphi, nz, nr)
-                assert Az.shape == (nphi, nz, nr)
+                if include_potential:
+                    Ar = f.variables["ar_001"][()]
+                    Aphi = f.variables["ap_001"][()]
+                    Az = f.variables["az_001"][()]
+                    assert Ar.shape == (nphi, nz, nr)
+                    assert Aphi.shape == (nphi, nz, nr)
+                    assert Az.shape == (nphi, nz, nr)
                 r = np.linspace(rmin, rmax, nr)
                 phi = np.linspace(0, 2 * np.pi / nfp, nphi, endpoint=False)
                 z = np.linspace(zmin, zmax, nz)
@@ -1084,9 +1085,15 @@ def test_to_mgrid(self, include_potential=True):
                             np.testing.assert_allclose(Br[jphi, jz, jr], bs.B_cyl()[0, 0])
                             np.testing.assert_allclose(Bphi[jphi, jz, jr], bs.B_cyl()[0, 1])
                             np.testing.assert_allclose(Bz[jphi, jz, jr], bs.B_cyl()[0, 2])
-                            np.testing.assert_allclose(Ar[jphi, jz, jr], bs.A_cyl()[0, 0])
-                            np.testing.assert_allclose(Aphi[jphi, jz, jr], bs.A_cyl()[0, 1])
-                            np.testing.assert_allclose(Az[jphi, jz, jr], bs.A_cyl()[0, 2])
+                            if include_potential:
+                                np.testing.assert_allclose(Ar[jphi, jz, jr], bs.A_cyl()[0, 0])
+                                np.testing.assert_allclose(Aphi[jphi, jz, jr], bs.A_cyl()[0, 1])
+                                np.testing.assert_allclose(Az[jphi, jz, jr], bs.A_cyl()[0, 2])
+
+    def test_to_mgrid(self):
+        for include_potential in [True, False]:
+            with self.subTest(include_potential=include_potential):
+                self.subtest_to_mgrid(include_potential)
 
     def test_poloidal_field(self):
         B0 = 1.1
diff --git a/tests/field/test_mgrid.py b/tests/field/test_mgrid.py
index 0799f051f..c4872a247 100644
--- a/tests/field/test_mgrid.py
+++ b/tests/field/test_mgrid.py
@@ -56,7 +56,7 @@ def test_add_field_cylinder(self):
         assert mgrid.coil_names[0] == '__________test_coil___________'
         assert np.allclose(mgrid.br_arr[0], br)
 
-    def test_write(self): 
+    def test_write(self):
         mgrid = MGrid.from_file(test_file)
         with tempfile.TemporaryDirectory() as tmpdir:
             filename = Path(tmpdir) / 'mgrid.test.nc'