diff --git a/Scripts_Matthieu/AnalysisFunctions.py b/Scripts_Matthieu/AnalysisFunctions.py
index b1a25f030c2c065e33aadc37f6ea1a31b7725c3e..03a2706a1089bc4b8866fdf383a8dd745f7b000d 100644
--- a/Scripts_Matthieu/AnalysisFunctions.py
+++ b/Scripts_Matthieu/AnalysisFunctions.py
@@ -29,8 +29,8 @@ def extractValues(positions, drifts, efield, gradb, drgradb, edescr, bdescr, gra
         dbn = [float(x) for x in drgradb[i].split()]
 
         # filter our boundary cells => These have zero field and drift and are not useful for interpretation!
-        if dn != [0,0,0] and not math.isnan(dn[0]):
-            denum.append(dn)
+        if den != [0,0,0] and not math.isnan(den[0]):
+            denum.append(den)
             pos.append(dpn)
             enum.append(en)
             gnum.append(gbn)
@@ -39,7 +39,7 @@ def extractValues(positions, drifts, efield, gradb, drgradb, edescr, bdescr, gra
             efa = edescr(dpn[0], dpn[2], dpn[1], Constants.Rmajor)
             bfa = bdescr(dpn[0], dpn[2], dpn[1], Constants.Rmajor, Constants.Rminor)
             gfa = gradbdescr(dpn[0], dpn[2], dpn[1], Constants.Rmajor, Constants.Rminor)
-            dga.append(gradBdrift(bga,gfa).tolist())
+            dga.append(gradBdrift(bfa,gfa).tolist())
             ga.append(gfa)
             ea.append(efa)
             ba.append(bfa)
@@ -93,8 +93,8 @@ def getAllSlices(positions,allInfo,direction,coordinate):
     return result
 
 def getAllInfo(positions,drifts,efield, gradb, drgradb,edescr,bdescr, gradbdescr, ExB, gradBdrift):
-    [pos, enum, denum,gnum,dgnum, enumperp,denumperp, gnumperp ,dgnumperp,ea , ba, dea, ga,dga] =
-    extractValues(positions, drifts, efield, gradb, drgradb, edescr, bdescr,gradbdescr, ExB, gradBdrift)
+    [pos, enum, denum,gnum,dgnum, enumperp,denumperp, gnumperp ,dgnumperp,ea , ba, dea, ga,dga] = extractValues(positions,
+                        drifts, efield, gradb, drgradb, edescr, bdescr,gradbdescr, ExB, gradBdrift)
 
     [deres_mag, dean_mag, deer_mag, deer_vec, deer_rel_mag, deer_rel_vec] =  extractVectorInfo(denum,dea)
     [eres_mag, ean_mag, eer_mag, eer_vec, eer_rel_mag, eer_rel_vec] = extractVectorInfo(enum, ea)
@@ -106,12 +106,11 @@ def getAllInfo(positions,drifts,efield, gradb, drgradb,edescr,bdescr, gradbdescr
     torSliceDict = {}
     keys = ['edrifts', 'efield','edrifts_perp','efield_perp' ,'ef_mag', 'ef_an_mag', 'edrift_mag', 'edrift_an_mag', 'ef_vec_er', 'ef_mag_er',
             'edrift_vec_er', 'edrift_er_mag','efield_er_rel_mag', 'efield_er_rel_vec', 'edrift_er_rel_mag', 'edrift_er_rel_vec',
-	    'gdrifts', 'gradb','gdrifts_perp','gradb_perp', 'gb_mag','gb_an_mag', 'gdrift_mag','gdrift_an_mag','gb_vec_er', 'gb_mag_er',
-	     'gdrift_vec_er','gdrift_er_mag','gb_er_rel_mag','gb_er_rel_vec','gdrift_er_rel_mag', 'gb_er_rel_vec', 'pos']
-    #torslices = getAllSlices(pos,[dres_mag, dan_mag, der_mag, der_vec, der_rel_mag, der_rel_vec,eres_mag, ean_mag, eer_mag, eer_vec, eer_rel_mag, eer_rel_vec],1,torsurf)
+            'gdrifts', 'gradb','gdrifts_perp','gradb_perp', 'gb_mag','gb_an_mag', 'gdrift_mag','gdrift_an_mag','gb_vec_er', 'gb_mag_er',
+            'gdrift_vec_er','gdrift_er_mag','gb_er_rel_mag','gb_er_rel_vec','gdrift_er_rel_mag', 'gb_er_rel_vec', 'pos']
 
-    values = [denum,enum,denumperp,enumperp,eres_mag,ean_mag,deres_mag,dean_mag,eer_vec,eer_mag,deer_vec,der_mag,eer_rel_mag,eer_rel_vec,deer_rel_mag,deer_rel_vec,
-             dgnum,enum,dgnumperp,enumperp,eres_mag,ean_mag,dgres_mag,dgan_mag,eer_vec,eer_mag,dger_vec,dgr_mag,eer_rel_mag,eer_rel_vec,dger_rel_mag,dger_rel_vec,pos]
+    values = [denum,enum,denumperp,enumperp,eres_mag,ean_mag,deres_mag,dean_mag,eer_vec,eer_mag,deer_vec,deer_mag,eer_rel_mag,eer_rel_vec,deer_rel_mag,deer_rel_vec,
+              dgnum,enum,dgnumperp,enumperp,eres_mag,ean_mag,dgres_mag,dgan_mag,eer_vec,eer_mag,dger_vec,dger_mag,eer_rel_mag,eer_rel_vec,dger_rel_mag,dger_rel_vec,pos]
     for i in range(len(values)):
         torSliceDict[keys[i]] = values[i]
 
diff --git a/Scripts_Matthieu/AnalyticalPrescriptions.py b/Scripts_Matthieu/AnalyticalPrescriptions.py
index 36233eca8f8eaf27f699ec659b226519e0920c18..7bc8e37ec674fdd0e0e874063570dfde4abd4b1e 100644
--- a/Scripts_Matthieu/AnalyticalPrescriptions.py
+++ b/Scripts_Matthieu/AnalyticalPrescriptions.py
@@ -36,7 +36,7 @@ def potential_description2(r,z,theta,Rmajor,Rminor):
 def potential_description3(r,z,theta,Rmajor,Rminor):
     r = (r-Rmajor)/100
     z = z/100
-    potval = 500*math.cos(4*math.atan2(z,r))
+    potval = 500*math.cos(math.atan2(z,r))
     return potval
 # Electric Fields
 def e_descr1(r,z,theta, Rmajor):
diff --git a/Scripts_Matthieu/AnalyzeNumericalResults.py b/Scripts_Matthieu/AnalyzeNumericalResults.py
index 1742b0c79afc096a92288601e9ad0560d797c62a..adf9695ccf61b6119532b7988e444270af15eed2 100644
--- a/Scripts_Matthieu/AnalyzeNumericalResults.py
+++ b/Scripts_Matthieu/AnalyzeNumericalResults.py
@@ -51,22 +51,23 @@ for line in gradb_file:
 
 
 # Take slices in different directions to allow for extra analysis
-#print(drift_pos_num)
-torsurf = drift_pos_num[0][1]
-radsurf = drift_pos_num[150][0] # Based on inspection
+
+torsurf = pos[0][1]
+radsurf = pos[150][0] # Based on inspection
 
 # Create a dictionary containing the different slices
 torSliceDict = {}
 keys = ['edrifts', 'efield','edrifts_perp','efield_perp' ,'ef_mag', 'ef_an_mag', 'edrift_mag', 'edrift_an_mag', 'ef_vec_er', 'ef_mag_er',
             'edrift_vec_er', 'edrift_er_mag','efield_er_rel_mag', 'efield_er_rel_vec', 'edrift_er_rel_mag', 'edrift_er_rel_vec',
-	    'gdrifts', 'gradb','gdrifts_perp','gradb_perp', 'gb_mag','gb_an_mag', 'gdrift_mag','gdrift_an_mag','gb_vec_er', 'gb_mag_er',
-	     'gdrift_vec_er','gdrift_er_mag','gb_er_rel_mag','gb_er_rel_vec','gdrift_er_rel_mag', 'gb_er_rel_vec', 'pos']
+        'gdrifts', 'gradb','gdrifts_perp','gradb_perp', 'gb_mag','gb_an_mag', 'gdrift_mag','gdrift_an_mag','gb_vec_er', 'gb_mag_er',
+        'gdrift_vec_er','gdrift_er_mag','gb_er_rel_mag','gb_er_rel_vec','gdrift_er_rel_mag', 'gb_er_rel_vec', 'pos']
+
 allInfo = [enum, denum,gnum,dgnum, enumperp,denumperp, gnumperp ,dgnumperp,ea , ba, dea, ga,dga,
            edrift_res_mag,edrift_an_mag,edrift_er_mag,edrift_er_vec,edrift_er_rel_mag,edrift_er_rel_vec,
-	   efield_res_mag,efield_an_mag,efield_er_mag,efield_er_vec,efield_er_rel_mag,efield_er_rel_vec,
-	   gdrift_res_mag,gdrift_an_mag,gdrift_er_mag,gdrift_er_vec,gdrift_er_rel_mag,gdrift_er_rel_vec,
-	   gradb_res_mag,gradb_an_mag,gradb_er_mag,gradb_er_vec,gradb_er_rel_mag,gradb_er_rel_vec]
-					    	     
+           efield_res_mag,efield_an_mag,efield_er_mag,efield_er_vec,efield_er_rel_mag,efield_er_rel_vec,
+           gdrift_res_mag,gdrift_an_mag,gdrift_er_mag,gdrift_er_vec,gdrift_er_rel_mag,gdrift_er_rel_vec,
+           gradb_res_mag,gradb_an_mag,gradb_er_mag,gradb_er_vec,gradb_er_rel_mag,gradb_er_rel_vec]
+
 torslices = AnalysisFunctions.getAllSlices(pos,allInfo,1,torsurf)
 for i in range(len(torslices)):
     torSliceDict[keys[i]] = torslices[i]
diff --git a/Scripts_Matthieu/PlotCase1.py b/Scripts_Matthieu/PlotCase1.py
index 1535885e6d73404034b36dd2f6c12354fee7d47a..485aca539bca911ea9c194837bb73fea3bdb5e14 100644
--- a/Scripts_Matthieu/PlotCase1.py
+++ b/Scripts_Matthieu/PlotCase1.py
@@ -4,7 +4,7 @@ import matplotlib
 import matplotlib.pyplot as plt
 from MeshFunctions import *
 import math
-
+import AnalyticalPrescriptions
 from driftFunctions import *
 
 
@@ -109,4 +109,4 @@ plt.show()
 # Plot 2d contour of Bfield
 ns = 100
 plot2DContour(ns,ns,Rmajor-rminor,Rmajor+rminor,-rminor,rminor,Rmajor,rminor,b_field_description,'X [cm]', 'Y [cm]','Magnetic field strength [T]',True)
-plot2DContour(ns,ns,Rmajor-rminor,Rmajor+rminor,-rminor,rminor,Rmajor,rminor,potential,'X [cm]', 'Y [cm]','Potential [V]',True)
\ No newline at end of file
+plot2DContour(ns,ns,Rmajor-rminor,Rmajor+rminor,-rminor,rminor,Rmajor,rminor,AnalyticalPrescriptions.potential_description3,'X [cm]', 'Y [cm]','Potential [V]',True)
\ No newline at end of file
diff --git a/Scripts_Matthieu/PrepareGeometryDescription1.py b/Scripts_Matthieu/PrepareGeometryDescription1.py
index f2af80e2bded541d817b26c0d4b5660a3613c58a..6e8a83219059ec7da992c01a667aaf15135f2be4 100644
--- a/Scripts_Matthieu/PrepareGeometryDescription1.py
+++ b/Scripts_Matthieu/PrepareGeometryDescription1.py
@@ -4,7 +4,7 @@ import shutil
 import AnalyticalPrescriptions
 # The folder where the different cases should be stored is
 basepath = "D:\Matthieu\Documents\\1.Master\Thesis\Data\\"
-folder = "D:\Matthieu\Documents\\1.Master\Thesis\Data\Description2"
+folder = "D:\Matthieu\Documents\\1.Master\Thesis\Data\TorSurf"
 if os.path.isdir(folder):
     Flag = int(input("Do you want to remove" + folder + " Answer should be 1/0"))
     if Flag:
@@ -17,14 +17,14 @@ os.mkdir(folder)
 # Number of zones considered
 nz = 1
 # Loop over different cases to prepare
-nt = 10
-nrs = [10,50,100,200,400,1000]
-nps = nrs
-cases = [nrs,nps]
-caseNames = ['Case10', 'Case50', 'Case100','Case200', 'Case400', 'Case1000']
+nts = [5,10,20,40,80,160]
+nr = 200
+np= nr
+#cases = [nrs,nps]
+caseNames = ['Case5','Case10', 'Case20', 'Case40','Case80', 'Case160']
 # Number of toroidal, radial, poloidal surfaces
 
-for i in range(len(cases[0])):
+for i in range(len(nts)):
     # A specific folder for the case should be defined
     path = folder + "\\" + caseNames[i]
     if os.path.isdir(path):
@@ -36,8 +36,7 @@ for i in range(len(cases[0])):
     os.mkdir(path)
     path = path + "\\Inputs"
     os.mkdir(path)
-    nr = cases[0][i]
-    np = cases[1][i]
+    nt = nts[i]
     potentialdescription = AnalyticalPrescriptions.potential_description1
     bfielddescription = AnalyticalPrescriptions.bfieldstrength1
     # toroidal range (degrees)