diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..7b54a3d
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,27 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
\ No newline at end of file
diff --git a/.gitignore~ b/.gitignore~
new file mode 100644
index 0000000..b123bf7
--- /dev/null
+++ b/.gitignore~
@@ -0,0 +1 @@
+*.egg
\ No newline at end of file
diff --git a/INGRID.egg-info/SOURCES.txt b/INGRID.egg-info/SOURCES.txt
index 7cbb69f..255dcb5 100644
--- a/INGRID.egg-info/SOURCES.txt
+++ b/INGRID.egg-info/SOURCES.txt
@@ -1,5 +1,6 @@
 README.md
 setup.py
+INGRID/DistribFunc.py
 INGRID/ImportEqFiles.py
 INGRID/__init__.py
 INGRID/geometry.py
diff --git a/INGRID/DistribFunc.py b/INGRID/DistribFunc.py
new file mode 100755
index 0000000..277e644
--- /dev/null
+++ b/INGRID/DistribFunc.py
@@ -0,0 +1,147 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Sep 17 10:04:20 2020
+
+@author: jguterl
+"""
+import types
+import numpy as np
+def DistribFunc(F:dict or str or types.FunctionType='x,x',**kwargs)->types.FunctionType:
+        
+        def get_lambdafunc( func: str) -> types.FunctionType:
+            """
+            Create a function from a string input.
+    
+            Will be used to generate a poloidal or radial transformation
+            function.
+    
+            Parameters
+            ----------
+            func : str
+                An expression to generate a function from.
+    
+            Returns
+            -------
+                A function generated from the str input.
+    
+            Examples
+            --------
+            When calling method ``get_func`` the user must provide a **string**
+            with the following general format:
+    
+            .. math::
+    
+                x, f(x)
+    
+            That is, the dependent variable and expression to evaluate are
+            separated with a comma character.
+    
+            The `Sympy` library supports most expressions that can be generated with
+            Python natively. See the `Sympy` documentation for advanced features.
+    
+            Defining a function representing f(x) = x ^ 2:
+    
+            >>> func = 'x, x ** 2'
+            >>> f = MyTopologyUtils.get_func(func)
+            >>> f(np.array([0, 1, 2, 3]))
+            array([0, 1, 4, 9])
+    
+            Defining a function representing f(x) = exp(x)
+    
+            >>> func = 'x, exp(x)'
+            >>> f = MyTopologyUtils.get_func(func)
+            >>> f(np.array([0, 1, 2, 3]))
+            array([ 1.        ,  2.71828183,  7.3890561 , 20.08553692])
+    
+            """
+    
+            def make_sympy_func(var, expression):
+                import sympy as sp
+                _f = sp.lambdify(var, expression, 'numpy')
+                return _f
+    
+            f_str_raw = func
+    
+            f_str_raw = f_str_raw.replace(' ', '')
+            delim = f_str_raw.index(',')
+    
+            var = f_str_raw[0: delim]
+            expression = f_str_raw[delim + 1:]
+    
+            func = make_sympy_func(var, expression)
+            # TODO: Check Normalization of the function to 1
+            return func
+    
+        def CheckLambdaFunction(expression: str, Verbose: bool = False) -> bool:
+            """
+            Check if a str is in the correct format for method ``get_func``
+    
+            Parameters
+            ----------
+            expression : str
+                Expression to check.
+    
+            Verbose : bool, optional
+                Print full output to terminal. Default to False
+    
+            Returns
+            -------
+                True if expression is valid. False otherwise.
+            """
+            ExpValid = False
+            try:
+                com = 'lambda {}:{}'.format(expression.split(',')[0], expression.split(',')[1])
+                if Verbose:
+                    print(com)
+                eval(com)
+                ExpValid = True
+            except:
+                ExpValid=False
+            return ExpValid
+        
+        if type(F)==dict:
+            Func=DistribFunc(**F)
+        elif type(F)==str:
+            if F=='exp':
+                if kwargs.get('alpha') is None:
+                    alpha=1
+                else:
+                    alpha=kwargs.get('alpha')
+                def Func(x):
+                    return (1-np.exp(-(x)/alpha))/(1-np.exp(-1/alpha))
+            elif F=='mexp':
+                if kwargs.get('alpha') is None:
+                    alpha=1
+                else:
+                    alpha=kwargs.get('alpha')
+                def Func(x):
+                    return 1-((1-np.exp(-((1-x)/alpha)))/(1-np.exp(-1/alpha)))    
+        
+            elif F=='pow':
+                if kwargs.get('alpha') is None:
+                    alpha=1
+                else:
+                    alpha=kwargs.get('alpha')
+                def Func(x):
+                    return x**alpha
+            elif F=='mpow':
+                if kwargs.get('alpha') is None:
+                    alpha=1
+                else:
+                    alpha=kwargs.get('alpha')
+                def Func(x):
+                    return 1-x**alpha
+            else:
+                if CheckLambdaFunction(F):    
+                    Func=get_lambdafunc(F)
+                else:
+                    raise ValueError('Unable to parse expression entry "{}".'.format(F))
+                    Func=None
+        else:
+            Func=F
+        
+        if not callable(Func):
+            Func=None
+            
+        return Func    
\ No newline at end of file
diff --git a/INGRID/geometry.py b/INGRID/geometry.py
index 2fc6e7d..03d112e 100644
--- a/INGRID/geometry.py
+++ b/INGRID/geometry.py
@@ -504,6 +504,19 @@ def plot_border(self, color: str = 'red', ax: 'matplotlib.axes.Axes' = None) ->
                 [self.vertices['SW'].y, self.vertices['NW'].y],
                 linewidth=1, color=color, label='cell')
 
+    def CollectBorder(self):
+        return np.array([(self.vertices[P].x, self.vertices[P].y) for P in ['NE','NW','SW','SE','NE']])
+
+    def CollectEdges(self):
+        Line1=[(self.vertices['NW'].x, self.vertices['NW'].y),(self.vertices['NE'].x, self.vertices['NE'].y)]
+        Line2=[(self.vertices['NE'].x, self.vertices['NE'].y),(self.vertices['SE'].x, self.vertices['SE'].y)]
+        Line3=[(self.vertices['SE'].x, self.vertices['SE'].y),(self.vertices['SW'].x, self.vertices['SW'].y)]
+        Line4=[(self.vertices['SW'].x, self.vertices['SW'].y),(self.vertices['NW'].x, self.vertices['NW'].y)]
+        return [Line1,Line2,Line3,Line4]
+    
+    
+
+        
     def plot_center(self, color='black', ax: 'matplotlib.axes.Axes' = None) -> None:
         """
         Plot the center of a Cell.
@@ -565,8 +578,16 @@ def __init__(self, lines: 'array-like', patch_name: str = '', PatchTagMap: dict
         self.E = lines[1]
         self.S = lines[2]
         self.W = lines[3]
+        self.E_vertices=[]
+        self.W_vertices=[]
+        self.N_vertices=[]
+        self.S_vertices=[]
         self.BoundaryPoints = {}
-
+        self.N_spl=None
+        self.E_spl=None
+        self.W_spl=None
+        self.N_spl=None
+        self.radial_spl=[]
         # This is the border for the fill function
         # It need to only include N and S lines
         self.RemoveDuplicatePoints()
@@ -576,7 +597,10 @@ def __init__(self, lines: 'array-like', patch_name: str = '', PatchTagMap: dict
         self.plate_location = plate_location
         self.patch_name = patch_name
         self.color = color
+        self.alpha=1
         self.Verbose = True
+        self.RadOrientation=1
+        self.PolOrientation=1
 
         # TODO: Populate these attributes when generating subgrid.
         #       Currently this is being done in the concat_grid method.
@@ -595,10 +619,17 @@ def __init__(self, lines: 'array-like', patch_name: str = '', PatchTagMap: dict
             'S': 'SOL',
             'C': 'Core'
         }
-
-        self.ax = None
         self.cell_grid = None
-
+    def ResetSplines(self):
+        self.radial_spl=[]
+        
+    def GetVertices(self):
+        return [getattr(self,'{}_vertices'.format(O)) for O in ['N','S','E','W']]
+    
+    def SetVertices(self,VerticesList):
+        for V,O in zip(VerticesList,['N','S','E','W']):
+            setattr(self,'{}_vertices'.format(O),V)
+            
     def plot_border(self, color='red', ax=None):
         """
         Draw solid borders around the patch.
@@ -625,7 +656,8 @@ def fill(self, color='lightsalmon', ax=None, alpha=1.0):
         x = np.array([p.x for p in self.p])
         y = np.array([p.y for p in self.p])
         arr = np.column_stack((x, y))
-        patch = Polygon(arr, fill=True, closed=True, color=color, label=self.get_tag(), alpha=alpha)
+        PatchLabel = self.patch_name + ' (' + UnfoldLabel(self.PatchLabelDoc, self.patch_name) + ')'
+        patch = Polygon(arr, fill=True, closed=True, color=color, label=PatchLabel, alpha=alpha)
         _ax = plt.gca() if ax is None else ax
         _ax.add_patch(patch)
         _ax.plot()
@@ -636,8 +668,71 @@ def plot_subgrid(self, ax=None, color='blue'):
         for row in self.cell_grid:
             for cell in row:
                 cell.plot_border(color=color, ax=ax)
-        plt.pause(0.1)
-
+        #plt.pause(0.1)
+        
+    def PlotSubGrid(self,ax=None,color=None):
+        if ax is None:
+            ax = plt.gca()
+        ListBorder=[]
+        ListEdges=[]
+        DataEdges=[]
+        
+        for row in self.cell_grid:
+            for cell in row:
+                Border=cell.CollectBorder()
+                
+                p=matplotlib.patches.Polygon(Border,edgecolor=color,fill=False)
+                ListBorder.append(p)
+                Edges=cell.CollectEdges()
+                ListEdges.extend(Edges)
+                for E in Edges:
+                    DataEdges.append(GetLength(E))
+                  
+
+        #Collection border of cell
+        C=matplotlib.collections.PatchCollection(ListBorder)
+        C.set_picker(True)
+        C.set_facecolor('')
+        if color is None:
+            C.set_edgecolor(self.color)
+        else:
+            C.set_edgecolor(color)
+        #ax.add_collection(C) 
+        
+        #Collection edges of cells
+        CE=matplotlib.collections.LineCollection(ListEdges,array=np.array(range(0,len(ListEdges))))
+        CE.set_picker(1)
+        ax.add_collection(CE) 
+        if color is None:
+            CE.set_color(self.color)
+            CE.set_edgecolors(self.color)
+            CE.set_facecolors(self.color)
+        else:
+            CE.set_color(color)
+            CE.set_edgecolors(color)
+            CE.set_facecolors(color)
+            
+        CE.set_alpha(self.alpha)
+        
+            
+        def onpick(evt):           
+            if evt.artist==CE:
+                annot.set_visible(False)
+                annot.xy = (evt.mouseevent.xdata, evt.mouseevent.ydata)
+                annot.set_text('length={:.2e} m'.format(DataEdges[evt.artist.get_array()[evt.ind[0]]]))
+                annot.set_visible(True)
+            if evt.mouseevent.button == 3:
+                annot.set_visible(False)  
+            ax.figure.canvas.draw_idle()                 
+        annot = ax.annotate("", xy=(0,0), xytext=(-20,20),textcoords="offset points",
+                                bbox=dict(boxstyle="round", fc="w"),
+                                arrowprops=dict(arrowstyle="->"))
+        annot.set_visible(False)
+        ax.figure.canvas.mpl_connect('pick_event', onpick)
+        
+                             
+        ax.autoscale_view()
+                
     def RemoveDuplicatePoints(self):
         for line in [self.N, self.E, self.S, self.W]:
             line.RemoveDuplicatePoints()
@@ -711,95 +806,29 @@ def as_np(self):
             patch_data.append(Z)
         patch_data = np.array(patch_data)
         cell_data = self.cell_grid_as_np()
+        cell_grid = self.cell_grid
         patch_settings = self.get_settings()
-        return np.array([patch_data, cell_data, patch_settings])
-
-    def make_subgrid(self, grid, np_cells=2, nr_cells=2, _poloidal_f=lambda x: x, _radial_f=lambda x: x, verbose=False, visual=False, ShowVertices=False):
-        """
-        Generate a refined grid within a patch.
-        This 'refined-grid' within a Patch is a collection
-        of num x num Cell objects
-
-        Parameters
-        ----------
-        grid : Ingrid
-                To be used for obtaining Efit data and all
-                other class information.
-        num  : int, optional
-                Number to be used to generate num x num
-                cells within our Patch.
-        """
-
-        def psi_parameterize(grid, r, z):
-            """
-            Helper function to be used to generate a
-            list of values to parameterize a spline
-            in Psi. Returns a list to be used for splprep only
-            """
-            vmax = grid.PsiNorm.get_psi(r[-1], z[-1])
-            vmin = grid.PsiNorm.get_psi(r[0], z[0])
-
-            vparameterization = np.empty(0)
-            for i in range(len(r)):
-                vcurr = grid.PsiNorm.get_psi(r[i], z[i])
-                vparameterization = np.append(vparameterization, abs((vcurr - vmin) / (vmax - vmin)))
-
-            return vparameterization
-
-        def psi_test(leg, grid):
-            """
-            Determine whether a leg's orientation is increasing in psi
-            by checking the endpoints of a given leg.
-
-            Returns True if final Point entry is greater in psi than the
-            first Point entry. Returns False otherwise.
-            """
-            p1 = leg.p[0]
-            p2 = leg.p[-1]
-
-            return True if grid.PsiNorm.get_psi(p2.x, p2.y) > grid.PsiNorm.get_psi(p1.x, p1.y) else False
-
-        def set_dynamic_step(ratio=0.01):
-            """
-            Compute a dt value for line_tracing that occurs during subgrid generation of
-            this Patch object. The value is taken to be a ratio of the average length of
-            a Patch in the radial direction.
-            """
-            d = 0
-            for i in range(nr_lines):
-                d += np.sqrt((self.E_vertices[i].x - self.W_vertices[i].x)**2 + (self.E_vertices[i].y - self.W_vertices[i].y)**2)
-            dynamic_step = d / nr_lines
-            # print('Dynamic-Step value came out to be: {}\n'.format(dynamic_step * ratio))
-            return dynamic_step * ratio
-
-        # Allocate space for collection of cell objects.
-        # Arbitrary 2D container for now.
-        cell_grid = []
-
-        if verbose:
-            print('Constructing grid for patch "{}" with dimensions (np, nr) = ({}, {})'.format(self.patch_name, np_cells, nr_cells))
-            print(np_cells)
-            print(nr_cells)
-        np_lines = np_cells + 1
-        nr_lines = nr_cells + 1
-        if verbose: print(' # Create B-Splines along the North and South boundaries.')
+        return np.array([patch_data, cell_data, patch_settings,self.GetVertices(),cell_grid])
+    
+    
+    def CreateBoundarySplines(self,grid):
         # Create B-Splines along the North and South boundaries.
+        if self.Verbose: print(' # Create B-Splines along the North and South boundaries.')
         N_vals = self.N.fluff()
-
         self.N_spl, uN = splprep([N_vals[0], N_vals[1]], s=0)
         # Reverse the orientation of the South line to line up with the North.
 
         S_vals = self.S.reverse_copy().fluff()
         self.S_spl, uS = splprep([S_vals[0], S_vals[1]], s=0)
-        if verbose: print(' # Create B-Splines along West boundaries.')
+        if self.Verbose: print(' # Create B-Splines along West boundaries.')
         # Create B-Splines along the East and West boundaries.
         # Parameterize EW splines in Psi
         try:
             #Cannot fluff with too many points
             n = 500 if len(self.W.p) < 50 else 100
            # W_vals = self.W.reverse_copy().fluff(num = n)
-            W_vals = self.W.reverse_copy().fluff(n, verbose=verbose)
-            Psi = psi_parameterize(grid, W_vals[0], W_vals[1])
+            W_vals = self.W.reverse_copy().fluff(n, verbose=self.Verbose)
+            Psi = self.psi_parameterize(grid, W_vals[0], W_vals[1])
             self.W_spl, uW = splprep([W_vals[0], W_vals[1]], u=Psi, s=10)
         except Exception as e:
             exc_type, exc_obj, tb = sys.exc_info()
@@ -810,16 +839,17 @@ def set_dynamic_step(ratio=0.01):
             line = linecache.getline(filename, lineno, f.f_globals)
             print('EXCEPTION IN ({}, LINE {} "{}"): {}'.format(filename, lineno, line.strip(), exc_obj))
 
-        if verbose: print(' # Create B-Splines along the East boundaries.')
+        if self.Verbose: print(' # Create B-Splines along the East boundaries.')
         try:
             n = 500 if len(self.E.p) < 50 else 100
             E_vals = self.E.fluff(num=n)
-            self.E_spl, uE = splprep([E_vals[0], E_vals[1]], u=psi_parameterize(grid, E_vals[0], E_vals[1]), s=10)
+            self.E_spl, uE = splprep([E_vals[0], E_vals[1]], u=self.psi_parameterize(grid, E_vals[0], E_vals[1]), s=10)
         except Exception as e:
             print(' Number of points on the boundary:', len(self.E.p))
             plt.plot(E_vals[0], E_vals[1], '.', color='black')
             print(repr(e))
-        if verbose: print(' #check plate_patch')
+        if self.Verbose: print(' #check plate_patch')
+        # ACCURACY ON PSI_MIN SIDE OF W IDL LINE ISSUE.
 
         if self.plate_patch:
 
@@ -845,7 +875,7 @@ def find_nearest(array, value):
                 plate_north = [N_vals[0][-1], N_vals[1][-1]]
                 plate_south = [S_vals[0][-1], S_vals[1][-1]]
 
-            if verbose:
+            if self.Verbose:
                 print('=' * 80 + '\n')
                 print('Parameterization of Target Plate in PSI:')
                 print(_u)
@@ -860,19 +890,19 @@ def find_nearest(array, value):
                 try:
                     plate_north_index = lookup[find_nearest(_u, fsolve(f, 0, args=(U_spl, plate_north[0], plate_north[1])))]
                 except:
-                    if verbose: print('NorthIndex: ERROR IN PARAMETERIZATION IN PSI')
+                    print('NorthIndex: ERROR IN PARAMETERIZATION IN PSI')
             try:
                 plate_south_index = lookup[find_nearest(_u, brentq(f, _u[0], _u[-1], args=(U_spl, plate_south[0], plate_south[1])))]
             except ValueError:
                 try:
                     plate_south_index = lookup[find_nearest(_u, fsolve(f, 0, args=(U_spl, plate_south[0], plate_south[1])))]
                 except:
-                    if verbose: print('SouthIndex: ERROR IN PARAMETERIZATION IN PSI')
+                    print('SouthIndex: ERROR IN PARAMETERIZATION IN PSI')
             if plate_south_index > plate_north_index:
                 plate_north_index, plate_south_index = plate_south_index, plate_north_index
 
             U_vals = [U_vals[0][plate_south_index:plate_north_index + 1], U_vals[1][plate_south_index:plate_north_index + 1]]
-            U_spl, _u = splprep([U_vals[0], U_vals[1]], u=psi_parameterize(grid, U_vals[0], U_vals[1]), s=0)
+            U_spl, _u = splprep([U_vals[0], U_vals[1]], u=self.psi_parameterize(grid, U_vals[0], U_vals[1]), s=0)
 
             if self.plate_location == 'W':
                 W_vals = U_vals
@@ -882,9 +912,86 @@ def find_nearest(array, value):
                 E_vals = U_vals
                 self.E_spl = U_spl
                 uE = _u
+                
+                
+    @staticmethod
+    def psi_parameterize(grid, r, z):
+            """
+            Helper function to be used to generate a
+            list of values to parameterize a spline
+            in Psi. Returns a list to be used for splprep only
+            """
+            vmax = grid.PsiNorm.get_psi(r[-1], z[-1])
+            vmin = grid.PsiNorm.get_psi(r[0], z[0])
+
+            vparameterization = np.empty(0)
+            for i in range(len(r)):
+                vcurr = grid.PsiNorm.get_psi(r[i], z[i])
+                vparameterization = np.append(vparameterization, abs((vcurr - vmin) / (vmax - vmin)))
+
+            return vparameterization
+        
+    def make_subgrid(self, grid, np_cells=2, nr_cells=2, _poloidal_f=lambda x: x, _radial_f=lambda x: x, verbose=False, visual=False, ShowVertices=False):
+        """
+        Generate a refined grid within a patch.
+        This 'refined-grid' within a Patch is a collection
+        of num x num Cell objects
+
+        Parameters
+        ----------
+        grid : Ingrid
+                To be used for obtaining Efit data and all
+                other class information.
+        num  : int, optional
+                Number to be used to generate num x num
+                cells within our Patch.
+        """
+
+        
+
+        def psi_test(leg, grid):
+            """
+            Determine whether a leg's orientation is increasing in psi
+            by checking the endpoints of a given leg.
+
+            Returns True if final Point entry is greater in psi than the
+            first Point entry. Returns False otherwise.
+            """
+            p1 = leg.p[0]
+            p2 = leg.p[-1]
+
+            return True if grid.PsiNorm.get_psi(p2.x, p2.y) > grid.PsiNorm.get_psi(p1.x, p1.y) else False
+
+        def set_dynamic_step(ratio=0.01):
+            """
+            Compute a dt value for line_tracing that occurs during subgrid generation of
+            this Patch object. The value is taken to be a ratio of the average length of
+            a Patch in the radial direction.
+            """
+            d = 0
+            for i in range(nr_lines):
+                d += np.sqrt((self.E_vertices[i].x - self.W_vertices[i].x)**2 + (self.E_vertices[i].y - self.W_vertices[i].y)**2)
+            dynamic_step = d / nr_lines
+            # print('Dynamic-Step value came out to be: {}\n'.format(dynamic_step * ratio))
+            return dynamic_step * ratio
+
+        # Allocate space for collection of cell objects.
+        # Arbitrary 2D container for now.
+        cell_grid = []
+
+        if verbose:
+            print('Constructing grid for patch "{}" with dimensions (np, nr) = ({}, {})'.format(self.patch_name, np_cells, nr_cells))
+            print(np_cells)
+            print(nr_cells)
+        np_lines = np_cells + 1
+        nr_lines = nr_cells + 1
+        
+        #self.CreatBoundarySpleen()
 
         # Generate our sub-grid anchor points along the North
         # and South boundaries of our patch.
+        self.W_vertices_old=self.W_vertices
+        
         self.N_vertices = []
         self.S_vertices = []
         self.E_vertices = []
@@ -898,7 +1005,7 @@ def find_nearest(array, value):
                 _n = splev(_poloidal_f(i / (np_lines - 1)), self.N_spl)
                 self.N_vertices.append(Point((float(_n[0]), float(_n[1]))))
         else:
-            if verbose: print('Find boundary points at face "N" for {}:{}'.format(self.patch_name, self.BoundaryPoints.get('N')))
+            if self.Verbose: print('Find boundary points at face "N" for {}:{}'.format(self.patch_name, self.BoundaryPoints.get('N')))
             self.N_vertices = self.BoundaryPoints.get('N')
 
         if self.BoundaryPoints.get('S') is None:
@@ -911,15 +1018,15 @@ def find_nearest(array, value):
         u = [_radial_f(i / (nr_lines - 1)) for i in range(nr_lines)]
 
         Npts = 1000
-        xy = splev(np.linspace(0, 1, Npts), self.E_spl)
-
+        #xy = splev(np.linspace(0, 1, Npts), self.E_spl)
+        
         if self.BoundaryPoints.get('W') is None:
             for i in range(nr_lines):
                 _w = splev(u[i], self.W_spl)
                 self.W_vertices.append(Point((float(_w[0]), float(_w[1]))))
         else:
             self.W_vertices = self.BoundaryPoints.get('W')
-
+        
         if self.BoundaryPoints.get('E') is None:
             for i in range(nr_lines):
                 _e = splev(u[i], self.E_spl)
@@ -943,11 +1050,11 @@ def find_nearest(array, value):
 
             else:
                 color = self.color
-                ms = 6
-                markers = ['o', 'X', 's', 'D']
+                ms = 10
+                markers = ['o', 'P', 's', 'D']
             for vertices, mark in zip([self.W_vertices, self.E_vertices, self.N_vertices, self.S_vertices], markers):
                 for p in vertices:
-                    plt.plot(p.x, p.y, '.', color=color, markersize=ms, marker=mark, markeredgecolor='black')
+                    plt.plot(p.x, p.y, '.', color=color, markersize=ms, marker=mark, markeredgecolor=color)
         # Radial lines of Psi surfaces. Ordered with increasing magnitude, starting with
         # the South boundary of the current Patch, and ending with the North boundary of
         # this current Patch. These will serve as delimiters when constructing cells.
@@ -956,32 +1063,38 @@ def find_nearest(array, value):
         dynamic_step = set_dynamic_step()
         if verbose: print('# Interpolate radial lines between North and South patch boundaries..')
         # Interpolate radial lines between North and South patch boundaries.
+        if hasattr(self,'radial_spl'):
+            self.radial_spl_old=self.radial_spl
+        else:
+            self.radial_spl_old=[]
         self.radial_spl = []
         temp_vertices = []
         temp_vertices.append(self.N.p[-1])
         for i in range(len(self.W_vertices) - 2):
             #TODO: parallelize tracing of radial lines (draw_line function must be "externalized" in the scope of the script)
+            # if (len(self.W_vertices)!=len(self.W_vertices_old) or self.W_vertices[i+1]!=self.W_vertices_old[i+1]) or len(self.radial_spl_old)!=len(self.W_vertices)-2:
+    
             radial_lines.append(grid.LineTracer.draw_line(self.W_vertices[i + 1], {'line': self.E}, option='theta',
-                direction='cw', show_plot=visual, dynamic_step=dynamic_step, text=verbose))
+                direction='cw', show_plot=0, dynamic_step=dynamic_step))
             temp_vertices.append(radial_lines[-1].p[-1])
             radial_vals = radial_lines[i + 1].fluff(1000)
             Radial_spl, uR = splprep([radial_vals[0], radial_vals[1]], s=0)
             self.radial_spl.append(Radial_spl)
+            # else:
+            #     self.radial_spl.append(self.radial_spl_old[i])
+                
             vertex_list = []
             for j in range(np_lines):
                 u = _poloidal_f(j / (np_lines - 1))
                 _r = splev(u, self.radial_spl[i])
                 Pt = Point((float(_r[0]), float(_r[1])))
                 if self.distortion_correction['active'] and j > 0 and j < np_lines - 1:
-                    Res = self.distortion_correction['resolution']
-                    ThetaMin = self.distortion_correction['theta_min']
-                    ThetaMax = self.distortion_correction['theta_max']
                     umin = _poloidal_f((j - 1) / (np_lines - 1))
                     umax = _poloidal_f((j + 1) / (np_lines - 1))
                     Pt1 = radial_vertices[i][j]
                     Pt2 = radial_vertices[i][j - 1]
                     Tag = '---- Correcting points: {},{}'.format(i, j)
-                    Pt = CorrectDistortion(u, Pt, Pt1, Pt2, self.radial_spl[i], ThetaMin, ThetaMax, umin, umax, Res, visual, Tag, verbose)
+                    Pt = CorrectDistortion(u, Pt, Pt1, Pt2, self.radial_spl[i],umin,umax,Tag,**self.distortion_correction)
 
                 vertex_list.append(Pt)
                 if visual:
@@ -998,15 +1111,12 @@ def find_nearest(array, value):
             for j in range(1, np_lines - 1):
                 u = _poloidal_f(j / (np_lines - 1))
                 Pt = self.S_vertices[j]
-                Res = self.distortion_correction['resolution']
-                ThetaMin = self.distortion_correction['theta_min']
-                ThetaMax = self.distortion_correction['theta_max']
                 umin = _poloidal_f((j - 1) / (np_lines - 1))
                 umax = _poloidal_f((j + 1) / (np_lines - 1))
                 Pt1 = radial_vertices[-1][j]
                 Pt2 = radial_vertices[-1][j - 1]
                 Tag = '---- Correcting south boundary points:{}'.format(j)
-                self.S_vertices[j] = CorrectDistortion(u, Pt, Pt1, Pt2, self.S_spl, ThetaMin, ThetaMax, umin, umax, Res, visual, Tag, verbose)
+                self.S_vertices[j] = CorrectDistortion(u, Pt, Pt1, Pt2, self.S_spl, umin, umax,  Tag,**self.distortion_correction)
         radial_vertices.append(self.S_vertices)
         if verbose: print('# Create Cells: South boundary -> North Boundary')
         # Create Cells: South boundary -> North Boundary
@@ -1298,34 +1408,36 @@ def trim_geometry(geoline, start, end):
     return trim
 
 
-def CorrectDistortion(u, Pt, Pt1, Pt2, spl, ThetaMin, ThetaMax, umin, umax, Resolution, visual, Tag, MinTol=1.02, MaxTol=0.98, Verbose=False):
-    MaxIter = Resolution * 10
-    dumax = (umax - u) / Resolution
-    dumin = (u - umin) / Resolution
+def CorrectDistortion(u, Pt, Pt1, Pt2, spl, umin, umax , Tag, theta_min=60, theta_max=120,resolution=1000,min_tol=1.05, max_tol=0.95, zero_span=0.90,max_span=1.0,verbose=False,visual=False,**kwargs):
+    MaxIter = resolution * 10
+    dumax = (umax - u) / resolution
+    dumin = (u - umin) / resolution
     Theta = Line([Pt1, Pt]).GetAngle(Line([Pt1, Pt2]))
-    if Verbose: print(f'Current theta value: {Theta}')
-    if Theta < ThetaMin or Theta > ThetaMax:
-        if Verbose: print('{}: u={};Theta={};ThetaMin={};ThetaMax={}'.format(Tag, u, Theta, ThetaMin, ThetaMax))
+    if verbose: print(f'Current theta value: {Theta}')
+    if Theta < theta_min or Theta > theta_max:
+        if verbose: print('{}: u={};Theta={};ThetaMin={};ThetaMax={}'.format(Tag, u, Theta, theta_min, theta_max))
         if visual:
                 plt.plot(Pt.x, Pt.y, '.', color='red', markersize=8, marker='o')
                 plt.show()
                 plt.draw()
         icount = 0
         color = 'purple'
-        while Theta < ThetaMin or Theta > ThetaMax:
+        u0=u
+
+        while Theta < theta_min or Theta > theta_max:
             icount += 1
-            if Theta < ThetaMin:
+            if Theta < theta_min:
                u = u + dumax
-            elif Theta > ThetaMax:
+            elif Theta > theta_max:
                u = u - dumin
-            if (u > umax * MaxTol or u < umin * MinTol) or icount >= MaxIter:
-               if Verbose: print('[{}]>>>> umax={} umin={};u:{};Theta={}'.format(icount, umax, umin, u, Theta))
+            if (u-u0) > (umax* max_tol-u0)*max_span or (u-u0) < (umin * min_tol-u0)*max_span or icount >= MaxIter or (abs(umin)<1e-12 and u<u0*zero_span):
+               if verbose: print('[{}]>>>> umax={} umin={};u:{};Theta={}'.format(icount, umax, umin, u, Theta))
                color = 'gray'
                break
             _r = splev(u, spl)
             Pt = Point((_r[0], _r[1]))
             Theta = Line([Pt1, Pt]).GetAngle(Line([Pt1, Pt2]))
-        if Verbose: print('[{}]>>>> u:{};Theta={}'.format(icount, u, Theta))
+        if verbose: print('[{}]>>>> u:{};Theta={}'.format(icount, u, Theta))
         if visual:
             plt.plot(Pt.x, Pt.y, '.', color=color, markersize=8, marker='s')
     return Pt
@@ -1490,7 +1602,7 @@ def segment_intersect(line1, line2, verbose=False):
 
 def UnfoldLabel(Dic: dict, Name: str) -> str:
     """
-    Unfold Patch label (e.g. "C1" -> "Inner Core Top")
+    Unfold Patch label (e.g. "ICT" -> "Inner Core Top")
 
     Parameters
     ----------
@@ -1515,3 +1627,6 @@ def UnfoldLabel(Dic: dict, Name: str) -> str:
         return ''.join(Output)
     else:
         return ''
+    
+def GetLength(E):
+    return np.sqrt((E[1][0]-E[0][0])**2+(E[1][1]-E[1][1])**2)
diff --git a/INGRID/ingrid.py b/INGRID/ingrid.py
index 7aad095..c13b861 100644
--- a/INGRID/ingrid.py
+++ b/INGRID/ingrid.py
@@ -17,18 +17,17 @@
     pass
 import matplotlib.pyplot as plt
 import pathlib
-import inspect
-from scipy.optimize import root, minimize
+
 
 import yaml as yml
 import os
 from pathlib import Path
 from time import time
-from collections import OrderedDict
 
-from omfit_eqdsk import OMFITgeqdsk
+
+from INGRID.DistribFunc import *
 from INGRID.interpol import EfitData
-from INGRID.interpol import Bicubic
+
 from INGRID.utils import IngridUtils
 from INGRID.topologies.snl import SNL
 from INGRID.topologies.sf15 import SF15
@@ -38,8 +37,10 @@
 from INGRID.topologies.sf135 import SF135
 from INGRID.topologies.sf165 import SF165
 from INGRID.topologies.udn import UDN
-from INGRID.line_tracing import LineTracing
 from INGRID.geometry import Point, Line
+import types
+
+
 
 
 def QuickStart() -> None:
@@ -193,6 +194,45 @@ def SaveSettingsFile(self, fname: str = '', settings: dict = {}) -> 'Path':
             yml.dump(self.settings, f)
 
         return fpath
+    
+    
+    def StoreSettings(self, fname: str = '', settings: dict = {}) -> 'Path':
+        """
+        Save a new settings .yml file.
+
+        Parameters
+        ----------
+        fname : optional
+            Name of new settings '.yml' file.
+            If default value of '', then Ingrid will generate
+            a '.yml' file named 'INGRID_Session' appended with
+            a timestamp.
+
+        settings : optional
+            Ingrid settings dictionary to dump into
+            the '.yml' file. Defaults to empty dict which produces a
+            template settings file.
+
+        Returns
+        -------
+            A Path instance representing the saved YAML file.
+        """
+
+        if fname == '':
+            fname = 'INGRID_Session' + str(int(time())) + '.yml'
+
+        fpath = Path(fname)
+
+        if fpath.suffix == '':
+            fname += '.yml'
+        elif fpath.suffix != '.yml':
+            fname = fname[:-4] + '.yml'
+
+        with open(fname, 'w') as f:
+            yml.dump(self.settings, f)
+
+        return fpath
+
 
     def PrintSummaryInput(self) -> None:
         """
@@ -778,9 +818,8 @@ def PlotTargetPlate(self, plate_key: str, color: str = 'red', ax: object = None)
             ax = plt.gca()
 
         if plate_key in [k for k in self.PlateData.keys()]:
-            self.RemovePlotLine(label=plate_key, ax=ax)
             if type(self.PlateData[plate_key]) is Line:
-                print(f"# Plotting plate '{plate_key}'")
+                self.RemovePlotLine(label=plate_key, ax=ax)
                 self.PlateData[plate_key].plot(label=plate_key, color=color)
             else:
                 pass
@@ -1050,7 +1089,7 @@ def PlotPatches(self) -> None:
         if self.settings['grid_settings']['patch_generation']['strike_pt_loc'] == 'target_plates':
             self.RemovePlotLine(label='limiter', ax=self.PatchAx)
 
-    def PlotSubgrid(self) -> None:
+    def PlotSubgrid(self,NewFig=False,**kwargs) -> None:
         """
         Plot the grid that was generated with method 'CreateSubgrid'.
         """
@@ -1058,9 +1097,13 @@ def PlotSubgrid(self) -> None:
             plt.close(self._SubgridFig)
         except:
             pass
-        self._SubgridFig = plt.figure('INGRID: ' + self.CurrentTopology.config + ' Grid', figsize=(6, 10))
-        self._SubgridAx = self._SubgridFig.add_subplot(111)
-        self.CurrentTopology.grid_diagram(fig=self._SubgridFig, ax=self._SubgridAx)
+        if NewFig:
+            self._SubgridFig = plt.figure('INGRID: ' + self.CurrentTopology.config + ' Grid', figsize=(6, 10))
+            self._SubgridAx = self._SubgridFig.add_subplot(111)
+        else:
+            self._SubgridFig = plt.gcf()
+            self._SubgridAx = plt.gca()
+        self.CurrentTopology.grid_diagram(fig=self._SubgridFig, ax=self._SubgridAx,**kwargs)
 
     def AutoRefineMagAxis(self) -> None:
         """
@@ -1268,11 +1311,12 @@ def ConstructPatches(self) -> None:
         self.OrderTargetPlates()
         self.PrepLineTracing()
         self.CurrentTopology.construct_patches()
-        self.CurrentTopology.GroupPatches()
+        self.CurrentTopology.GroupPatches() 
         self.CheckPatches()
-
-        if self.settings['patch_data']['preferences']['new_file']:
-            self.SavePatches(self.settings['patch_data']['preferences']['new_fname'])
+        self.CurrentTopology.SetBoundarySplines()
+        #not practical in a workflow
+        #if self.settings['patch_data']['preferences']['new_file']:
+        #   self.SavePatches(self.settings['patch_data']['preferences']['new_fname'])
 
     def SavePatches(self, fname: str = '') -> None:
         """
@@ -1299,6 +1343,35 @@ def SavePatches(self, fname: str = '') -> None:
             tail = '.npy'
         print(f"# Saved patch data for file {Path(fname).name + tail}")
 
+        
+    def SuperSavePatches(self, fname: str = '') -> None:
+        """
+        Save patches as '.npy' file for later reconstruction in Ingrid.
+
+        This file contains the information required to reconstruct patches
+        at a later time and bypass the line_tracing.
+
+        Parameters
+        ----------
+        fname : str, optional
+            Name of file/location for patch data.
+        """
+        if fname in ['', None]:
+            fname = self.CurrentTopology.config + '_patches_' + str(int(time()))
+        self.CurrentTopology.ResetHeavyAttr()
+        
+        Data={'patches':self.CurrentTopology.patches,
+              'config':self.CurrentTopology.config,
+              'xpt_data':self.CurrentTopology.LineTracer.NSEW_lookup
+              }
+        np.save(fname, Data)
+        
+        if Path(fname).suffix == '.npy':
+            tail = ''
+        else:
+            tail = '.npy'
+        print(f"# Saved patch data for file {Path(fname).name + tail}")
+
     def LoadPatches(self, fname: str = '') -> None:
         """
         Load patches stored in an Ingrid generated '.npy' file.
@@ -1318,7 +1391,7 @@ def LoadPatches(self, fname: str = '') -> None:
         if fname.strip() == '':  # Check if settings contains patch data.
             fname = self.settings['patch_data']['file']
         config, xpt_data, patches = self.ReconstructPatches(fname)
-        self.PopulateSettings(Ingrid.ReadYamlFile(self.InputFile))
+        #self.PopulateSettings(Ingrid.ReadYamlFile(self.InputFile))
         self.StartSetup()
         self.LineTracer.NSEW_lookup = xpt_data
         self.SetTopology(config)
@@ -1326,8 +1399,37 @@ def LoadPatches(self, fname: str = '') -> None:
         self.CurrentTopology.OrderPatches()
         self.CurrentTopology.SetupPatchMatrix()
         self.CheckPatches()
+        
+    def SuperLoadPatches(self, fname: str = '') -> None:
+        """
+        Load patches stored in an Ingrid generated '.npy' file.
+
+        Parameters
+        ----------
+        fname : str, optional
+            Path to patch data.
+                If no fname is provided to method 'LoadPatches', Ingrid code will check the settings
+                'dict' for a file under entry ``settings['patch_data']['file']``
+
+        """
+
+        if type(fname) is not str:
+            raise ValueError('# User did not provide a string to patch data.')
+
+        if fname.strip() == '':  # Check if settings contains patch data.
+            fname = self.settings['patch_data']['file']
+        Data=np.load(fname,allow_pickle=True).tolist()
+        #config, xpt_data, patches = self.ReconstructPatches(fname)
+        #self.PopulateSettings(Ingrid.ReadYamlFile(self.InputFile))
+        self.StartSetup()
+        self.LineTracer.NSEW_lookup = Data['xpt_data']
+        self.SetTopology(Data['config'])
+        self.CurrentTopology.patches = Data['patches']
+        self.CurrentTopology.OrderPatches()
+        self.CurrentTopology.SetupPatchMatrix()
+        self.CheckPatches()
 
-    def CreateSubgrid(self, NewFig: bool = True, ShowVertices: bool = False) -> None:
+    def CreateSubgrid(self, NewFig: bool = True, ShowVertices: bool = False,**kwargs) -> None:
         """
         Refine a generated patch map into a grid for exporting.
 
@@ -1350,9 +1452,9 @@ def CreateSubgrid(self, NewFig: bool = True, ShowVertices: bool = False) -> None
             the particular patch.
 
         """
-        self.CurrentTopology.construct_grid()
+        self.CurrentTopology.construct_grid( ShowVertices=ShowVertices,**kwargs)
 
-    def ExportGridue(self, fname: str = 'gridue') -> None:
+    def ExportGridue(self, fname: str = 'gridue',**kwargs) -> None:
         """
         Export a gridue file for the created grid.
 
@@ -1362,12 +1464,9 @@ def ExportGridue(self, fname: str = 'gridue') -> None:
             Name of gridue file to save.
 
         """
-        if type(self.CurrentTopology) in [SNL]:
-            if self.WriteGridueSNL(self.CurrentTopology.gridue_settings, fname):
-                print(f"# Successfully saved gridue file as '{fname}'")
-        elif type(self.CurrentTopology) in [SF15, SF45, SF75, SF105, SF135, SF165, UDN]:
-            if self.WriteGridueDNL(self.CurrentTopology.gridue_settings, fname):
-                print(f"# Successfully saved gridue file as '{fname}'")
+        self.CurrentTopology.PrepGridue(**kwargs)
+        self.CurrentTopology.WriteGridue(fname,**kwargs)
+        print(f"# Successfully saved gridue file as '{fname}'")
 
     @staticmethod
     def ImportGridue(self, fname: str = 'gridue') -> dict:
diff --git a/INGRID/line_tracing.py b/INGRID/line_tracing.py
index 79a9b6f..9846e9e 100644
--- a/INGRID/line_tracing.py
+++ b/INGRID/line_tracing.py
@@ -847,7 +847,6 @@ def fend_Z(y):
 
         if Verbose:
             print('Drew for {} seconds\n'.format(end - start))
-        print('')
         return Line(line)
 
     def PsiCostFunc(self, xy):
diff --git a/INGRID/topologies/snl.py b/INGRID/topologies/snl.py
index 3067e61..27cf268 100644
--- a/INGRID/topologies/snl.py
+++ b/INGRID/topologies/snl.py
@@ -340,7 +340,7 @@ def OrderPatches(self):
         patches = ['A2', 'B2', 'C2', 'D2', 'E2', 'F2', 'A1', 'F1', 'B1', 'C1', 'D1', 'E1']
         self.patches = OrderedDict([(pname, self.patches[pname]) for pname in patches])
 
-    def set_gridue(self):
+    def set_gridue(self,**kwargs):
         """
         Prepares a ``gridue_settings`` dictionary with required data
         for writing a gridue file.
@@ -394,8 +394,71 @@ def set_gridue(self):
                     b[i][j][k] = _b
 
         self.gridue_settings = {
-            'nxm': nxm, 'nym': nym, 'ixpt1': ixpt1, 'ixpt2': ixpt2, 'iyseptrx1': iyseparatrix1,
+            'nxm': nxm, 'nym': nym, 'ixpt1': ixpt1, 'ixpt2': ixpt2, 'iyseptrx1': iyseparatrix1,'simagxs':self.OMFIT_psi['SIMAG'] ,'sibdrys':self.OMFIT_psi['SIBRY'],
             'rm': self.rm, 'zm': self.zm, 'psi': psi, 'br': br, 'bz': bz, 'bpol': bpol, 'bphi': bphi, 'b': b
         }
+        
+    
+    
+    def WriteGridue(self, fname: str = 'gridue',NewFormat=False,**kwargs) -> bool:
+        """
+        Write a gridue file for a single-null configuration.
+
+        Parameters
+        ----------
+        gridue_settings : dict
+            A dictionary containing grid data to be written to the gridue file.
+
+        fname : str, optional
+            The file name/path to save the gridue file to.
+            Defaults to 'gridue'.
+
+        Returns
+        -------
+            True if file was written with no errors
+        """
+
+        def format_header(gridue):
+            iheader_items = ['nxm', 'nym', 'ixpt1', 'ixpt2', 'iyseptrx1']
+            rheader_items = ['simagxs','sibdrys']
+            
+            Str=''.join(['{}'.format(gridue[item]).rjust(4) for item in iheader_items])
+            if NewFormat:
+                Str=Str+''.join(['{:16.10f}'.format(gridue[item]).rjust(17) for item in rheader_items])
+            return Str+ '\n'
+
+        def format_body(data):
+
+            delim_val = 0
+            delim_char = ''
+            body = ''
+
+            for n in range(5):
+                for j in range(len(data[0])):
+                    for i in range(len(data)):
+                        delim_val += 1
+                        val = np.format_float_scientific(data[i][j][n], precision=15, unique=False).rjust(23).replace('e', 'D')
+                        if delim_val == 3:
+                            delim_val = 0
+                            delim_char = '\n'
+                        body += val + delim_char
+                        delim_char = ''
+
+            if delim_val % 3 != 0:
+                body += '\n'
+
+            return body
+
+        f = open(fname, mode='w')
+        f.write(format_header(self.gridue_settings) + '\n')
+
+        body_items = ['rm', 'zm', 'psi', 'br', 'bz', 'bpol', 'bphi', 'b']
+        for item in body_items:
+            f.write(format_body(self.gridue_settings[item]) + '\n')
+
+        runidg = 'iogridue'
+        f.write(runidg + '\n')
+
+        f.close()
 
-        return self.gridue_settings
+        return True
diff --git a/INGRID/topologies/udn.py b/INGRID/topologies/udn.py
index 6f07beb..cfc2027 100644
--- a/INGRID/topologies/udn.py
+++ b/INGRID/topologies/udn.py
@@ -618,10 +618,78 @@ def set_gridue(self) -> dict:
                     b[i][j][k] = _b
 
         self.gridue_settings = {
-            'nxm': nxm, 'nym': nym, 'iyseparatrix1': iyseparatrix1, 'iyseparatrix2': iyseparatrix2,
+            'nxm': nxm, 'nym': nym, 'iyseparatrix1': iyseparatrix1, 'iyseparatrix2': iyseparatrix2,'simagxs':self.OMFIT_psi['SIMAG'] ,'sibdrys':self.OMFIT_psi['SIBRY'],
             'ix_plate1': ix_plate1, 'ix_cut1': ix_cut1, 'ix_cut2': ix_cut2, 'ix_plate2': ix_plate2, 'iyseparatrix3': iyseparatrix3,
             'iyseparatrix4': iyseparatrix4, 'ix_plate3': ix_plate3, 'ix_cut3': ix_cut3, 'ix_cut4': ix_cut4, 'ix_plate4': ix_plate4,
             'rm': self.rm, 'zm': self.zm, 'psi': psi, 'br': br, 'bz': bz, 'bpol': bpol, 'bphi': bphi, 'b': b, '_FILLER_': -1
         }
+    
+    
+    def WriteGridue(self, gridue_settings: dict, fname: str = 'gridue',**kwargs) -> bool:
+        """
+        Write a gridue file for a double-null configuration.
+
+        Parameters
+        ----------
+        gridue_settings : dict
+            A dictionary containing grid data to be written to the gridue file.
+
+        fname : str, optional
+            The file name/path to save the gridue file to.
+            Defaults to 'gridue'.
+
+        Returns
+        -------
+            True if file was written with no errors
+        """
 
-        return self.gridue_settings
+        def format_header(gridue):
+            header_rows = [
+                ['nxm', 'nym'],
+                ['iyseparatrix1', 'iyseparatrix2'],
+                ['ix_plate1', 'ix_cut1', '_FILLER_', 'ix_cut2', 'ix_plate2'],
+                ['iyseparatrix3', 'iyseparatrix4'],
+                ['ix_plate3', 'ix_cut3', '_FILLER_', 'ix_cut4', 'ix_plate4']
+            ]
+
+            header = ''
+            for header_items in header_rows:
+                for item in header_items:
+                    header += '{}'.format(gridue[item]).rjust(4)
+                header += '\n'
+            return header
+
+        def format_body(data):
+
+            delim_val = 0
+            delim_char = ''
+            body = ''
+
+            for n in range(5):
+                for j in range(len(data[0])):
+                    for i in range(len(data)):
+                        delim_val += 1
+                        val = np.format_float_scientific(data[i][j][n], precision=15, unique=False).rjust(23).replace('e', 'D')
+                        if delim_val == 3:
+                            delim_val = 0
+                            delim_char = '\n'
+                        body += val + delim_char
+                        delim_char = ''
+
+            if delim_val % 3 != 0:
+                body += '\n'
+
+            return body
+
+        f = open(fname, mode='w')
+        f.write(format_header(gridue_settings) + '\n')
+
+        body_items = ['rm', 'zm', 'psi', 'br', 'bz', 'bpol', 'bphi', 'b']
+        for item in body_items:
+            f.write(format_body(gridue_settings[item]) + '\n')
+
+        runidg = 'iogridue'
+        f.write(runidg + '\n')
+
+        f.close()
+        return True
diff --git a/INGRID/utils.py b/INGRID/utils.py
index 569119a..91367a7 100644
--- a/INGRID/utils.py
+++ b/INGRID/utils.py
@@ -15,9 +15,10 @@
 except:
     pass
 import matplotlib.pyplot as plt
-import pathlib
 import inspect
-from scipy.optimize import root, minimize
+from scipy.optimize import root
+
+
 
 import yaml as yml
 import os
@@ -28,8 +29,9 @@
 from omfit_eqdsk import OMFITgeqdsk
 from INGRID.interpol import EfitData
 from INGRID.line_tracing import LineTracing
-from INGRID.geometry import Point, Line, Patch, segment_intersect, orientation_between
-
+from INGRID.DistribFunc import DistribFunc
+from INGRID.geometry import Point, Line, Patch, orientation_between
+import types
 
 class IngridUtils():
     """
@@ -1070,137 +1072,9 @@ def ClassifyTopology(self, visual=False) -> None:
 
         self.config = self.LineTracer.config
 
-    def WriteGridueSNL(self, gridue_settings: dict, fname: str = 'gridue') -> bool:
-        """
-        Write a gridue file for a single-null configuration.
-
-        Parameters
-        ----------
-        gridue_settings : dict
-            A dictionary containing grid data to be written to the gridue file.
-
-        fname : str, optional
-            The file name/path to save the gridue file to.
-            Defaults to 'gridue'.
-
-        Returns
-        -------
-            True if file was written with no errors
-        """
-
-        def format_header(gridue):
-            header_items = ['nxm', 'nym', 'ixpt1', 'ixpt2', 'iyseptrx1']
-            header = ''
-            for item in header_items:
-                header += '{}'.format(gridue[item]).rjust(4)
-
-            header += '\n'
-            return header
-
-        def format_body(data):
-
-            delim_val = 0
-            delim_char = ''
-            body = ''
-
-            for n in range(5):
-                for j in range(len(data[0])):
-                    for i in range(len(data)):
-                        delim_val += 1
-                        val = np.format_float_scientific(data[i][j][n], precision=15, unique=False).rjust(23).replace('e', 'D')
-                        if delim_val == 3:
-                            delim_val = 0
-                            delim_char = '\n'
-                        body += val + delim_char
-                        delim_char = ''
+    
 
-            if delim_val % 3 != 0:
-                body += '\n'
-
-            return body
-
-        f = open(fname, mode='w')
-        f.write(format_header(gridue_settings) + '\n')
-
-        body_items = ['rm', 'zm', 'psi', 'br', 'bz', 'bpol', 'bphi', 'b']
-        for item in body_items:
-            f.write(format_body(gridue_settings[item]) + '\n')
-
-        runidg = 'iogridue'
-        f.write(runidg + '\n')
-
-        f.close()
-
-        return True
-
-    def WriteGridueDNL(self, gridue_settings: dict, fname: str = 'gridue') -> bool:
-        """
-        Write a gridue file for a double-null configuration.
-
-        Parameters
-        ----------
-        gridue_settings : dict
-            A dictionary containing grid data to be written to the gridue file.
-
-        fname : str, optional
-            The file name/path to save the gridue file to.
-            Defaults to 'gridue'.
-
-        Returns
-        -------
-            True if file was written with no errors
-        """
-
-        def format_header(gridue):
-            header_rows = [
-                ['nxm', 'nym'],
-                ['iyseparatrix1', 'iyseparatrix2'],
-                ['ix_plate1', 'ix_cut1', '_FILLER_', 'ix_cut2', 'ix_plate2'],
-                ['iyseparatrix3', 'iyseparatrix4'],
-                ['ix_plate3', 'ix_cut3', '_FILLER_', 'ix_cut4', 'ix_plate4']
-            ]
-
-            header = ''
-            for header_items in header_rows:
-                for item in header_items:
-                    header += '{}'.format(gridue[item]).rjust(4)
-                header += '\n'
-            return header
-
-        def format_body(data):
-
-            delim_val = 0
-            delim_char = ''
-            body = ''
-
-            for n in range(5):
-                for j in range(len(data[0])):
-                    for i in range(len(data)):
-                        delim_val += 1
-                        val = np.format_float_scientific(data[i][j][n], precision=15, unique=False).rjust(23).replace('e', 'D')
-                        if delim_val == 3:
-                            delim_val = 0
-                            delim_char = '\n'
-                        body += val + delim_char
-                        delim_char = ''
-
-            if delim_val % 3 != 0:
-                body += '\n'
-
-            return body
-
-        f = open(fname, mode='w')
-        f.write(format_header(gridue_settings) + '\n')
-
-        body_items = ['rm', 'zm', 'psi', 'br', 'bz', 'bpol', 'bphi', 'b']
-        for item in body_items:
-            f.write(format_body(gridue_settings[item]) + '\n')
-
-        runidg = 'iogridue'
-        f.write(runidg + '\n')
-
-        f.close()
-        return True
+    
 
     def ReconstructPatches(self, fname: str) -> tuple:
         """
@@ -1227,7 +1101,7 @@ def ReconstructPatches(self, fname: str) -> tuple:
         patches = {}
 
         for raw_patch in data[2:]:
-            patch_data, cell_data, patch_settings = raw_patch
+            patch_data, cell_data, patch_settings,Vertices,cell_grid = raw_patch
             NR = patch_data[0]
             NZ = patch_data[1]
             ER = patch_data[2]
@@ -1246,7 +1120,8 @@ def ReconstructPatches(self, fname: str) -> tuple:
                           plate_patch=patch_settings['plate_patch'],
                           plate_location=patch_settings['plate_location'],
                           PatchTagMap=patch_settings['PatchTagMap'])
-
+            patch.SetVertices(Vertices)
+            patch.cell_grid=cell_grid
             patches[patch.patch_name] = patch
 
         patches = OrderedDict([(k, v) for k, v in patches.items()])
@@ -1267,16 +1142,16 @@ def CheckPatches(self, verbose: bool = False) -> None:
         """
         self.CurrentTopology.CheckPatches(verbose=verbose)
 
-    def PrepGridue(self, guard_cell_eps=1e-3) -> None:
-        """
-        Prepare the gridue for writing.
+    # def PrepGridue(self, guard_cell_eps=1e-3) -> None:
+    #     """
+    #     Prepare the gridue for writing.
 
-        This method calls topology specific implementations of methods that
-        concatenate the Patch object subgrids into a global grid.
-        """
-        self.CurrentTopology.SetupPatchMatrix()
-        self.CurrentTopology.concat_grid(guard_cell_eps=guard_cell_eps)
-        self.CurrentTopology.set_gridue()
+    #     This method calls topology specific implementations of methods that
+    #     concatenate the Patch object subgrids into a global grid.
+    #     """
+    #     self.CurrentTopology.SetupPatchMatrix()
+    #     self.CurrentTopology.concat_grid(guard_cell_eps=guard_cell_eps)
+    #     self.CurrentTopology.set_gridue()
 
     @classmethod
     def _CheckOverlapCells(Grid, Verbose=False):
@@ -1359,6 +1234,7 @@ def __init__(self, Ingrid_obj: object, config: str):
         self.PlateData = Ingrid_obj.PlateData
         self.PatchTagMap = self.parent.GetPatchTagMap(config=config)
         self.LineTracer = Ingrid_obj.LineTracer
+        self.OMFIT_psi = Ingrid_obj.OMFIT_psi
         self.PsiUNorm = Ingrid_obj.PsiUNorm
         self.PsiNorm = Ingrid_obj.PsiNorm
         self.CurrentListPatch = {}
@@ -1404,6 +1280,7 @@ def patch_diagram(self, fig: object = None, ax: object = None) -> None:
         for i, patch in enumerate(self.patches.values()):
             patch.plot_border(color='black', ax=a)
             patch.fill(colors[patch.get_tag()[0]], ax=a, alpha=alpha[patch.get_tag()[-1]])
+            patch.alpha=alpha[patch.get_tag()[-1]]
             patch.color = colors[patch.get_tag()[0]]
         handles, labels = a.get_legend_handles_labels()
         lookup = {label: handle for label, handle in zip(labels, handles)}
@@ -1412,7 +1289,7 @@ def patch_diagram(self, fig: object = None, ax: object = None) -> None:
                  ncol=1)
         f.show()
 
-    def grid_diagram(self, fig: object = None, ax: object = None) -> None:
+    def grid_diagram(self, fig: object = None, ax: object = None,color=None,**kwargs) -> None:
         """
         Generates the grid diagram for a given configuration.
 
@@ -1441,7 +1318,8 @@ def grid_diagram(self, fig: object = None, ax: object = None) -> None:
         ax.set_title(f'{self.config} Subgrid')
 
         for patch in self.patches.values():
-            patch.plot_subgrid(ax=ax)
+            #patch.plot_subgrid(ax=ax)
+            patch.PlotSubGrid(ax=ax,color=color)
             print(f'# Plotting subgrid {patch.patch_name}')
 
         fig.show()
@@ -1480,96 +1358,37 @@ def get_config(self) -> str:
         """
         return self.config
 
-    def get_func(self, func: str) -> 'function':
-        """
-        Create a function from a string input.
-
-        Will be used to generate a poloidal or radial transformation
-        function.
-
-        Parameters
-        ----------
-        func : str
-            An expression to generate a function from.
-
-        Returns
-        -------
-            A function generated from the str input.
-
-        Examples
-        --------
-        When calling method ``get_func`` the user must provide a **string**
-        with the following general format:
-
-        .. math::
-
-            x, f(x)
-
-        That is, the dependent variable and expression to evaluate are
-        separated with a comma character.
-
-        The `Sympy` library supports most expressions that can be generated with
-        Python natively. See the `Sympy` documentation for advanced features.
-
-        Defining a function representing f(x) = x ^ 2:
-
-        >>> func = 'x, x ** 2'
-        >>> f = MyTopologyUtils.get_func(func)
-        >>> f(np.array([0, 1, 2, 3]))
-        array([0, 1, 4, 9])
-
-        Defining a function representing f(x) = exp(x)
-
-        >>> func = 'x, exp(x)'
-        >>> f = MyTopologyUtils.get_func(func)
-        >>> f(np.array([0, 1, 2, 3]))
-        array([ 1.        ,  2.71828183,  7.3890561 , 20.08553692])
-
-        """
-
-        def make_sympy_func(var, expression):
-            import sympy as sp
-            _f = sp.lambdify(var, expression, 'numpy')
-            return _f
-
-        f_str_raw = func
-
-        f_str_raw = f_str_raw.replace(' ', '')
-        delim = f_str_raw.index(',')
-
-        var = f_str_raw[0: delim]
-        expression = f_str_raw[delim + 1:]
-
-        func = make_sympy_func(var, expression)
-        # TODO: Check Normalization of the function to 1
-        return func
-
-    def CheckFunction(self, expression: str, Verbose: bool = False) -> bool:
-        """
-        Check if a str is in the correct format for method ``get_func``
-
-        Parameters
-        ----------
-        expression : str
-            Expression to check.
-
-        Verbose : bool, optional
-            Print full output to terminal. Default to False
-
-        Returns
-        -------
-            True if expression is valid. False otherwise.
-        """
-        ExpValid = False
-        try:
-            com = 'lambda {}:{}'.format(expression.split(',')[0], expression.split(',')[1])
-            if Verbose:
-                print(com)
-            eval(com)
-            ExpValid = True
-        except:
-            raise ValueError('Unable to parse expression entry "{}".'.format(expression))
-        return ExpValid
+    
+    def SetPatchesOrientation(self):
+        
+        
+        # set radial orientation negative if patch is in PFR
+        for name, patch in self.patches.items():
+            if  patch.plate_patch and self.GetPatchLetterNumber(patch.get_tag())[1]=='1':
+                patch.RadOrientation=-1
+            else:
+                patch.RadOrientation=1
+        
+        #set poloidal orientation negative if patches in high field side region (must be checked for other config)
+        Letters=list(dict.fromkeys([L for L,N in  self.GetPatchLetterNumber(self.GetListAllPatches())]))
+        Nhalf=int(len(Letters)/2)
+        ListLFS=Letters[0:Nhalf]
+        for name, patch in self.patches.items():
+            if self.GetPatchLetterNumber(patch.get_tag())[0] in ListLFS:
+                patch.PolOrientation=-1
+            else:
+                patch.PolOrientation=1
+        
+        
+        
+        
+                
+        
+    
+    
+    def ResetHeavyAttr(self):
+        for p in self.patches.values():
+            p.ResetSplines()
 
     def GetFunctions(self, Patch: Patch, Verbose: bool = False) -> tuple:
         """
@@ -1592,49 +1411,43 @@ def GetFunctions(self, Patch: Patch, Verbose: bool = False) -> tuple:
         -------
             2-tuple containing functions for radial and poloidal direction respectively.
         """
-
+        
         poloidal_tag, radial_tag = Patch.get_tag()
         p_f = 'poloidal_f_' + poloidal_tag
         r_f = 'radial_f_' + radial_tag
-
-        try:
-            _poloidal_f = self.settings['grid_settings']['grid_generation'][p_f]
-            valid_function = self.CheckFunction(_poloidal_f, Verbose)
-            if valid_function:
-                _poloidal_f = self.get_func(_poloidal_f)
-            else:
-                raise ValueError('# Invalid function entry. Applying default poloidal function.')
-        except:
-            _poloidal_f = self.settings['grid_settings']['grid_generation']['poloidal_f_default']
-            valid_function = self.CheckFunction(_poloidal_f, Verbose)
-            if valid_function:
-                _poloidal_f = self.get_func(_poloidal_f)
-            else:
-                _poloidal_f = lambda x: x
-
-        try:
-
-            # Adding CORE radial_f support for SNL cases via entry 'radial_f_3'
-            if self.config in ['USN', 'LSN'] \
-                and self.settings['grid_settings']['grid_generation'].get('radial_f_3') is not None \
-                    and poloidal_tag + radial_tag in ['B1', 'C1', 'D1', 'E1']:
-                _radial_f = self.settings['grid_settings']['grid_generation']['radial_f_3']
-            else:
-                _radial_f = self.settings['grid_settings']['grid_generation'][r_f]
-            valid_function = self.CheckFunction(_radial_f, Verbose)
-            if valid_function:
-                _radial_f = self.get_func(_radial_f)
-            else:
-                raise ValueError('# Invalid function entry. Applying default radial function.')
-        except:
-            _radial_f = self.settings['grid_settings']['grid_generation']['radial_f_default']
-            valid_function = self.CheckFunction(_radial_f, Verbose)
-            if valid_function:
-                _radial_f = self.get_func(_radial_f)
-            else:
-                _radial_f = lambda x: x
-
-        return (_radial_f, _poloidal_f)
+        
+        _poloidal_f = DistribFunc(self.settings['grid_settings']['grid_generation'].get(p_f))
+        if _poloidal_f is None:
+            _poloidal_f = DistribFunc(self.settings['grid_settings']['grid_generation'].get('poloidal_f_default'))
+        if _poloidal_f is None:
+            _poloidal_f = DistribFunc('x:x')
+
+    
+        (L,N)=self.GetPatchLetterNumber(Patch.get_tag())
+        # Adding CORE radial_f support for SNL cases via entry 'radial_f_3'
+        if self.config in ['USN', 'LSN'] and N=='1' and Patch.plate_patch:
+            _radial_f = DistribFunc(self.settings['grid_settings']['grid_generation'].get('radial_f_3'))
+        else:
+            _radial_f = DistribFunc(self.settings['grid_settings']['grid_generation'].get(r_f))
+        if _radial_f is None:
+            _radial_f = DistribFunc(self.settings['grid_settings']['grid_generation'].get('radial_f_default'))
+        if _radial_f is None:
+            _radial_f = lambda x: x
+            
+        if hasattr(Patch,'RadOrientation') and Patch.RadOrientation==-1:
+            _radial_f_=lambda x:1-_radial_f(1-x)
+        else:
+            _radial_f_=_radial_f
+        if hasattr(Patch,'PolOrientation') and Patch.PolOrientation==-1:
+            _poloidal_f_=lambda x:1-_poloidal_f(1-x)
+        else:
+            _poloidal_f_=_poloidal_f
+            
+            
+        
+        return (_radial_f_, _poloidal_f_)
+    
+    
 
     def GetNpoints(self, Patch: Patch) -> tuple:
         """
@@ -1680,15 +1493,32 @@ def GetDistortionCorrectionSettings(self) -> dict:
             A dictionary containing ``CorrectDistortion`` settings.
         """
         if self.settings['grid_settings']['grid_generation'].get('distortion_correction') is not None:
-            CD = self.settings['grid_settings']['grid_generation']['distortion_correction']
+            self.distortion_correction = self.settings['grid_settings']['grid_generation']['distortion_correction']
         else:
-            CD = {}
-
-        self.distortion_correction = CD
-        return CD
-
+            self.distortion_correction = {}
+
+        if hasattr(self,'patches'):
+            for name, patch in self.patches.items():
+                (L,N)=self.GetPatchLetterNumber(patch.get_tag())
+                if self.distortion_correction.get(name) is not None:
+                    patch.distortion_correction = self.distortion_correction.get(name)
+                elif self.distortion_correction.get(patch.get_tag()) is not None:
+                    patch.distortion_correction = self.distortion_correction.get(patch.get_tag())
+                elif self.distortion_correction.get(N) is not None:
+                    patch.distortion_correction = self.distortion_correction.get(N)
+                elif self.distortion_correction.get(L) is not None:
+                    patch.distortion_correction = self.distortion_correction.get(L)    
+                elif self.distortion_correction.get('all') is not None:
+                    patch.distortion_correction = self.distortion_correction.get('all')
+                else:
+                    patch.distortion_correction = {'active': False}
+        
+
+    def GetListAllPatches(self):
+        return [v.get_tag() for k,v in self.patches.items()]
+    
     def construct_grid(self, np_cells: int = 1, nr_cells: int = 1, Verbose: bool = False,
-                       ShowVertices: bool = False, RestartScratch: bool = False, ListPatches: str = 'all') -> None:
+                       ShowVertices: bool = False, RestartScratch: bool = False, ListPatches: str = 'all',Visual=True) -> None:
         """
         Construct a grid by refining a Patch map.
 
@@ -1756,45 +1586,100 @@ def construct_grid(self, np_cells: int = 1, nr_cells: int = 1, Verbose: bool = F
             visual = self.settings['DEBUG']['visual']['subgrid']
         except:
             visual = False
+        if Visual:
+            visual=True
         try:
             verbose = self.settings['DEBUG']['verbose']['grid_generation']
         except:
             verbose = False
 
-        verbose = Verbose or verbose
-
+        if Verbose:
+            verbose=True
+        ListAllPatches=self.GetListAllPatches()
+        print('>>> Existing Patches :',ListAllPatches)
+        ListPatches=self.ProcessListPatches(ListPatches)
+        print('>>> Generating grids for patches:',ListPatches)
+        self.SetPatchesOrientation()
         self.GetDistortionCorrectionSettings()
-
-        print('>>> Patches:', [k for k in self.patches.keys()])
+        
+        
         if RestartScratch:
             self.CurrentListPatch = {}
 
         for name, patch in self.patches.items():
-
-            if self.distortion_correction.get(name) is not None:
-                patch.distortion_correction = self.distortion_correction.get(name)
-            elif self.distortion_correction.get(patch.get_tag()) is not None:
-                patch.distortion_correction = self.distortion_correction.get(patch.get_tag())
-            elif self.distortion_correction.get('all') is not None:
-                patch.distortion_correction = self.distortion_correction.get('all')
-            else:
-                patch.distortion_correction = {'Active': False}
-            if (ListPatches == 'all' and patch not in self.CurrentListPatch) or (ListPatches != 'all' and name in ListPatches):
-                self.SetPatchBoundaryPoints(patch, verbose)
+            self.SetPatchBoundaryPoints(patch)
+            
+        for name, patch in self.patches.items():
+            if (ListPatches == 'all' and patch not in self.CurrentListPatch) or (ListPatches != 'all' and patch.get_tag() in ListPatches):
+                self.SetPatchBoundaryPoints(patch)
                 (nr_cells, np_cells) = self.GetNpoints(patch)
                 (_radial_f, _poloidal_f) = self.GetFunctions(patch)
-                print(f'>>> Making subgrid in patch {name}:')
-                print(f'    np = {np_cells}, nr = {nr_cells}')
-                print(f'    fp = {inspect.getsource(_poloidal_f)}')
-                print(f'    fr = {inspect.getsource(_radial_f)}', end='')
+                print('>>> Making subgrid in patch:{} with np={},nr={},fp={},fr={}'.format(patch.get_tag(), np_cells, nr_cells, inspect.getsource(_poloidal_f), inspect.getsource(_radial_f)))
                 patch.RemoveDuplicatePoints()
                 patch.make_subgrid(self, np_cells, nr_cells, _poloidal_f=_poloidal_f, _radial_f=_radial_f, verbose=verbose, visual=visual, ShowVertices=ShowVertices)
 
+                
                 self.CurrentListPatch[name] = patch
-                print(f'    {name} subgrid complete.\n\n')
+                
         self.AdjustGrid()
-
-    def SetPatchBoundaryPoints(self, Patch: Patch, verbose: bool = False) -> None:
+    
+    @staticmethod
+    def GetPatchLetterNumber(PatchName:str or str)->tuple:
+        
+        if type(PatchName)==str:
+            if len(PatchName)<1 or len(PatchName)>2:
+                raise IOError('PatchName must be a string of one or two characters')
+            L=''
+            N=None
+            
+            for C in PatchName:
+                if C.isnumeric():
+                    N=C
+                elif C.isalpha():
+                    L=C
+                else:
+                    raise IOError('PatchName must be a string of one or two characters with one letter max and one digit max')
+            return (L,N) 
+        
+        elif type(PatchName)==list:
+            Out=[]
+            for P in PatchName:
+                Out.append(TopologyUtils.GetPatchLetterNumber(P))
+            return Out
+        else:
+            raise ValueError('PatchName must be a list of a str')
+        
+    def ProcessListPatches(self,ListPatches:list or str):
+        if type(ListPatches)==str:
+            ListPatches=[ListPatches]
+        OrderedList=self.GetListAllPatches()
+  
+        Numbers=list(dict.fromkeys([N for L,N in  self.GetPatchLetterNumber(OrderedList) if N is not None]))
+        Letters=list(dict.fromkeys([L for L,N in  self.GetPatchLetterNumber(OrderedList) if L!='']))
+        if self.Verbose:
+            print('Numbers:',Numbers)
+            print('Letters:',Letters)
+        if 'all' in ListPatches:
+            return OrderedList
+        else:
+            
+            OutList=ListPatches.copy()
+            for P in ListPatches:
+                (L,N)=self.GetPatchLetterNumber(P)
+                if L=='' and N is None:
+                    raise ValueError('Incorrect patch name pattern:{}'.format(P))
+                elif N is None:   
+                    OutList.remove(P)
+                    OutList.extend([L+Num for Num in Numbers])
+                elif L=='':
+                    OutList.remove(P)
+                    OutList.extend([Let+N for Let in Letters])
+            # now let order them 
+            OutList=[P for P in OrderedList if P in list(dict.fromkeys(OutList))]
+            return OutList
+        
+
+    def SetPatchBoundaryPoints(self, Patch: Patch) -> None:
         """
         Set the Patch ``BoundaryPoints`` dict based off TopologyUtils ``ConnexionMap``.
 
@@ -1802,8 +1687,6 @@ def SetPatchBoundaryPoints(self, Patch: Patch, verbose: bool = False) -> None:
         ----------
         Patch : Patch
             The Patch to set the boundary points for.
-        verbose: bool
-            Print full output to terminal.
 
         Notes
         -----
@@ -1812,12 +1695,12 @@ def SetPatchBoundaryPoints(self, Patch: Patch, verbose: bool = False) -> None:
         """
         Patch.TerminatesLoop = False
         if self.ConnexionMap.get(Patch.get_tag()) is not None:
-            if verbose:
+            if self.Verbose:
                 print('Find connexion map for patch {}'.format(Patch.patch_name))
             for v in self.ConnexionMap.get(Patch.get_tag()).items():
                 Boundary, AdjacentPatch = v
                 Patch.BoundaryPoints[Boundary] = self.GetBoundaryPoints(AdjacentPatch)
-                if verbose:
+                if self.Verbose:
                     print('Find Boundaries points for {}'.format(Patch.patch_name))
             if self.ConnexionMap.get(Patch.get_tag()).get('E') is not None:
                 Patch.TerminatesLoop = True
@@ -2189,3 +2072,18 @@ def set_guard(cell_map, ix, iy, eps, boundary):
 
             if debug:
                 self._animate_grid()
+                
+    def SetBoundarySplines(self):
+        for n,p in self.patches.items():
+            p.CreateBoundarySplines(self)
+    
+    def PrepGridue(self, guard_cell_eps=1e-3,**kwargs) -> None:
+        """
+        Prepare the gridue for writing.
+
+        This method calls topology specific implementations of methods that
+        concatenate the Patch object subgrids into a global grid.
+        """
+        self.SetupPatchMatrix()
+        self.concat_grid(guard_cell_eps=guard_cell_eps)
+        self.set_gridue()
diff --git a/build/lib/INGRID/geometry.py b/build/lib/INGRID/geometry.py
index 2fc6e7d..03d112e 100644
--- a/build/lib/INGRID/geometry.py
+++ b/build/lib/INGRID/geometry.py
@@ -504,6 +504,19 @@ def plot_border(self, color: str = 'red', ax: 'matplotlib.axes.Axes' = None) ->
                 [self.vertices['SW'].y, self.vertices['NW'].y],
                 linewidth=1, color=color, label='cell')
 
+    def CollectBorder(self):
+        return np.array([(self.vertices[P].x, self.vertices[P].y) for P in ['NE','NW','SW','SE','NE']])
+
+    def CollectEdges(self):
+        Line1=[(self.vertices['NW'].x, self.vertices['NW'].y),(self.vertices['NE'].x, self.vertices['NE'].y)]
+        Line2=[(self.vertices['NE'].x, self.vertices['NE'].y),(self.vertices['SE'].x, self.vertices['SE'].y)]
+        Line3=[(self.vertices['SE'].x, self.vertices['SE'].y),(self.vertices['SW'].x, self.vertices['SW'].y)]
+        Line4=[(self.vertices['SW'].x, self.vertices['SW'].y),(self.vertices['NW'].x, self.vertices['NW'].y)]
+        return [Line1,Line2,Line3,Line4]
+    
+    
+
+        
     def plot_center(self, color='black', ax: 'matplotlib.axes.Axes' = None) -> None:
         """
         Plot the center of a Cell.
@@ -565,8 +578,16 @@ def __init__(self, lines: 'array-like', patch_name: str = '', PatchTagMap: dict
         self.E = lines[1]
         self.S = lines[2]
         self.W = lines[3]
+        self.E_vertices=[]
+        self.W_vertices=[]
+        self.N_vertices=[]
+        self.S_vertices=[]
         self.BoundaryPoints = {}
-
+        self.N_spl=None
+        self.E_spl=None
+        self.W_spl=None
+        self.N_spl=None
+        self.radial_spl=[]
         # This is the border for the fill function
         # It need to only include N and S lines
         self.RemoveDuplicatePoints()
@@ -576,7 +597,10 @@ def __init__(self, lines: 'array-like', patch_name: str = '', PatchTagMap: dict
         self.plate_location = plate_location
         self.patch_name = patch_name
         self.color = color
+        self.alpha=1
         self.Verbose = True
+        self.RadOrientation=1
+        self.PolOrientation=1
 
         # TODO: Populate these attributes when generating subgrid.
         #       Currently this is being done in the concat_grid method.
@@ -595,10 +619,17 @@ def __init__(self, lines: 'array-like', patch_name: str = '', PatchTagMap: dict
             'S': 'SOL',
             'C': 'Core'
         }
-
-        self.ax = None
         self.cell_grid = None
-
+    def ResetSplines(self):
+        self.radial_spl=[]
+        
+    def GetVertices(self):
+        return [getattr(self,'{}_vertices'.format(O)) for O in ['N','S','E','W']]
+    
+    def SetVertices(self,VerticesList):
+        for V,O in zip(VerticesList,['N','S','E','W']):
+            setattr(self,'{}_vertices'.format(O),V)
+            
     def plot_border(self, color='red', ax=None):
         """
         Draw solid borders around the patch.
@@ -625,7 +656,8 @@ def fill(self, color='lightsalmon', ax=None, alpha=1.0):
         x = np.array([p.x for p in self.p])
         y = np.array([p.y for p in self.p])
         arr = np.column_stack((x, y))
-        patch = Polygon(arr, fill=True, closed=True, color=color, label=self.get_tag(), alpha=alpha)
+        PatchLabel = self.patch_name + ' (' + UnfoldLabel(self.PatchLabelDoc, self.patch_name) + ')'
+        patch = Polygon(arr, fill=True, closed=True, color=color, label=PatchLabel, alpha=alpha)
         _ax = plt.gca() if ax is None else ax
         _ax.add_patch(patch)
         _ax.plot()
@@ -636,8 +668,71 @@ def plot_subgrid(self, ax=None, color='blue'):
         for row in self.cell_grid:
             for cell in row:
                 cell.plot_border(color=color, ax=ax)
-        plt.pause(0.1)
-
+        #plt.pause(0.1)
+        
+    def PlotSubGrid(self,ax=None,color=None):
+        if ax is None:
+            ax = plt.gca()
+        ListBorder=[]
+        ListEdges=[]
+        DataEdges=[]
+        
+        for row in self.cell_grid:
+            for cell in row:
+                Border=cell.CollectBorder()
+                
+                p=matplotlib.patches.Polygon(Border,edgecolor=color,fill=False)
+                ListBorder.append(p)
+                Edges=cell.CollectEdges()
+                ListEdges.extend(Edges)
+                for E in Edges:
+                    DataEdges.append(GetLength(E))
+                  
+
+        #Collection border of cell
+        C=matplotlib.collections.PatchCollection(ListBorder)
+        C.set_picker(True)
+        C.set_facecolor('')
+        if color is None:
+            C.set_edgecolor(self.color)
+        else:
+            C.set_edgecolor(color)
+        #ax.add_collection(C) 
+        
+        #Collection edges of cells
+        CE=matplotlib.collections.LineCollection(ListEdges,array=np.array(range(0,len(ListEdges))))
+        CE.set_picker(1)
+        ax.add_collection(CE) 
+        if color is None:
+            CE.set_color(self.color)
+            CE.set_edgecolors(self.color)
+            CE.set_facecolors(self.color)
+        else:
+            CE.set_color(color)
+            CE.set_edgecolors(color)
+            CE.set_facecolors(color)
+            
+        CE.set_alpha(self.alpha)
+        
+            
+        def onpick(evt):           
+            if evt.artist==CE:
+                annot.set_visible(False)
+                annot.xy = (evt.mouseevent.xdata, evt.mouseevent.ydata)
+                annot.set_text('length={:.2e} m'.format(DataEdges[evt.artist.get_array()[evt.ind[0]]]))
+                annot.set_visible(True)
+            if evt.mouseevent.button == 3:
+                annot.set_visible(False)  
+            ax.figure.canvas.draw_idle()                 
+        annot = ax.annotate("", xy=(0,0), xytext=(-20,20),textcoords="offset points",
+                                bbox=dict(boxstyle="round", fc="w"),
+                                arrowprops=dict(arrowstyle="->"))
+        annot.set_visible(False)
+        ax.figure.canvas.mpl_connect('pick_event', onpick)
+        
+                             
+        ax.autoscale_view()
+                
     def RemoveDuplicatePoints(self):
         for line in [self.N, self.E, self.S, self.W]:
             line.RemoveDuplicatePoints()
@@ -711,95 +806,29 @@ def as_np(self):
             patch_data.append(Z)
         patch_data = np.array(patch_data)
         cell_data = self.cell_grid_as_np()
+        cell_grid = self.cell_grid
         patch_settings = self.get_settings()
-        return np.array([patch_data, cell_data, patch_settings])
-
-    def make_subgrid(self, grid, np_cells=2, nr_cells=2, _poloidal_f=lambda x: x, _radial_f=lambda x: x, verbose=False, visual=False, ShowVertices=False):
-        """
-        Generate a refined grid within a patch.
-        This 'refined-grid' within a Patch is a collection
-        of num x num Cell objects
-
-        Parameters
-        ----------
-        grid : Ingrid
-                To be used for obtaining Efit data and all
-                other class information.
-        num  : int, optional
-                Number to be used to generate num x num
-                cells within our Patch.
-        """
-
-        def psi_parameterize(grid, r, z):
-            """
-            Helper function to be used to generate a
-            list of values to parameterize a spline
-            in Psi. Returns a list to be used for splprep only
-            """
-            vmax = grid.PsiNorm.get_psi(r[-1], z[-1])
-            vmin = grid.PsiNorm.get_psi(r[0], z[0])
-
-            vparameterization = np.empty(0)
-            for i in range(len(r)):
-                vcurr = grid.PsiNorm.get_psi(r[i], z[i])
-                vparameterization = np.append(vparameterization, abs((vcurr - vmin) / (vmax - vmin)))
-
-            return vparameterization
-
-        def psi_test(leg, grid):
-            """
-            Determine whether a leg's orientation is increasing in psi
-            by checking the endpoints of a given leg.
-
-            Returns True if final Point entry is greater in psi than the
-            first Point entry. Returns False otherwise.
-            """
-            p1 = leg.p[0]
-            p2 = leg.p[-1]
-
-            return True if grid.PsiNorm.get_psi(p2.x, p2.y) > grid.PsiNorm.get_psi(p1.x, p1.y) else False
-
-        def set_dynamic_step(ratio=0.01):
-            """
-            Compute a dt value for line_tracing that occurs during subgrid generation of
-            this Patch object. The value is taken to be a ratio of the average length of
-            a Patch in the radial direction.
-            """
-            d = 0
-            for i in range(nr_lines):
-                d += np.sqrt((self.E_vertices[i].x - self.W_vertices[i].x)**2 + (self.E_vertices[i].y - self.W_vertices[i].y)**2)
-            dynamic_step = d / nr_lines
-            # print('Dynamic-Step value came out to be: {}\n'.format(dynamic_step * ratio))
-            return dynamic_step * ratio
-
-        # Allocate space for collection of cell objects.
-        # Arbitrary 2D container for now.
-        cell_grid = []
-
-        if verbose:
-            print('Constructing grid for patch "{}" with dimensions (np, nr) = ({}, {})'.format(self.patch_name, np_cells, nr_cells))
-            print(np_cells)
-            print(nr_cells)
-        np_lines = np_cells + 1
-        nr_lines = nr_cells + 1
-        if verbose: print(' # Create B-Splines along the North and South boundaries.')
+        return np.array([patch_data, cell_data, patch_settings,self.GetVertices(),cell_grid])
+    
+    
+    def CreateBoundarySplines(self,grid):
         # Create B-Splines along the North and South boundaries.
+        if self.Verbose: print(' # Create B-Splines along the North and South boundaries.')
         N_vals = self.N.fluff()
-
         self.N_spl, uN = splprep([N_vals[0], N_vals[1]], s=0)
         # Reverse the orientation of the South line to line up with the North.
 
         S_vals = self.S.reverse_copy().fluff()
         self.S_spl, uS = splprep([S_vals[0], S_vals[1]], s=0)
-        if verbose: print(' # Create B-Splines along West boundaries.')
+        if self.Verbose: print(' # Create B-Splines along West boundaries.')
         # Create B-Splines along the East and West boundaries.
         # Parameterize EW splines in Psi
         try:
             #Cannot fluff with too many points
             n = 500 if len(self.W.p) < 50 else 100
            # W_vals = self.W.reverse_copy().fluff(num = n)
-            W_vals = self.W.reverse_copy().fluff(n, verbose=verbose)
-            Psi = psi_parameterize(grid, W_vals[0], W_vals[1])
+            W_vals = self.W.reverse_copy().fluff(n, verbose=self.Verbose)
+            Psi = self.psi_parameterize(grid, W_vals[0], W_vals[1])
             self.W_spl, uW = splprep([W_vals[0], W_vals[1]], u=Psi, s=10)
         except Exception as e:
             exc_type, exc_obj, tb = sys.exc_info()
@@ -810,16 +839,17 @@ def set_dynamic_step(ratio=0.01):
             line = linecache.getline(filename, lineno, f.f_globals)
             print('EXCEPTION IN ({}, LINE {} "{}"): {}'.format(filename, lineno, line.strip(), exc_obj))
 
-        if verbose: print(' # Create B-Splines along the East boundaries.')
+        if self.Verbose: print(' # Create B-Splines along the East boundaries.')
         try:
             n = 500 if len(self.E.p) < 50 else 100
             E_vals = self.E.fluff(num=n)
-            self.E_spl, uE = splprep([E_vals[0], E_vals[1]], u=psi_parameterize(grid, E_vals[0], E_vals[1]), s=10)
+            self.E_spl, uE = splprep([E_vals[0], E_vals[1]], u=self.psi_parameterize(grid, E_vals[0], E_vals[1]), s=10)
         except Exception as e:
             print(' Number of points on the boundary:', len(self.E.p))
             plt.plot(E_vals[0], E_vals[1], '.', color='black')
             print(repr(e))
-        if verbose: print(' #check plate_patch')
+        if self.Verbose: print(' #check plate_patch')
+        # ACCURACY ON PSI_MIN SIDE OF W IDL LINE ISSUE.
 
         if self.plate_patch:
 
@@ -845,7 +875,7 @@ def find_nearest(array, value):
                 plate_north = [N_vals[0][-1], N_vals[1][-1]]
                 plate_south = [S_vals[0][-1], S_vals[1][-1]]
 
-            if verbose:
+            if self.Verbose:
                 print('=' * 80 + '\n')
                 print('Parameterization of Target Plate in PSI:')
                 print(_u)
@@ -860,19 +890,19 @@ def find_nearest(array, value):
                 try:
                     plate_north_index = lookup[find_nearest(_u, fsolve(f, 0, args=(U_spl, plate_north[0], plate_north[1])))]
                 except:
-                    if verbose: print('NorthIndex: ERROR IN PARAMETERIZATION IN PSI')
+                    print('NorthIndex: ERROR IN PARAMETERIZATION IN PSI')
             try:
                 plate_south_index = lookup[find_nearest(_u, brentq(f, _u[0], _u[-1], args=(U_spl, plate_south[0], plate_south[1])))]
             except ValueError:
                 try:
                     plate_south_index = lookup[find_nearest(_u, fsolve(f, 0, args=(U_spl, plate_south[0], plate_south[1])))]
                 except:
-                    if verbose: print('SouthIndex: ERROR IN PARAMETERIZATION IN PSI')
+                    print('SouthIndex: ERROR IN PARAMETERIZATION IN PSI')
             if plate_south_index > plate_north_index:
                 plate_north_index, plate_south_index = plate_south_index, plate_north_index
 
             U_vals = [U_vals[0][plate_south_index:plate_north_index + 1], U_vals[1][plate_south_index:plate_north_index + 1]]
-            U_spl, _u = splprep([U_vals[0], U_vals[1]], u=psi_parameterize(grid, U_vals[0], U_vals[1]), s=0)
+            U_spl, _u = splprep([U_vals[0], U_vals[1]], u=self.psi_parameterize(grid, U_vals[0], U_vals[1]), s=0)
 
             if self.plate_location == 'W':
                 W_vals = U_vals
@@ -882,9 +912,86 @@ def find_nearest(array, value):
                 E_vals = U_vals
                 self.E_spl = U_spl
                 uE = _u
+                
+                
+    @staticmethod
+    def psi_parameterize(grid, r, z):
+            """
+            Helper function to be used to generate a
+            list of values to parameterize a spline
+            in Psi. Returns a list to be used for splprep only
+            """
+            vmax = grid.PsiNorm.get_psi(r[-1], z[-1])
+            vmin = grid.PsiNorm.get_psi(r[0], z[0])
+
+            vparameterization = np.empty(0)
+            for i in range(len(r)):
+                vcurr = grid.PsiNorm.get_psi(r[i], z[i])
+                vparameterization = np.append(vparameterization, abs((vcurr - vmin) / (vmax - vmin)))
+
+            return vparameterization
+        
+    def make_subgrid(self, grid, np_cells=2, nr_cells=2, _poloidal_f=lambda x: x, _radial_f=lambda x: x, verbose=False, visual=False, ShowVertices=False):
+        """
+        Generate a refined grid within a patch.
+        This 'refined-grid' within a Patch is a collection
+        of num x num Cell objects
+
+        Parameters
+        ----------
+        grid : Ingrid
+                To be used for obtaining Efit data and all
+                other class information.
+        num  : int, optional
+                Number to be used to generate num x num
+                cells within our Patch.
+        """
+
+        
+
+        def psi_test(leg, grid):
+            """
+            Determine whether a leg's orientation is increasing in psi
+            by checking the endpoints of a given leg.
+
+            Returns True if final Point entry is greater in psi than the
+            first Point entry. Returns False otherwise.
+            """
+            p1 = leg.p[0]
+            p2 = leg.p[-1]
+
+            return True if grid.PsiNorm.get_psi(p2.x, p2.y) > grid.PsiNorm.get_psi(p1.x, p1.y) else False
+
+        def set_dynamic_step(ratio=0.01):
+            """
+            Compute a dt value for line_tracing that occurs during subgrid generation of
+            this Patch object. The value is taken to be a ratio of the average length of
+            a Patch in the radial direction.
+            """
+            d = 0
+            for i in range(nr_lines):
+                d += np.sqrt((self.E_vertices[i].x - self.W_vertices[i].x)**2 + (self.E_vertices[i].y - self.W_vertices[i].y)**2)
+            dynamic_step = d / nr_lines
+            # print('Dynamic-Step value came out to be: {}\n'.format(dynamic_step * ratio))
+            return dynamic_step * ratio
+
+        # Allocate space for collection of cell objects.
+        # Arbitrary 2D container for now.
+        cell_grid = []
+
+        if verbose:
+            print('Constructing grid for patch "{}" with dimensions (np, nr) = ({}, {})'.format(self.patch_name, np_cells, nr_cells))
+            print(np_cells)
+            print(nr_cells)
+        np_lines = np_cells + 1
+        nr_lines = nr_cells + 1
+        
+        #self.CreatBoundarySpleen()
 
         # Generate our sub-grid anchor points along the North
         # and South boundaries of our patch.
+        self.W_vertices_old=self.W_vertices
+        
         self.N_vertices = []
         self.S_vertices = []
         self.E_vertices = []
@@ -898,7 +1005,7 @@ def find_nearest(array, value):
                 _n = splev(_poloidal_f(i / (np_lines - 1)), self.N_spl)
                 self.N_vertices.append(Point((float(_n[0]), float(_n[1]))))
         else:
-            if verbose: print('Find boundary points at face "N" for {}:{}'.format(self.patch_name, self.BoundaryPoints.get('N')))
+            if self.Verbose: print('Find boundary points at face "N" for {}:{}'.format(self.patch_name, self.BoundaryPoints.get('N')))
             self.N_vertices = self.BoundaryPoints.get('N')
 
         if self.BoundaryPoints.get('S') is None:
@@ -911,15 +1018,15 @@ def find_nearest(array, value):
         u = [_radial_f(i / (nr_lines - 1)) for i in range(nr_lines)]
 
         Npts = 1000
-        xy = splev(np.linspace(0, 1, Npts), self.E_spl)
-
+        #xy = splev(np.linspace(0, 1, Npts), self.E_spl)
+        
         if self.BoundaryPoints.get('W') is None:
             for i in range(nr_lines):
                 _w = splev(u[i], self.W_spl)
                 self.W_vertices.append(Point((float(_w[0]), float(_w[1]))))
         else:
             self.W_vertices = self.BoundaryPoints.get('W')
-
+        
         if self.BoundaryPoints.get('E') is None:
             for i in range(nr_lines):
                 _e = splev(u[i], self.E_spl)
@@ -943,11 +1050,11 @@ def find_nearest(array, value):
 
             else:
                 color = self.color
-                ms = 6
-                markers = ['o', 'X', 's', 'D']
+                ms = 10
+                markers = ['o', 'P', 's', 'D']
             for vertices, mark in zip([self.W_vertices, self.E_vertices, self.N_vertices, self.S_vertices], markers):
                 for p in vertices:
-                    plt.plot(p.x, p.y, '.', color=color, markersize=ms, marker=mark, markeredgecolor='black')
+                    plt.plot(p.x, p.y, '.', color=color, markersize=ms, marker=mark, markeredgecolor=color)
         # Radial lines of Psi surfaces. Ordered with increasing magnitude, starting with
         # the South boundary of the current Patch, and ending with the North boundary of
         # this current Patch. These will serve as delimiters when constructing cells.
@@ -956,32 +1063,38 @@ def find_nearest(array, value):
         dynamic_step = set_dynamic_step()
         if verbose: print('# Interpolate radial lines between North and South patch boundaries..')
         # Interpolate radial lines between North and South patch boundaries.
+        if hasattr(self,'radial_spl'):
+            self.radial_spl_old=self.radial_spl
+        else:
+            self.radial_spl_old=[]
         self.radial_spl = []
         temp_vertices = []
         temp_vertices.append(self.N.p[-1])
         for i in range(len(self.W_vertices) - 2):
             #TODO: parallelize tracing of radial lines (draw_line function must be "externalized" in the scope of the script)
+            # if (len(self.W_vertices)!=len(self.W_vertices_old) or self.W_vertices[i+1]!=self.W_vertices_old[i+1]) or len(self.radial_spl_old)!=len(self.W_vertices)-2:
+    
             radial_lines.append(grid.LineTracer.draw_line(self.W_vertices[i + 1], {'line': self.E}, option='theta',
-                direction='cw', show_plot=visual, dynamic_step=dynamic_step, text=verbose))
+                direction='cw', show_plot=0, dynamic_step=dynamic_step))
             temp_vertices.append(radial_lines[-1].p[-1])
             radial_vals = radial_lines[i + 1].fluff(1000)
             Radial_spl, uR = splprep([radial_vals[0], radial_vals[1]], s=0)
             self.radial_spl.append(Radial_spl)
+            # else:
+            #     self.radial_spl.append(self.radial_spl_old[i])
+                
             vertex_list = []
             for j in range(np_lines):
                 u = _poloidal_f(j / (np_lines - 1))
                 _r = splev(u, self.radial_spl[i])
                 Pt = Point((float(_r[0]), float(_r[1])))
                 if self.distortion_correction['active'] and j > 0 and j < np_lines - 1:
-                    Res = self.distortion_correction['resolution']
-                    ThetaMin = self.distortion_correction['theta_min']
-                    ThetaMax = self.distortion_correction['theta_max']
                     umin = _poloidal_f((j - 1) / (np_lines - 1))
                     umax = _poloidal_f((j + 1) / (np_lines - 1))
                     Pt1 = radial_vertices[i][j]
                     Pt2 = radial_vertices[i][j - 1]
                     Tag = '---- Correcting points: {},{}'.format(i, j)
-                    Pt = CorrectDistortion(u, Pt, Pt1, Pt2, self.radial_spl[i], ThetaMin, ThetaMax, umin, umax, Res, visual, Tag, verbose)
+                    Pt = CorrectDistortion(u, Pt, Pt1, Pt2, self.radial_spl[i],umin,umax,Tag,**self.distortion_correction)
 
                 vertex_list.append(Pt)
                 if visual:
@@ -998,15 +1111,12 @@ def find_nearest(array, value):
             for j in range(1, np_lines - 1):
                 u = _poloidal_f(j / (np_lines - 1))
                 Pt = self.S_vertices[j]
-                Res = self.distortion_correction['resolution']
-                ThetaMin = self.distortion_correction['theta_min']
-                ThetaMax = self.distortion_correction['theta_max']
                 umin = _poloidal_f((j - 1) / (np_lines - 1))
                 umax = _poloidal_f((j + 1) / (np_lines - 1))
                 Pt1 = radial_vertices[-1][j]
                 Pt2 = radial_vertices[-1][j - 1]
                 Tag = '---- Correcting south boundary points:{}'.format(j)
-                self.S_vertices[j] = CorrectDistortion(u, Pt, Pt1, Pt2, self.S_spl, ThetaMin, ThetaMax, umin, umax, Res, visual, Tag, verbose)
+                self.S_vertices[j] = CorrectDistortion(u, Pt, Pt1, Pt2, self.S_spl, umin, umax,  Tag,**self.distortion_correction)
         radial_vertices.append(self.S_vertices)
         if verbose: print('# Create Cells: South boundary -> North Boundary')
         # Create Cells: South boundary -> North Boundary
@@ -1298,34 +1408,36 @@ def trim_geometry(geoline, start, end):
     return trim
 
 
-def CorrectDistortion(u, Pt, Pt1, Pt2, spl, ThetaMin, ThetaMax, umin, umax, Resolution, visual, Tag, MinTol=1.02, MaxTol=0.98, Verbose=False):
-    MaxIter = Resolution * 10
-    dumax = (umax - u) / Resolution
-    dumin = (u - umin) / Resolution
+def CorrectDistortion(u, Pt, Pt1, Pt2, spl, umin, umax , Tag, theta_min=60, theta_max=120,resolution=1000,min_tol=1.05, max_tol=0.95, zero_span=0.90,max_span=1.0,verbose=False,visual=False,**kwargs):
+    MaxIter = resolution * 10
+    dumax = (umax - u) / resolution
+    dumin = (u - umin) / resolution
     Theta = Line([Pt1, Pt]).GetAngle(Line([Pt1, Pt2]))
-    if Verbose: print(f'Current theta value: {Theta}')
-    if Theta < ThetaMin or Theta > ThetaMax:
-        if Verbose: print('{}: u={};Theta={};ThetaMin={};ThetaMax={}'.format(Tag, u, Theta, ThetaMin, ThetaMax))
+    if verbose: print(f'Current theta value: {Theta}')
+    if Theta < theta_min or Theta > theta_max:
+        if verbose: print('{}: u={};Theta={};ThetaMin={};ThetaMax={}'.format(Tag, u, Theta, theta_min, theta_max))
         if visual:
                 plt.plot(Pt.x, Pt.y, '.', color='red', markersize=8, marker='o')
                 plt.show()
                 plt.draw()
         icount = 0
         color = 'purple'
-        while Theta < ThetaMin or Theta > ThetaMax:
+        u0=u
+
+        while Theta < theta_min or Theta > theta_max:
             icount += 1
-            if Theta < ThetaMin:
+            if Theta < theta_min:
                u = u + dumax
-            elif Theta > ThetaMax:
+            elif Theta > theta_max:
                u = u - dumin
-            if (u > umax * MaxTol or u < umin * MinTol) or icount >= MaxIter:
-               if Verbose: print('[{}]>>>> umax={} umin={};u:{};Theta={}'.format(icount, umax, umin, u, Theta))
+            if (u-u0) > (umax* max_tol-u0)*max_span or (u-u0) < (umin * min_tol-u0)*max_span or icount >= MaxIter or (abs(umin)<1e-12 and u<u0*zero_span):
+               if verbose: print('[{}]>>>> umax={} umin={};u:{};Theta={}'.format(icount, umax, umin, u, Theta))
                color = 'gray'
                break
             _r = splev(u, spl)
             Pt = Point((_r[0], _r[1]))
             Theta = Line([Pt1, Pt]).GetAngle(Line([Pt1, Pt2]))
-        if Verbose: print('[{}]>>>> u:{};Theta={}'.format(icount, u, Theta))
+        if verbose: print('[{}]>>>> u:{};Theta={}'.format(icount, u, Theta))
         if visual:
             plt.plot(Pt.x, Pt.y, '.', color=color, markersize=8, marker='s')
     return Pt
@@ -1490,7 +1602,7 @@ def segment_intersect(line1, line2, verbose=False):
 
 def UnfoldLabel(Dic: dict, Name: str) -> str:
     """
-    Unfold Patch label (e.g. "C1" -> "Inner Core Top")
+    Unfold Patch label (e.g. "ICT" -> "Inner Core Top")
 
     Parameters
     ----------
@@ -1515,3 +1627,6 @@ def UnfoldLabel(Dic: dict, Name: str) -> str:
         return ''.join(Output)
     else:
         return ''
+    
+def GetLength(E):
+    return np.sqrt((E[1][0]-E[0][0])**2+(E[1][1]-E[1][1])**2)
diff --git a/build/lib/INGRID/ingrid.py b/build/lib/INGRID/ingrid.py
index e2427d7..c13b861 100644
--- a/build/lib/INGRID/ingrid.py
+++ b/build/lib/INGRID/ingrid.py
@@ -17,18 +17,17 @@
     pass
 import matplotlib.pyplot as plt
 import pathlib
-import inspect
-from scipy.optimize import root, minimize
+
 
 import yaml as yml
 import os
 from pathlib import Path
 from time import time
-from collections import OrderedDict
 
-from omfit_eqdsk import OMFITgeqdsk
+
+from INGRID.DistribFunc import *
 from INGRID.interpol import EfitData
-from INGRID.interpol import Bicubic
+
 from INGRID.utils import IngridUtils
 from INGRID.topologies.snl import SNL
 from INGRID.topologies.sf15 import SF15
@@ -38,8 +37,10 @@
 from INGRID.topologies.sf135 import SF135
 from INGRID.topologies.sf165 import SF165
 from INGRID.topologies.udn import UDN
-from INGRID.line_tracing import LineTracing
 from INGRID.geometry import Point, Line
+import types
+
+
 
 
 def QuickStart() -> None:
@@ -193,6 +194,45 @@ def SaveSettingsFile(self, fname: str = '', settings: dict = {}) -> 'Path':
             yml.dump(self.settings, f)
 
         return fpath
+    
+    
+    def StoreSettings(self, fname: str = '', settings: dict = {}) -> 'Path':
+        """
+        Save a new settings .yml file.
+
+        Parameters
+        ----------
+        fname : optional
+            Name of new settings '.yml' file.
+            If default value of '', then Ingrid will generate
+            a '.yml' file named 'INGRID_Session' appended with
+            a timestamp.
+
+        settings : optional
+            Ingrid settings dictionary to dump into
+            the '.yml' file. Defaults to empty dict which produces a
+            template settings file.
+
+        Returns
+        -------
+            A Path instance representing the saved YAML file.
+        """
+
+        if fname == '':
+            fname = 'INGRID_Session' + str(int(time())) + '.yml'
+
+        fpath = Path(fname)
+
+        if fpath.suffix == '':
+            fname += '.yml'
+        elif fpath.suffix != '.yml':
+            fname = fname[:-4] + '.yml'
+
+        with open(fname, 'w') as f:
+            yml.dump(self.settings, f)
+
+        return fpath
+
 
     def PrintSummaryInput(self) -> None:
         """
@@ -778,8 +818,8 @@ def PlotTargetPlate(self, plate_key: str, color: str = 'red', ax: object = None)
             ax = plt.gca()
 
         if plate_key in [k for k in self.PlateData.keys()]:
-            self.RemovePlotLine(label=plate_key, ax=ax)
             if type(self.PlateData[plate_key]) is Line:
+                self.RemovePlotLine(label=plate_key, ax=ax)
                 self.PlateData[plate_key].plot(label=plate_key, color=color)
             else:
                 pass
@@ -1049,7 +1089,7 @@ def PlotPatches(self) -> None:
         if self.settings['grid_settings']['patch_generation']['strike_pt_loc'] == 'target_plates':
             self.RemovePlotLine(label='limiter', ax=self.PatchAx)
 
-    def PlotSubgrid(self) -> None:
+    def PlotSubgrid(self,NewFig=False,**kwargs) -> None:
         """
         Plot the grid that was generated with method 'CreateSubgrid'.
         """
@@ -1057,9 +1097,13 @@ def PlotSubgrid(self) -> None:
             plt.close(self._SubgridFig)
         except:
             pass
-        self._SubgridFig = plt.figure('INGRID: ' + self.CurrentTopology.config + ' Grid', figsize=(6, 10))
-        self._SubgridAx = self._SubgridFig.add_subplot(111)
-        self.CurrentTopology.grid_diagram(fig=self._SubgridFig, ax=self._SubgridAx)
+        if NewFig:
+            self._SubgridFig = plt.figure('INGRID: ' + self.CurrentTopology.config + ' Grid', figsize=(6, 10))
+            self._SubgridAx = self._SubgridFig.add_subplot(111)
+        else:
+            self._SubgridFig = plt.gcf()
+            self._SubgridAx = plt.gca()
+        self.CurrentTopology.grid_diagram(fig=self._SubgridFig, ax=self._SubgridAx,**kwargs)
 
     def AutoRefineMagAxis(self) -> None:
         """
@@ -1267,11 +1311,12 @@ def ConstructPatches(self) -> None:
         self.OrderTargetPlates()
         self.PrepLineTracing()
         self.CurrentTopology.construct_patches()
-        self.CurrentTopology.GroupPatches()
+        self.CurrentTopology.GroupPatches() 
         self.CheckPatches()
-
-        if self.settings['patch_data']['preferences']['new_file']:
-            self.SavePatches(self.settings['patch_data']['preferences']['new_fname'])
+        self.CurrentTopology.SetBoundarySplines()
+        #not practical in a workflow
+        #if self.settings['patch_data']['preferences']['new_file']:
+        #   self.SavePatches(self.settings['patch_data']['preferences']['new_fname'])
 
     def SavePatches(self, fname: str = '') -> None:
         """
@@ -1298,6 +1343,35 @@ def SavePatches(self, fname: str = '') -> None:
             tail = '.npy'
         print(f"# Saved patch data for file {Path(fname).name + tail}")
 
+        
+    def SuperSavePatches(self, fname: str = '') -> None:
+        """
+        Save patches as '.npy' file for later reconstruction in Ingrid.
+
+        This file contains the information required to reconstruct patches
+        at a later time and bypass the line_tracing.
+
+        Parameters
+        ----------
+        fname : str, optional
+            Name of file/location for patch data.
+        """
+        if fname in ['', None]:
+            fname = self.CurrentTopology.config + '_patches_' + str(int(time()))
+        self.CurrentTopology.ResetHeavyAttr()
+        
+        Data={'patches':self.CurrentTopology.patches,
+              'config':self.CurrentTopology.config,
+              'xpt_data':self.CurrentTopology.LineTracer.NSEW_lookup
+              }
+        np.save(fname, Data)
+        
+        if Path(fname).suffix == '.npy':
+            tail = ''
+        else:
+            tail = '.npy'
+        print(f"# Saved patch data for file {Path(fname).name + tail}")
+
     def LoadPatches(self, fname: str = '') -> None:
         """
         Load patches stored in an Ingrid generated '.npy' file.
@@ -1317,7 +1391,7 @@ def LoadPatches(self, fname: str = '') -> None:
         if fname.strip() == '':  # Check if settings contains patch data.
             fname = self.settings['patch_data']['file']
         config, xpt_data, patches = self.ReconstructPatches(fname)
-        self.PopulateSettings(Ingrid.ReadYamlFile(self.InputFile))
+        #self.PopulateSettings(Ingrid.ReadYamlFile(self.InputFile))
         self.StartSetup()
         self.LineTracer.NSEW_lookup = xpt_data
         self.SetTopology(config)
@@ -1325,8 +1399,37 @@ def LoadPatches(self, fname: str = '') -> None:
         self.CurrentTopology.OrderPatches()
         self.CurrentTopology.SetupPatchMatrix()
         self.CheckPatches()
+        
+    def SuperLoadPatches(self, fname: str = '') -> None:
+        """
+        Load patches stored in an Ingrid generated '.npy' file.
+
+        Parameters
+        ----------
+        fname : str, optional
+            Path to patch data.
+                If no fname is provided to method 'LoadPatches', Ingrid code will check the settings
+                'dict' for a file under entry ``settings['patch_data']['file']``
+
+        """
+
+        if type(fname) is not str:
+            raise ValueError('# User did not provide a string to patch data.')
+
+        if fname.strip() == '':  # Check if settings contains patch data.
+            fname = self.settings['patch_data']['file']
+        Data=np.load(fname,allow_pickle=True).tolist()
+        #config, xpt_data, patches = self.ReconstructPatches(fname)
+        #self.PopulateSettings(Ingrid.ReadYamlFile(self.InputFile))
+        self.StartSetup()
+        self.LineTracer.NSEW_lookup = Data['xpt_data']
+        self.SetTopology(Data['config'])
+        self.CurrentTopology.patches = Data['patches']
+        self.CurrentTopology.OrderPatches()
+        self.CurrentTopology.SetupPatchMatrix()
+        self.CheckPatches()
 
-    def CreateSubgrid(self, NewFig: bool = True, ShowVertices: bool = False) -> None:
+    def CreateSubgrid(self, NewFig: bool = True, ShowVertices: bool = False,**kwargs) -> None:
         """
         Refine a generated patch map into a grid for exporting.
 
@@ -1349,9 +1452,9 @@ def CreateSubgrid(self, NewFig: bool = True, ShowVertices: bool = False) -> None
             the particular patch.
 
         """
-        self.CurrentTopology.construct_grid()
+        self.CurrentTopology.construct_grid( ShowVertices=ShowVertices,**kwargs)
 
-    def ExportGridue(self, fname: str = 'gridue') -> None:
+    def ExportGridue(self, fname: str = 'gridue',**kwargs) -> None:
         """
         Export a gridue file for the created grid.
 
@@ -1361,12 +1464,9 @@ def ExportGridue(self, fname: str = 'gridue') -> None:
             Name of gridue file to save.
 
         """
-        if type(self.CurrentTopology) in [SNL]:
-            if self.WriteGridueSNL(self.CurrentTopology.gridue_settings, fname):
-                print(f"# Successfully saved gridue file as '{fname}'")
-        elif type(self.CurrentTopology) in [SF15, SF45, SF75, SF105, SF135, SF165, UDN]:
-            if self.WriteGridueDNL(self.CurrentTopology.gridue_settings, fname):
-                print(f"# Successfully saved gridue file as '{fname}'")
+        self.CurrentTopology.PrepGridue(**kwargs)
+        self.CurrentTopology.WriteGridue(fname,**kwargs)
+        print(f"# Successfully saved gridue file as '{fname}'")
 
     @staticmethod
     def ImportGridue(self, fname: str = 'gridue') -> dict:
diff --git a/build/lib/INGRID/line_tracing.py b/build/lib/INGRID/line_tracing.py
index 79a9b6f..9846e9e 100644
--- a/build/lib/INGRID/line_tracing.py
+++ b/build/lib/INGRID/line_tracing.py
@@ -847,7 +847,6 @@ def fend_Z(y):
 
         if Verbose:
             print('Drew for {} seconds\n'.format(end - start))
-        print('')
         return Line(line)
 
     def PsiCostFunc(self, xy):
diff --git a/build/lib/INGRID/topologies/snl.py b/build/lib/INGRID/topologies/snl.py
index 3067e61..27cf268 100644
--- a/build/lib/INGRID/topologies/snl.py
+++ b/build/lib/INGRID/topologies/snl.py
@@ -340,7 +340,7 @@ def OrderPatches(self):
         patches = ['A2', 'B2', 'C2', 'D2', 'E2', 'F2', 'A1', 'F1', 'B1', 'C1', 'D1', 'E1']
         self.patches = OrderedDict([(pname, self.patches[pname]) for pname in patches])
 
-    def set_gridue(self):
+    def set_gridue(self,**kwargs):
         """
         Prepares a ``gridue_settings`` dictionary with required data
         for writing a gridue file.
@@ -394,8 +394,71 @@ def set_gridue(self):
                     b[i][j][k] = _b
 
         self.gridue_settings = {
-            'nxm': nxm, 'nym': nym, 'ixpt1': ixpt1, 'ixpt2': ixpt2, 'iyseptrx1': iyseparatrix1,
+            'nxm': nxm, 'nym': nym, 'ixpt1': ixpt1, 'ixpt2': ixpt2, 'iyseptrx1': iyseparatrix1,'simagxs':self.OMFIT_psi['SIMAG'] ,'sibdrys':self.OMFIT_psi['SIBRY'],
             'rm': self.rm, 'zm': self.zm, 'psi': psi, 'br': br, 'bz': bz, 'bpol': bpol, 'bphi': bphi, 'b': b
         }
+        
+    
+    
+    def WriteGridue(self, fname: str = 'gridue',NewFormat=False,**kwargs) -> bool:
+        """
+        Write a gridue file for a single-null configuration.
+
+        Parameters
+        ----------
+        gridue_settings : dict
+            A dictionary containing grid data to be written to the gridue file.
+
+        fname : str, optional
+            The file name/path to save the gridue file to.
+            Defaults to 'gridue'.
+
+        Returns
+        -------
+            True if file was written with no errors
+        """
+
+        def format_header(gridue):
+            iheader_items = ['nxm', 'nym', 'ixpt1', 'ixpt2', 'iyseptrx1']
+            rheader_items = ['simagxs','sibdrys']
+            
+            Str=''.join(['{}'.format(gridue[item]).rjust(4) for item in iheader_items])
+            if NewFormat:
+                Str=Str+''.join(['{:16.10f}'.format(gridue[item]).rjust(17) for item in rheader_items])
+            return Str+ '\n'
+
+        def format_body(data):
+
+            delim_val = 0
+            delim_char = ''
+            body = ''
+
+            for n in range(5):
+                for j in range(len(data[0])):
+                    for i in range(len(data)):
+                        delim_val += 1
+                        val = np.format_float_scientific(data[i][j][n], precision=15, unique=False).rjust(23).replace('e', 'D')
+                        if delim_val == 3:
+                            delim_val = 0
+                            delim_char = '\n'
+                        body += val + delim_char
+                        delim_char = ''
+
+            if delim_val % 3 != 0:
+                body += '\n'
+
+            return body
+
+        f = open(fname, mode='w')
+        f.write(format_header(self.gridue_settings) + '\n')
+
+        body_items = ['rm', 'zm', 'psi', 'br', 'bz', 'bpol', 'bphi', 'b']
+        for item in body_items:
+            f.write(format_body(self.gridue_settings[item]) + '\n')
+
+        runidg = 'iogridue'
+        f.write(runidg + '\n')
+
+        f.close()
 
-        return self.gridue_settings
+        return True
diff --git a/build/lib/INGRID/topologies/udn.py b/build/lib/INGRID/topologies/udn.py
index 6f07beb..cfc2027 100644
--- a/build/lib/INGRID/topologies/udn.py
+++ b/build/lib/INGRID/topologies/udn.py
@@ -618,10 +618,78 @@ def set_gridue(self) -> dict:
                     b[i][j][k] = _b
 
         self.gridue_settings = {
-            'nxm': nxm, 'nym': nym, 'iyseparatrix1': iyseparatrix1, 'iyseparatrix2': iyseparatrix2,
+            'nxm': nxm, 'nym': nym, 'iyseparatrix1': iyseparatrix1, 'iyseparatrix2': iyseparatrix2,'simagxs':self.OMFIT_psi['SIMAG'] ,'sibdrys':self.OMFIT_psi['SIBRY'],
             'ix_plate1': ix_plate1, 'ix_cut1': ix_cut1, 'ix_cut2': ix_cut2, 'ix_plate2': ix_plate2, 'iyseparatrix3': iyseparatrix3,
             'iyseparatrix4': iyseparatrix4, 'ix_plate3': ix_plate3, 'ix_cut3': ix_cut3, 'ix_cut4': ix_cut4, 'ix_plate4': ix_plate4,
             'rm': self.rm, 'zm': self.zm, 'psi': psi, 'br': br, 'bz': bz, 'bpol': bpol, 'bphi': bphi, 'b': b, '_FILLER_': -1
         }
+    
+    
+    def WriteGridue(self, gridue_settings: dict, fname: str = 'gridue',**kwargs) -> bool:
+        """
+        Write a gridue file for a double-null configuration.
+
+        Parameters
+        ----------
+        gridue_settings : dict
+            A dictionary containing grid data to be written to the gridue file.
+
+        fname : str, optional
+            The file name/path to save the gridue file to.
+            Defaults to 'gridue'.
+
+        Returns
+        -------
+            True if file was written with no errors
+        """
 
-        return self.gridue_settings
+        def format_header(gridue):
+            header_rows = [
+                ['nxm', 'nym'],
+                ['iyseparatrix1', 'iyseparatrix2'],
+                ['ix_plate1', 'ix_cut1', '_FILLER_', 'ix_cut2', 'ix_plate2'],
+                ['iyseparatrix3', 'iyseparatrix4'],
+                ['ix_plate3', 'ix_cut3', '_FILLER_', 'ix_cut4', 'ix_plate4']
+            ]
+
+            header = ''
+            for header_items in header_rows:
+                for item in header_items:
+                    header += '{}'.format(gridue[item]).rjust(4)
+                header += '\n'
+            return header
+
+        def format_body(data):
+
+            delim_val = 0
+            delim_char = ''
+            body = ''
+
+            for n in range(5):
+                for j in range(len(data[0])):
+                    for i in range(len(data)):
+                        delim_val += 1
+                        val = np.format_float_scientific(data[i][j][n], precision=15, unique=False).rjust(23).replace('e', 'D')
+                        if delim_val == 3:
+                            delim_val = 0
+                            delim_char = '\n'
+                        body += val + delim_char
+                        delim_char = ''
+
+            if delim_val % 3 != 0:
+                body += '\n'
+
+            return body
+
+        f = open(fname, mode='w')
+        f.write(format_header(gridue_settings) + '\n')
+
+        body_items = ['rm', 'zm', 'psi', 'br', 'bz', 'bpol', 'bphi', 'b']
+        for item in body_items:
+            f.write(format_body(gridue_settings[item]) + '\n')
+
+        runidg = 'iogridue'
+        f.write(runidg + '\n')
+
+        f.close()
+        return True
diff --git a/build/lib/INGRID/utils.py b/build/lib/INGRID/utils.py
index 569119a..91367a7 100644
--- a/build/lib/INGRID/utils.py
+++ b/build/lib/INGRID/utils.py
@@ -15,9 +15,10 @@
 except:
     pass
 import matplotlib.pyplot as plt
-import pathlib
 import inspect
-from scipy.optimize import root, minimize
+from scipy.optimize import root
+
+
 
 import yaml as yml
 import os
@@ -28,8 +29,9 @@
 from omfit_eqdsk import OMFITgeqdsk
 from INGRID.interpol import EfitData
 from INGRID.line_tracing import LineTracing
-from INGRID.geometry import Point, Line, Patch, segment_intersect, orientation_between
-
+from INGRID.DistribFunc import DistribFunc
+from INGRID.geometry import Point, Line, Patch, orientation_between
+import types
 
 class IngridUtils():
     """
@@ -1070,137 +1072,9 @@ def ClassifyTopology(self, visual=False) -> None:
 
         self.config = self.LineTracer.config
 
-    def WriteGridueSNL(self, gridue_settings: dict, fname: str = 'gridue') -> bool:
-        """
-        Write a gridue file for a single-null configuration.
-
-        Parameters
-        ----------
-        gridue_settings : dict
-            A dictionary containing grid data to be written to the gridue file.
-
-        fname : str, optional
-            The file name/path to save the gridue file to.
-            Defaults to 'gridue'.
-
-        Returns
-        -------
-            True if file was written with no errors
-        """
-
-        def format_header(gridue):
-            header_items = ['nxm', 'nym', 'ixpt1', 'ixpt2', 'iyseptrx1']
-            header = ''
-            for item in header_items:
-                header += '{}'.format(gridue[item]).rjust(4)
-
-            header += '\n'
-            return header
-
-        def format_body(data):
-
-            delim_val = 0
-            delim_char = ''
-            body = ''
-
-            for n in range(5):
-                for j in range(len(data[0])):
-                    for i in range(len(data)):
-                        delim_val += 1
-                        val = np.format_float_scientific(data[i][j][n], precision=15, unique=False).rjust(23).replace('e', 'D')
-                        if delim_val == 3:
-                            delim_val = 0
-                            delim_char = '\n'
-                        body += val + delim_char
-                        delim_char = ''
+    
 
-            if delim_val % 3 != 0:
-                body += '\n'
-
-            return body
-
-        f = open(fname, mode='w')
-        f.write(format_header(gridue_settings) + '\n')
-
-        body_items = ['rm', 'zm', 'psi', 'br', 'bz', 'bpol', 'bphi', 'b']
-        for item in body_items:
-            f.write(format_body(gridue_settings[item]) + '\n')
-
-        runidg = 'iogridue'
-        f.write(runidg + '\n')
-
-        f.close()
-
-        return True
-
-    def WriteGridueDNL(self, gridue_settings: dict, fname: str = 'gridue') -> bool:
-        """
-        Write a gridue file for a double-null configuration.
-
-        Parameters
-        ----------
-        gridue_settings : dict
-            A dictionary containing grid data to be written to the gridue file.
-
-        fname : str, optional
-            The file name/path to save the gridue file to.
-            Defaults to 'gridue'.
-
-        Returns
-        -------
-            True if file was written with no errors
-        """
-
-        def format_header(gridue):
-            header_rows = [
-                ['nxm', 'nym'],
-                ['iyseparatrix1', 'iyseparatrix2'],
-                ['ix_plate1', 'ix_cut1', '_FILLER_', 'ix_cut2', 'ix_plate2'],
-                ['iyseparatrix3', 'iyseparatrix4'],
-                ['ix_plate3', 'ix_cut3', '_FILLER_', 'ix_cut4', 'ix_plate4']
-            ]
-
-            header = ''
-            for header_items in header_rows:
-                for item in header_items:
-                    header += '{}'.format(gridue[item]).rjust(4)
-                header += '\n'
-            return header
-
-        def format_body(data):
-
-            delim_val = 0
-            delim_char = ''
-            body = ''
-
-            for n in range(5):
-                for j in range(len(data[0])):
-                    for i in range(len(data)):
-                        delim_val += 1
-                        val = np.format_float_scientific(data[i][j][n], precision=15, unique=False).rjust(23).replace('e', 'D')
-                        if delim_val == 3:
-                            delim_val = 0
-                            delim_char = '\n'
-                        body += val + delim_char
-                        delim_char = ''
-
-            if delim_val % 3 != 0:
-                body += '\n'
-
-            return body
-
-        f = open(fname, mode='w')
-        f.write(format_header(gridue_settings) + '\n')
-
-        body_items = ['rm', 'zm', 'psi', 'br', 'bz', 'bpol', 'bphi', 'b']
-        for item in body_items:
-            f.write(format_body(gridue_settings[item]) + '\n')
-
-        runidg = 'iogridue'
-        f.write(runidg + '\n')
-
-        f.close()
-        return True
+    
 
     def ReconstructPatches(self, fname: str) -> tuple:
         """
@@ -1227,7 +1101,7 @@ def ReconstructPatches(self, fname: str) -> tuple:
         patches = {}
 
         for raw_patch in data[2:]:
-            patch_data, cell_data, patch_settings = raw_patch
+            patch_data, cell_data, patch_settings,Vertices,cell_grid = raw_patch
             NR = patch_data[0]
             NZ = patch_data[1]
             ER = patch_data[2]
@@ -1246,7 +1120,8 @@ def ReconstructPatches(self, fname: str) -> tuple:
                           plate_patch=patch_settings['plate_patch'],
                           plate_location=patch_settings['plate_location'],
                           PatchTagMap=patch_settings['PatchTagMap'])
-
+            patch.SetVertices(Vertices)
+            patch.cell_grid=cell_grid
             patches[patch.patch_name] = patch
 
         patches = OrderedDict([(k, v) for k, v in patches.items()])
@@ -1267,16 +1142,16 @@ def CheckPatches(self, verbose: bool = False) -> None:
         """
         self.CurrentTopology.CheckPatches(verbose=verbose)
 
-    def PrepGridue(self, guard_cell_eps=1e-3) -> None:
-        """
-        Prepare the gridue for writing.
+    # def PrepGridue(self, guard_cell_eps=1e-3) -> None:
+    #     """
+    #     Prepare the gridue for writing.
 
-        This method calls topology specific implementations of methods that
-        concatenate the Patch object subgrids into a global grid.
-        """
-        self.CurrentTopology.SetupPatchMatrix()
-        self.CurrentTopology.concat_grid(guard_cell_eps=guard_cell_eps)
-        self.CurrentTopology.set_gridue()
+    #     This method calls topology specific implementations of methods that
+    #     concatenate the Patch object subgrids into a global grid.
+    #     """
+    #     self.CurrentTopology.SetupPatchMatrix()
+    #     self.CurrentTopology.concat_grid(guard_cell_eps=guard_cell_eps)
+    #     self.CurrentTopology.set_gridue()
 
     @classmethod
     def _CheckOverlapCells(Grid, Verbose=False):
@@ -1359,6 +1234,7 @@ def __init__(self, Ingrid_obj: object, config: str):
         self.PlateData = Ingrid_obj.PlateData
         self.PatchTagMap = self.parent.GetPatchTagMap(config=config)
         self.LineTracer = Ingrid_obj.LineTracer
+        self.OMFIT_psi = Ingrid_obj.OMFIT_psi
         self.PsiUNorm = Ingrid_obj.PsiUNorm
         self.PsiNorm = Ingrid_obj.PsiNorm
         self.CurrentListPatch = {}
@@ -1404,6 +1280,7 @@ def patch_diagram(self, fig: object = None, ax: object = None) -> None:
         for i, patch in enumerate(self.patches.values()):
             patch.plot_border(color='black', ax=a)
             patch.fill(colors[patch.get_tag()[0]], ax=a, alpha=alpha[patch.get_tag()[-1]])
+            patch.alpha=alpha[patch.get_tag()[-1]]
             patch.color = colors[patch.get_tag()[0]]
         handles, labels = a.get_legend_handles_labels()
         lookup = {label: handle for label, handle in zip(labels, handles)}
@@ -1412,7 +1289,7 @@ def patch_diagram(self, fig: object = None, ax: object = None) -> None:
                  ncol=1)
         f.show()
 
-    def grid_diagram(self, fig: object = None, ax: object = None) -> None:
+    def grid_diagram(self, fig: object = None, ax: object = None,color=None,**kwargs) -> None:
         """
         Generates the grid diagram for a given configuration.
 
@@ -1441,7 +1318,8 @@ def grid_diagram(self, fig: object = None, ax: object = None) -> None:
         ax.set_title(f'{self.config} Subgrid')
 
         for patch in self.patches.values():
-            patch.plot_subgrid(ax=ax)
+            #patch.plot_subgrid(ax=ax)
+            patch.PlotSubGrid(ax=ax,color=color)
             print(f'# Plotting subgrid {patch.patch_name}')
 
         fig.show()
@@ -1480,96 +1358,37 @@ def get_config(self) -> str:
         """
         return self.config
 
-    def get_func(self, func: str) -> 'function':
-        """
-        Create a function from a string input.
-
-        Will be used to generate a poloidal or radial transformation
-        function.
-
-        Parameters
-        ----------
-        func : str
-            An expression to generate a function from.
-
-        Returns
-        -------
-            A function generated from the str input.
-
-        Examples
-        --------
-        When calling method ``get_func`` the user must provide a **string**
-        with the following general format:
-
-        .. math::
-
-            x, f(x)
-
-        That is, the dependent variable and expression to evaluate are
-        separated with a comma character.
-
-        The `Sympy` library supports most expressions that can be generated with
-        Python natively. See the `Sympy` documentation for advanced features.
-
-        Defining a function representing f(x) = x ^ 2:
-
-        >>> func = 'x, x ** 2'
-        >>> f = MyTopologyUtils.get_func(func)
-        >>> f(np.array([0, 1, 2, 3]))
-        array([0, 1, 4, 9])
-
-        Defining a function representing f(x) = exp(x)
-
-        >>> func = 'x, exp(x)'
-        >>> f = MyTopologyUtils.get_func(func)
-        >>> f(np.array([0, 1, 2, 3]))
-        array([ 1.        ,  2.71828183,  7.3890561 , 20.08553692])
-
-        """
-
-        def make_sympy_func(var, expression):
-            import sympy as sp
-            _f = sp.lambdify(var, expression, 'numpy')
-            return _f
-
-        f_str_raw = func
-
-        f_str_raw = f_str_raw.replace(' ', '')
-        delim = f_str_raw.index(',')
-
-        var = f_str_raw[0: delim]
-        expression = f_str_raw[delim + 1:]
-
-        func = make_sympy_func(var, expression)
-        # TODO: Check Normalization of the function to 1
-        return func
-
-    def CheckFunction(self, expression: str, Verbose: bool = False) -> bool:
-        """
-        Check if a str is in the correct format for method ``get_func``
-
-        Parameters
-        ----------
-        expression : str
-            Expression to check.
-
-        Verbose : bool, optional
-            Print full output to terminal. Default to False
-
-        Returns
-        -------
-            True if expression is valid. False otherwise.
-        """
-        ExpValid = False
-        try:
-            com = 'lambda {}:{}'.format(expression.split(',')[0], expression.split(',')[1])
-            if Verbose:
-                print(com)
-            eval(com)
-            ExpValid = True
-        except:
-            raise ValueError('Unable to parse expression entry "{}".'.format(expression))
-        return ExpValid
+    
+    def SetPatchesOrientation(self):
+        
+        
+        # set radial orientation negative if patch is in PFR
+        for name, patch in self.patches.items():
+            if  patch.plate_patch and self.GetPatchLetterNumber(patch.get_tag())[1]=='1':
+                patch.RadOrientation=-1
+            else:
+                patch.RadOrientation=1
+        
+        #set poloidal orientation negative if patches in high field side region (must be checked for other config)
+        Letters=list(dict.fromkeys([L for L,N in  self.GetPatchLetterNumber(self.GetListAllPatches())]))
+        Nhalf=int(len(Letters)/2)
+        ListLFS=Letters[0:Nhalf]
+        for name, patch in self.patches.items():
+            if self.GetPatchLetterNumber(patch.get_tag())[0] in ListLFS:
+                patch.PolOrientation=-1
+            else:
+                patch.PolOrientation=1
+        
+        
+        
+        
+                
+        
+    
+    
+    def ResetHeavyAttr(self):
+        for p in self.patches.values():
+            p.ResetSplines()
 
     def GetFunctions(self, Patch: Patch, Verbose: bool = False) -> tuple:
         """
@@ -1592,49 +1411,43 @@ def GetFunctions(self, Patch: Patch, Verbose: bool = False) -> tuple:
         -------
             2-tuple containing functions for radial and poloidal direction respectively.
         """
-
+        
         poloidal_tag, radial_tag = Patch.get_tag()
         p_f = 'poloidal_f_' + poloidal_tag
         r_f = 'radial_f_' + radial_tag
-
-        try:
-            _poloidal_f = self.settings['grid_settings']['grid_generation'][p_f]
-            valid_function = self.CheckFunction(_poloidal_f, Verbose)
-            if valid_function:
-                _poloidal_f = self.get_func(_poloidal_f)
-            else:
-                raise ValueError('# Invalid function entry. Applying default poloidal function.')
-        except:
-            _poloidal_f = self.settings['grid_settings']['grid_generation']['poloidal_f_default']
-            valid_function = self.CheckFunction(_poloidal_f, Verbose)
-            if valid_function:
-                _poloidal_f = self.get_func(_poloidal_f)
-            else:
-                _poloidal_f = lambda x: x
-
-        try:
-
-            # Adding CORE radial_f support for SNL cases via entry 'radial_f_3'
-            if self.config in ['USN', 'LSN'] \
-                and self.settings['grid_settings']['grid_generation'].get('radial_f_3') is not None \
-                    and poloidal_tag + radial_tag in ['B1', 'C1', 'D1', 'E1']:
-                _radial_f = self.settings['grid_settings']['grid_generation']['radial_f_3']
-            else:
-                _radial_f = self.settings['grid_settings']['grid_generation'][r_f]
-            valid_function = self.CheckFunction(_radial_f, Verbose)
-            if valid_function:
-                _radial_f = self.get_func(_radial_f)
-            else:
-                raise ValueError('# Invalid function entry. Applying default radial function.')
-        except:
-            _radial_f = self.settings['grid_settings']['grid_generation']['radial_f_default']
-            valid_function = self.CheckFunction(_radial_f, Verbose)
-            if valid_function:
-                _radial_f = self.get_func(_radial_f)
-            else:
-                _radial_f = lambda x: x
-
-        return (_radial_f, _poloidal_f)
+        
+        _poloidal_f = DistribFunc(self.settings['grid_settings']['grid_generation'].get(p_f))
+        if _poloidal_f is None:
+            _poloidal_f = DistribFunc(self.settings['grid_settings']['grid_generation'].get('poloidal_f_default'))
+        if _poloidal_f is None:
+            _poloidal_f = DistribFunc('x:x')
+
+    
+        (L,N)=self.GetPatchLetterNumber(Patch.get_tag())
+        # Adding CORE radial_f support for SNL cases via entry 'radial_f_3'
+        if self.config in ['USN', 'LSN'] and N=='1' and Patch.plate_patch:
+            _radial_f = DistribFunc(self.settings['grid_settings']['grid_generation'].get('radial_f_3'))
+        else:
+            _radial_f = DistribFunc(self.settings['grid_settings']['grid_generation'].get(r_f))
+        if _radial_f is None:
+            _radial_f = DistribFunc(self.settings['grid_settings']['grid_generation'].get('radial_f_default'))
+        if _radial_f is None:
+            _radial_f = lambda x: x
+            
+        if hasattr(Patch,'RadOrientation') and Patch.RadOrientation==-1:
+            _radial_f_=lambda x:1-_radial_f(1-x)
+        else:
+            _radial_f_=_radial_f
+        if hasattr(Patch,'PolOrientation') and Patch.PolOrientation==-1:
+            _poloidal_f_=lambda x:1-_poloidal_f(1-x)
+        else:
+            _poloidal_f_=_poloidal_f
+            
+            
+        
+        return (_radial_f_, _poloidal_f_)
+    
+    
 
     def GetNpoints(self, Patch: Patch) -> tuple:
         """
@@ -1680,15 +1493,32 @@ def GetDistortionCorrectionSettings(self) -> dict:
             A dictionary containing ``CorrectDistortion`` settings.
         """
         if self.settings['grid_settings']['grid_generation'].get('distortion_correction') is not None:
-            CD = self.settings['grid_settings']['grid_generation']['distortion_correction']
+            self.distortion_correction = self.settings['grid_settings']['grid_generation']['distortion_correction']
         else:
-            CD = {}
-
-        self.distortion_correction = CD
-        return CD
-
+            self.distortion_correction = {}
+
+        if hasattr(self,'patches'):
+            for name, patch in self.patches.items():
+                (L,N)=self.GetPatchLetterNumber(patch.get_tag())
+                if self.distortion_correction.get(name) is not None:
+                    patch.distortion_correction = self.distortion_correction.get(name)
+                elif self.distortion_correction.get(patch.get_tag()) is not None:
+                    patch.distortion_correction = self.distortion_correction.get(patch.get_tag())
+                elif self.distortion_correction.get(N) is not None:
+                    patch.distortion_correction = self.distortion_correction.get(N)
+                elif self.distortion_correction.get(L) is not None:
+                    patch.distortion_correction = self.distortion_correction.get(L)    
+                elif self.distortion_correction.get('all') is not None:
+                    patch.distortion_correction = self.distortion_correction.get('all')
+                else:
+                    patch.distortion_correction = {'active': False}
+        
+
+    def GetListAllPatches(self):
+        return [v.get_tag() for k,v in self.patches.items()]
+    
     def construct_grid(self, np_cells: int = 1, nr_cells: int = 1, Verbose: bool = False,
-                       ShowVertices: bool = False, RestartScratch: bool = False, ListPatches: str = 'all') -> None:
+                       ShowVertices: bool = False, RestartScratch: bool = False, ListPatches: str = 'all',Visual=True) -> None:
         """
         Construct a grid by refining a Patch map.
 
@@ -1756,45 +1586,100 @@ def construct_grid(self, np_cells: int = 1, nr_cells: int = 1, Verbose: bool = F
             visual = self.settings['DEBUG']['visual']['subgrid']
         except:
             visual = False
+        if Visual:
+            visual=True
         try:
             verbose = self.settings['DEBUG']['verbose']['grid_generation']
         except:
             verbose = False
 
-        verbose = Verbose or verbose
-
+        if Verbose:
+            verbose=True
+        ListAllPatches=self.GetListAllPatches()
+        print('>>> Existing Patches :',ListAllPatches)
+        ListPatches=self.ProcessListPatches(ListPatches)
+        print('>>> Generating grids for patches:',ListPatches)
+        self.SetPatchesOrientation()
         self.GetDistortionCorrectionSettings()
-
-        print('>>> Patches:', [k for k in self.patches.keys()])
+        
+        
         if RestartScratch:
             self.CurrentListPatch = {}
 
         for name, patch in self.patches.items():
-
-            if self.distortion_correction.get(name) is not None:
-                patch.distortion_correction = self.distortion_correction.get(name)
-            elif self.distortion_correction.get(patch.get_tag()) is not None:
-                patch.distortion_correction = self.distortion_correction.get(patch.get_tag())
-            elif self.distortion_correction.get('all') is not None:
-                patch.distortion_correction = self.distortion_correction.get('all')
-            else:
-                patch.distortion_correction = {'Active': False}
-            if (ListPatches == 'all' and patch not in self.CurrentListPatch) or (ListPatches != 'all' and name in ListPatches):
-                self.SetPatchBoundaryPoints(patch, verbose)
+            self.SetPatchBoundaryPoints(patch)
+            
+        for name, patch in self.patches.items():
+            if (ListPatches == 'all' and patch not in self.CurrentListPatch) or (ListPatches != 'all' and patch.get_tag() in ListPatches):
+                self.SetPatchBoundaryPoints(patch)
                 (nr_cells, np_cells) = self.GetNpoints(patch)
                 (_radial_f, _poloidal_f) = self.GetFunctions(patch)
-                print(f'>>> Making subgrid in patch {name}:')
-                print(f'    np = {np_cells}, nr = {nr_cells}')
-                print(f'    fp = {inspect.getsource(_poloidal_f)}')
-                print(f'    fr = {inspect.getsource(_radial_f)}', end='')
+                print('>>> Making subgrid in patch:{} with np={},nr={},fp={},fr={}'.format(patch.get_tag(), np_cells, nr_cells, inspect.getsource(_poloidal_f), inspect.getsource(_radial_f)))
                 patch.RemoveDuplicatePoints()
                 patch.make_subgrid(self, np_cells, nr_cells, _poloidal_f=_poloidal_f, _radial_f=_radial_f, verbose=verbose, visual=visual, ShowVertices=ShowVertices)
 
+                
                 self.CurrentListPatch[name] = patch
-                print(f'    {name} subgrid complete.\n\n')
+                
         self.AdjustGrid()
-
-    def SetPatchBoundaryPoints(self, Patch: Patch, verbose: bool = False) -> None:
+    
+    @staticmethod
+    def GetPatchLetterNumber(PatchName:str or str)->tuple:
+        
+        if type(PatchName)==str:
+            if len(PatchName)<1 or len(PatchName)>2:
+                raise IOError('PatchName must be a string of one or two characters')
+            L=''
+            N=None
+            
+            for C in PatchName:
+                if C.isnumeric():
+                    N=C
+                elif C.isalpha():
+                    L=C
+                else:
+                    raise IOError('PatchName must be a string of one or two characters with one letter max and one digit max')
+            return (L,N) 
+        
+        elif type(PatchName)==list:
+            Out=[]
+            for P in PatchName:
+                Out.append(TopologyUtils.GetPatchLetterNumber(P))
+            return Out
+        else:
+            raise ValueError('PatchName must be a list of a str')
+        
+    def ProcessListPatches(self,ListPatches:list or str):
+        if type(ListPatches)==str:
+            ListPatches=[ListPatches]
+        OrderedList=self.GetListAllPatches()
+  
+        Numbers=list(dict.fromkeys([N for L,N in  self.GetPatchLetterNumber(OrderedList) if N is not None]))
+        Letters=list(dict.fromkeys([L for L,N in  self.GetPatchLetterNumber(OrderedList) if L!='']))
+        if self.Verbose:
+            print('Numbers:',Numbers)
+            print('Letters:',Letters)
+        if 'all' in ListPatches:
+            return OrderedList
+        else:
+            
+            OutList=ListPatches.copy()
+            for P in ListPatches:
+                (L,N)=self.GetPatchLetterNumber(P)
+                if L=='' and N is None:
+                    raise ValueError('Incorrect patch name pattern:{}'.format(P))
+                elif N is None:   
+                    OutList.remove(P)
+                    OutList.extend([L+Num for Num in Numbers])
+                elif L=='':
+                    OutList.remove(P)
+                    OutList.extend([Let+N for Let in Letters])
+            # now let order them 
+            OutList=[P for P in OrderedList if P in list(dict.fromkeys(OutList))]
+            return OutList
+        
+
+    def SetPatchBoundaryPoints(self, Patch: Patch) -> None:
         """
         Set the Patch ``BoundaryPoints`` dict based off TopologyUtils ``ConnexionMap``.
 
@@ -1802,8 +1687,6 @@ def SetPatchBoundaryPoints(self, Patch: Patch, verbose: bool = False) -> None:
         ----------
         Patch : Patch
             The Patch to set the boundary points for.
-        verbose: bool
-            Print full output to terminal.
 
         Notes
         -----
@@ -1812,12 +1695,12 @@ def SetPatchBoundaryPoints(self, Patch: Patch, verbose: bool = False) -> None:
         """
         Patch.TerminatesLoop = False
         if self.ConnexionMap.get(Patch.get_tag()) is not None:
-            if verbose:
+            if self.Verbose:
                 print('Find connexion map for patch {}'.format(Patch.patch_name))
             for v in self.ConnexionMap.get(Patch.get_tag()).items():
                 Boundary, AdjacentPatch = v
                 Patch.BoundaryPoints[Boundary] = self.GetBoundaryPoints(AdjacentPatch)
-                if verbose:
+                if self.Verbose:
                     print('Find Boundaries points for {}'.format(Patch.patch_name))
             if self.ConnexionMap.get(Patch.get_tag()).get('E') is not None:
                 Patch.TerminatesLoop = True
@@ -2189,3 +2072,18 @@ def set_guard(cell_map, ix, iy, eps, boundary):
 
             if debug:
                 self._animate_grid()
+                
+    def SetBoundarySplines(self):
+        for n,p in self.patches.items():
+            p.CreateBoundarySplines(self)
+    
+    def PrepGridue(self, guard_cell_eps=1e-3,**kwargs) -> None:
+        """
+        Prepare the gridue for writing.
+
+        This method calls topology specific implementations of methods that
+        concatenate the Patch object subgrids into a global grid.
+        """
+        self.SetupPatchMatrix()
+        self.concat_grid(guard_cell_eps=guard_cell_eps)
+        self.set_gridue()
diff --git a/dist/INGRID-1.0.1-py3.7.egg b/dist/INGRID-1.0.1-py3.7.egg
index ba73b54..6c8dc1e 100644
Binary files a/dist/INGRID-1.0.1-py3.7.egg and b/dist/INGRID-1.0.1-py3.7.egg differ