From 1cb86f9cc62923e1ae1addf0106df6592f4447af Mon Sep 17 00:00:00 2001
From: Matthieu Jacobs <matthieu.jacobs@student.kuleuven.be>
Date: Thu, 7 Apr 2022 09:41:28 +0200
Subject: [PATCH] update python files

---
 Scripts_Matthieu/Analysis2D.py                | 209 ++++++++++++++++++
 Scripts_Matthieu/AnalysisFunctions.py         |  98 +++++---
 Scripts_Matthieu/AnalyticalPrescriptions.py   |  54 +++--
 Scripts_Matthieu/AnalyzeNumericalResults.py   |  31 +--
 Scripts_Matthieu/DoAnalysis.sh                |  46 ++++
 Scripts_Matthieu/MeshFunctions.py             |   9 +-
 Scripts_Matthieu/PrepareGeometry.py           |  37 ++--
 .../PrepareGeometryCompareNodal.py            |  51 +++++
 .../PrepareGeometryDescription1.py            |  31 +--
 .../PrepareGeometryOptimalResolution.py       |  51 +++++
 Scripts_Matthieu/ResolutionStudyGeometries.py |  82 +++++++
 Scripts_Matthieu/driftFunctions.py            |   6 +-
 Scripts_Matthieu/test.py                      |  31 +++
 13 files changed, 625 insertions(+), 111 deletions(-)
 create mode 100644 Scripts_Matthieu/Analysis2D.py
 create mode 100755 Scripts_Matthieu/DoAnalysis.sh
 create mode 100644 Scripts_Matthieu/PrepareGeometryCompareNodal.py
 create mode 100644 Scripts_Matthieu/PrepareGeometryOptimalResolution.py
 create mode 100644 Scripts_Matthieu/ResolutionStudyGeometries.py
 create mode 100644 Scripts_Matthieu/test.py

diff --git a/Scripts_Matthieu/Analysis2D.py b/Scripts_Matthieu/Analysis2D.py
new file mode 100644
index 0000000..ea11838
--- /dev/null
+++ b/Scripts_Matthieu/Analysis2D.py
@@ -0,0 +1,209 @@
+import math
+import driftFunctions
+import AnalyticalPrescriptions
+import Constants
+import Utilities
+import pandas as pd
+import matplotlib.pyplot as plt
+import PlottingFunctions
+import AnalysisFunctions
+import sys
+import os
+import seaborn as sns
+import numpy as np
+
+method = sys.argv[1]
+anCase = sys.argv[2]
+#Case = sys.argv[3]
+path = '' # Case+'/'
+allInfo = []
+
+info = AnalysisFunctions.readAllInfo(method,anCase,'') #Case+'/'
+keys = list(info.keys())
+#print(keys)
+for i in range(len(keys)-1): # don't include pos
+    allInfo.append(info.get(keys[i]))
+
+pos = info.get('pos')
+angles = [r[1] for r in pos]
+uniqueangles = list(set(angles))    
+print(uniqueangles)
+#print(info.get('efield_er_rel_mag'))
+# Get a slice for 2D visualization
+
+torsurf = uniqueangles[2]
+torSliceDict = {}
+torslices = AnalysisFunctions.getAllSlices(pos,allInfo,1,torsurf)
+for i in range(len(torslices)):
+    torSliceDict[keys[i]] = torslices[i]
+    
+# Make Plots
+x = [r[0] for r in torSliceDict.get('pos')]
+y = [r[1] for r in torSliceDict.get('pos')]
+
+potentials = [AnalyticalPrescriptions.potential_description1(r[0],r[2],r[1],500,100) for r in pos]
+gLengths = AnalysisFunctions.gradientLength(torSliceDict.get('ef_an_mag'),potentials)
+
+plot1 = plt.figure(1)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'Magnitude of Electic Field (Numerical)'
+plt.title(title)
+sc1 = plt.scatter(x, y, s=20, c=torSliceDict.get('ef_mag'), cmap='plasma')
+plt.colorbar(sc1)
+plt.savefig(path+method+anCase+title+'.png')
+
+plot2 = plt.figure(2)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'Horizontal component of Electric Field'
+plt.title(title)
+sc2 = plt.scatter(x, y, s=20, c = [x[0] for x in torSliceDict.get('ef_an')], cmap='plasma')
+plt.colorbar(sc2)
+plt.savefig(path+method+anCase+title+'.png')
+
+plot3 = plt.figure(3)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'Relative magnitude error of Electric Field'
+plt.title(title)
+sc3 = plt.scatter(x, y, s=20, c=torSliceDict.get('efield_er_rel_mag'), cmap='plasma')
+plt.colorbar(sc3)
+plt.savefig(path+method+anCase+title+'.png')
+
+plot4 = plt.figure(4)
+title = 'Potential'
+plt.title(title)
+PlottingFunctions.plot2DContour(Constants.ns,Constants.ns,Constants.Rmajor-Constants.Rminor,Constants.Rmajor+Constants.Rminor,-Constants.Rminor,Constants.Rminor,Constants.Rmajor,Constants.Rminor,AnalyticalPrescriptions.potential_description1,'X [cm]', 'Y [cm]','Potential [V]',False)
+plt.savefig(path+method+anCase+title+'.png')
+
+plot5 = plt.figure(5)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'Magnitude of Electric Field (Analytical)'
+plt.title(title)
+sc5 = plt.scatter(x, y, s=20, c = torSliceDict.get('ef_an_mag'), cmap='plasma')
+plt.colorbar(sc5)
+plt.savefig(path+method+anCase+title+'.png')
+
+
+plot6 = plt.figure(6)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'Magnitude of ExB Drift (Analytical)'
+plt.title(title)
+sc6 = plt.scatter(x, y, s=20, c=torSliceDict.get('edrift_an_mag'), cmap='plasma')
+plt.colorbar(sc6)
+plt.savefig(path+method+anCase+title+'.png')
+
+
+plot7 = plt.figure(7)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'Relative magnitude error of ExB Drift'
+plt.title(title)
+sc7 = plt.scatter(x, y, s=20, c=torSliceDict.get('edrift_er_rel_mag'), cmap='plasma')
+plt.colorbar(sc7)
+plt.savefig(path+method+anCase+title+'.png')
+
+
+plot8 = plt.figure(8)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'Magnitude of Efield (radial)'
+plt.title(title)
+sc8 = plt.scatter(x, y, s=20, c= [x[0] for x in torSliceDict.get('efield')], cmap='plasma')
+plt.colorbar(sc8)
+plt.savefig(path+method+anCase+title+'.png')
+
+
+plot9 = plt.figure(9)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'Magnitude of Efield (poloid)'
+plt.title(title)
+sc9 = plt.scatter(x, y, s=20, c= [x[2] for x in torSliceDict.get('efield')], cmap='plasma')
+plt.colorbar(sc9)
+plt.savefig(path+method+anCase+title+'.png')
+
+plot10 = plt.figure(10)
+sns.violinplot(torSliceDict.get('efield_er_rel_mag'))
+title = 'ErrorDistribution'+method+anCase
+plt.savefig(path+method+anCase+title+'.png')
+
+plot11 = plt.figure(11)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'GradientLengths'
+plt.title(title)
+sc11 = plt.scatter(x, y, s=20, c=gLengths, cmap='plasma')
+plt.colorbar(sc11)
+plt.savefig(path+method+anCase+title+'.png')
+
+normErrors = np.multiply(np.array(torSliceDict.get('edrift_er_rel_mag')),np.array(gLengths))
+plot12 = plt.figure(12)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'NormalizedRelErrors'
+plt.title(title)
+sc12 = plt.scatter(x, y, s=20, c=normErrors, cmap='plasma')
+plt.colorbar(sc12)
+plt.savefig(path+method+anCase+title+'.png')
+
+plot13 = plt.figure(13)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'Relative magnitude error of GradB Drift'
+plt.title(title)
+sc13 = plt.scatter(x, y, s=20, c=torSliceDict.get('gdrift_er_rel_mag'), cmap='plasma')
+plt.colorbar(sc13)
+plt.savefig(path+method+anCase+title+'.png')
+
+plot14 = plt.figure(14)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'Magnitude of GradB Drift (Analytical)'
+plt.title(title)
+sc14 = plt.scatter(x, y, s=20, c=torSliceDict.get('gdrift_an_mag'), cmap='plasma')
+plt.colorbar(sc14)
+plt.savefig(path+method+anCase+title+'.png')
+
+plot15 = plt.figure(15)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'Magnitude of GradB Drift (Numerical)'
+plt.title(title)
+sc15 = plt.scatter(x, y, s=20, c=torSliceDict.get('gdrift_mag'), cmap='plasma')
+plt.colorbar(sc15)
+plt.savefig(path+method+anCase+title+'.png')
+
+plot16 = plt.figure(16)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'Relative magnitude error of GradB '
+plt.title(title)
+sc16 = plt.scatter(x, y, s=20, c=torSliceDict.get('gb_er_rel_mag'), cmap='plasma')
+plt.colorbar(sc16)
+plt.savefig(path+method+anCase+title+'.png')
+
+plot17 = plt.figure(17)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'Magnitude of GradB (Analytical)'
+plt.title(title)
+sc17 = plt.scatter(x, y, s=20, c=torSliceDict.get('gb_an_mag'), cmap='plasma')
+plt.colorbar(sc17)
+plt.savefig(path+method+anCase+title+'.png')
+
+plot18 = plt.figure(18)
+plt.xlabel('Radial Position [cm]')
+plt.ylabel('Vertical Position [cm]')
+title = 'Magnitude of GradB (Numerical)'
+plt.title(title)
+sc18 = plt.scatter(x, y, s=20, c=torSliceDict.get('gb_mag'), cmap='plasma')
+plt.colorbar(sc18)
+plt.savefig(path+method+anCase+title+'.png')
+
+
+plt.show()
+print('Done Analysis2D')
diff --git a/Scripts_Matthieu/AnalysisFunctions.py b/Scripts_Matthieu/AnalysisFunctions.py
index 03a2706..3e3da76 100644
--- a/Scripts_Matthieu/AnalysisFunctions.py
+++ b/Scripts_Matthieu/AnalysisFunctions.py
@@ -3,8 +3,11 @@
 import Constants
 import math
 import Utilities
+import pandas as pd
+import numpy as np
+import ast
 
-def extractValues(positions, drifts, efield, gradb, drgradb, edescr, bdescr, gradbdescr, ExB, gradBdrift):
+def extractValues(positions, drifts, efield, gradb, drgradb, bcells, edescr, bdescr, gradbdescr, ExB, gradBdrift):
     pos = []
     enum = []
     denum = []
@@ -19,6 +22,7 @@ def extractValues(positions, drifts, efield, gradb, drgradb, edescr, bdescr, gra
     dea = []
     dga = []
     ga = []
+    bound = []
     
     for i in range(len(drifts)):
         # Convert to float values to allow computations and interpretation
@@ -27,9 +31,10 @@ def extractValues(positions, drifts, efield, gradb, drgradb, edescr, bdescr, gra
         en = [float(x) for x in efield[i].split()]
         gbn = [float(x) for x in gradb[i].split()]
         dbn = [float(x) for x in drgradb[i].split()]
+        bnd = float(bcells[i])
 
         # filter our boundary cells => These have zero field and drift and are not useful for interpretation!
-        if den != [0,0,0] and not math.isnan(den[0]):
+        if  bnd != 1: #not math.isnan(en[0]):# and en != [0,0,0]:
             denum.append(den)
             pos.append(dpn)
             enum.append(en)
@@ -39,15 +44,16 @@ 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(bfa,gfa).tolist())
+            dgfa = gradBdrift(bfa,gfa)
+            dga.append(dgfa.tolist())
             ga.append(gfa)
             ea.append(efa)
             ba.append(bfa)
             dea.append(ExB(efa, bfa).tolist())
-            denumperp.append(transformCyltoRadPerp(denum[-1],pos[-1]))
-            enumperp.append(transformCyltoRadPerp(enum[-1],pos[-1]))
-            gnumperp.append(transformCyltoRadPerp(gnum[-1],pos[-1]))
-            dgnumperp.append(transformCyltoRadPerp(dgnum[-1],pos[-1]))
+            denumperp.append(transformRadPerptoCyl(denum[-1],pos[-1]))
+            enumperp.append(transformRadPerptoCyl(enum[-1],pos[-1]))
+            gnumperp.append(transformRadPerptoCyl(gnum[-1],pos[-1]))
+            dgnumperp.append(transformRadPerptoCyl(dgnum[-1],pos[-1]))
     return [pos, enum, denum,gnum,dgnum, enumperp,denumperp, gnumperp ,dgnumperp,ea , ba, dea, ga,dga]
 
 def extractVectorInfo(numeric,analytical):
@@ -60,10 +66,12 @@ def extractVectorInfo(numeric,analytical):
     for i in range(len(numeric)):
         res_mag.append(math.sqrt(sum([x**2 for x in numeric[i]])))
         an_mag.append(math.sqrt(sum([x**2 for x in analytical[i]])))
-        er_mag.append(math.sqrt((res_mag[i]-an_mag[i])**2))
         er_vec.append(Utilities.listdif(numeric[i],analytical[i]))
-        er_rel_vec.append([math.sqrt(x**2)/max(1,an_mag[i]) for x in er_vec[i]])
-        er_rel_mag.append(er_mag[i]/max(1,an_mag[i]))
+        er_mag.append(math.sqrt(sum([x**2 for x in er_vec[i]])))
+        er_rel_vec.append([math.sqrt(x**2)/an_mag[i] for x in er_vec[i]])
+        if an_mag[i] == 0:
+            print('ZERO!')
+        er_rel_mag.append(sum(er_rel_vec[i]))
 
     return [res_mag,an_mag,er_mag,er_vec,er_rel_mag,er_rel_vec]
 
@@ -92,29 +100,30 @@ def getAllSlices(positions,allInfo,direction,coordinate):
     result.append(pos_slice)
     return result
 
-def getAllInfo(positions,drifts,efield, gradb, drgradb,edescr,bdescr, gradbdescr, ExB, gradBdrift):
+def getAllInfo(positions,drifts,efield, gradb, drgradb,bcells,edescr,bdescr, gradbdescr, ExB, gradBdrift,name):
     [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)
+                        drifts, efield, gradb, drgradb, bcells, 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)
-    [gres_mag, gan_mag, ger_mag, ger_vec, ger_rel_mag, ger_rel_vec] = extractVectorInfo(gnum, ga)
-    [dgres_mag, dgan_mag, dger_mag, dger_vec, dger_rel_mag, dger_rel_vec] = extractVectorInfo(dgnum, dga)
+    [deres_mag, dean_mag, deer_mag, deer_vec, deer_rel_mag, deer_rel_vec] =  extractVectorInfo(denumperp,dea)
+    [eres_mag, ean_mag, eer_mag, eer_vec, eer_rel_mag, eer_rel_vec] = extractVectorInfo(enumperp, ea)
+    [gres_mag, gan_mag, ger_mag, ger_vec, ger_rel_mag, ger_rel_vec] = extractVectorInfo(gnumperp, ga)
+    [dgres_mag, dgan_mag, dger_mag, dger_vec, dger_rel_mag, dger_rel_vec] = extractVectorInfo(dgnumperp, dga)
+#    print('an',dgan_mag)
+ #   print('res', dgres_mag)
 
-    torsurf = pos[0][1]
-
-    torSliceDict = {}
-    keys = ['edrifts', 'efield','edrifts_perp','efield_perp' ,'ef_mag', 'ef_an_mag', 'edrift_mag', 'edrift_an_mag', 'ef_vec_er', 'ef_mag_er',
+    InfoDict = {}
+    keys = ['edrifts', 'efield','edrifts_perp','efield_perp' ,'ef_mag', 'ef_an_mag','ef_an', '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']
 
-    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]
+    values = [denum,enum,denumperp,enumperp,eres_mag,ean_mag,ea,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,gnum,dgnumperp,gnumperp,gres_mag,ean_mag,dgres_mag,dgan_mag,ger_vec,ger_mag,dger_vec,dger_mag,ger_rel_mag,ger_rel_vec,dger_rel_mag,dger_rel_vec,pos]
     for i in range(len(values)):
-        torSliceDict[keys[i]] = values[i]
-
-    return torSliceDict
+        InfoDict[keys[i]] = values[i]
+    pandasDF = pd.DataFrame.from_dict(data = InfoDict)
+    pandasDF.to_csv('AllInfo'+name+'.csv', header=True)
+    return InfoDict
 
 
 def readText(filepath):
@@ -130,14 +139,43 @@ def transformCyltoRadPerp(v,pos):
     phi = math.atan2(pos[2],pos[0]-500)
     c = math.cos(phi)
     s = math.sin(phi)
-    result.append(abs(v[0]*c+v[2]*s))
+    result.append(v[0]*c+v[2]*s)
     result.append(v[1])
-    result.append(abs(v[2]*c-v[0]*s))
-#    print('phi',phi,s,c)
- #   print('v',v)
-  #  print('res',result)
+    result.append(v[2]*c-v[0]*s)
+    return result
+    
+def readAllInfo(method,descr, folder):
+    result = {}
+#    print(folder+'AllInfo'+method+'Descr'+descr+'.csv')
+    df = pd.read_csv(folder+'AllInfo'+method+'Descr'+descr+'.csv')
+    keys = list(df.columns)
+    for i in range(1,len(keys)):
+        series = df.loc[:,keys[i]].to_list()
+        N = len(series)
+        res = np.array([None]*N)
+        for j in range(N):
+            if type(series[j]) is str:
+                res[j] = np.array(ast.literal_eval(series[j]))
+            else:
+                res[j] = series[j]
+        result[keys[i]] = res
+    return result
+
+def gradientLength(gradientmags,functionvalues):
+    result = []
+    for i in range(len(gradientmags)):
+        result.append(functionvalues[i]/(gradientmags[i]+1e-6))
     return result
 
+def transformRadPerptoCyl(v,pos):
+    result = []
+    phi = math.atan2(pos[2],pos[0]-500)
+    c = math.cos(phi)
+    s = math.sin(phi)
+    result.append(v[0]*c-v[2]*s)
+    result.append(v[1])
+    result.append(v[2]*c+v[0]*s)
+    return result
 
 
 
diff --git a/Scripts_Matthieu/AnalyticalPrescriptions.py b/Scripts_Matthieu/AnalyticalPrescriptions.py
index 7bc8e37..23048bd 100644
--- a/Scripts_Matthieu/AnalyticalPrescriptions.py
+++ b/Scripts_Matthieu/AnalyticalPrescriptions.py
@@ -4,21 +4,41 @@ import PlottingFunctions
 # Magnetic Fields
 def b_field_description1(r,z,theta, Rmajor, Rminor):
     # Radial B-field (linear and max at central surfaces of the torus
-    R = math.sqrt((r-Rmajor)**2 + z**2)+0.01+0*theta
-    bval = 0.01/R*100  # R in cm
+    r = r/100
+    Rmajor = Rmajor/100
+    z = z/100
+    R = math.sqrt((r-Rmajor)**2 + z**2)
+    bval = 1/R  # R in cm
     return [0,bval,0]
+    
 def bfieldstrength1(r,z,theta,Rmajor,Rminor):
     # Radial B-field (linear and max at central surfaces of the torus
-    R = math.sqrt((r-Rmajor)**2 + z**2)+0.01+0*theta
-    bval = 0.01/R*100  # R in cm
+    r = r/100
+    Rmajor = Rmajor/100
+    z = z/100
+    R = math.sqrt((r-Rmajor)**2 + z**2)
+    bval = 1/R  # R in cm
     return bval
 
 def b_field_description2(r,z,theta, Rmajor, Rminor):
     # Radial B-field (linear and max at central surfaces of the torus
-    R = math.sqrt((r-Rmajor)**2 + z**2)+0.01+0*theta
-    bval = 0.01/R*100  # R in cm
+    r = r/100
+    Rmajor = Rmajor/100
+    z = z/100
+    R = math.sqrt((r-Rmajor)**2 + z**2)
+    bval = 1/R # R in cm
     return [0,bval,0]
 
+def gradb_description1(r,z,theta,Rmajor,Rminor):
+    r = r/100
+    Rmajor = Rmajor/100
+    z = z/100
+    R = math.sqrt((r-Rmajor)**2 + z**2)
+    gBr = -1/R**2*(r-Rmajor)/R
+    gBtheta = 0
+    gBz = -1/R**2*z/R
+    return [gBr, gBtheta, gBz]
+
 # Potentials
 def potential_description1(r,z,theta,Rmajor,Rminor):
     r = r/100
@@ -28,15 +48,15 @@ def potential_description1(r,z,theta,Rmajor,Rminor):
     potval = math.exp(-D)*500
     return potval
 
-def potential_description2(r,z,theta,Rmajor,Rminor):
+def potential_description3(r,z,theta,Rmajor,Rminor):
     R = (r-Rmajor)/100+0.001
     potval = math.exp(-R)*100
     return potval
 
-def potential_description3(r,z,theta,Rmajor,Rminor):
+def potential_description2(r,z,theta,Rmajor,Rminor):
     r = (r-Rmajor)/100
     z = z/100
-    potval = 500*math.cos(math.atan2(z,r))
+    potval = 100*(math.cos(2*math.atan2(z,r)) + math.sqrt(r**2+z**2))
     return potval
 # Electric Fields
 def e_descr1(r,z,theta, Rmajor):
@@ -49,23 +69,19 @@ def e_descr1(r,z,theta, Rmajor):
     Ez = -500*math.exp(-D)*1/D*z
     return [Er,Ephi, Ez]
 
-def e_descr2(r,z,theta, Rmajor):
+def e_descr3(r,z,theta, Rmajor):
     D = (r - Rmajor)/100+0.001
     Er = -100*math.exp(-D)
     Ephi = 0
     Ez = 0
     return [Er,Ephi, Ez]
 
-def e_descr3(r,z,theta, Rmajor,rminor):
+def e_descr2(r,z,theta, Rmajor):
     r = (r-Rmajor)/100
     z = z/100
-    Er = 2000*math.sin(4*math.atan2(z,r))*z/(z**2+r**2)
+    Er = 200*math.sin(2*math.atan2(z,r))*z/(z**2+r**2) +100*r/math.sqrt(r**2+z**2)
     Ephi = 0
-    Ez = 2000*math.sin(4*math.atan2(z,r))*-r/(z**2+r**2)
-#    return [Er,Ephi,Ez]
-    return Ez
+    Ez = 200*math.sin(2*math.atan2(z,r))*-r/(z**2+r**2) + 100*z/math.sqrt(r**2+z**2)
+    return [Er,Ephi,Ez]
+#    return Ez
 
-Rmajor = 500
-rminor = 100
-ns = 100
-PlottingFunctions.plot2DContour(ns,ns,Rmajor-rminor,Rmajor+rminor,-rminor,rminor,Rmajor,rminor,e_descr3,'X [cm]', 'Y [cm]','Potential [V]',True)
\ No newline at end of file
diff --git a/Scripts_Matthieu/AnalyzeNumericalResults.py b/Scripts_Matthieu/AnalyzeNumericalResults.py
index adf9695..e7beb07 100644
--- a/Scripts_Matthieu/AnalyzeNumericalResults.py
+++ b/Scripts_Matthieu/AnalyzeNumericalResults.py
@@ -9,14 +9,14 @@ import PlottingFunctions
 import AnalysisFunctions
 
 # Constants
-path = '/u/mbj/Drifts/data-for-drift-computations/Description2/Case50/OutputGG/'
+path = '/u/mbj/Drifts/data-for-drift-computations/ResolutionStudies/RadRes/Case50/OutputLSQDescr1/'
 edrifts_file = open(path+'ExB_RESULTS','r')
 driftpos_file = open(path+'DRIFT_POS','r')
 efield_file = open(path+'EFIELD','r')
 gradb_file = open(path+'GRADB', 'r')
 gdrifts_file = open(path+'GRADB_DRIFTS', 'r')
 edrifts = []
-edrift_pos = []
+drift_pos = []
 efield = []
 gradb = []
 gdrifts = []
@@ -41,7 +41,9 @@ for line in gradb_file:
     gradb.append(line)
 
 # Convert values from fortran to float, filter boundary cells and compute analytical profiles for cell center locations
-[pos, enum, denum,gnum,dgnum, enumperp,denumperp, gnumperp ,dgnumperp,ea , ba, dea, ga,dga] = AnalysisFunctions.extractValues(drift_pos,drifts, efield,AnalyticalPrescriptions.e_descr3,AnalyticalPrescriptions.b_field_description1,driftFunctions.ExBdrift)
+[pos, enum, denum,gnum,dgnum, enumperp,denumperp, gnumperp ,dgnumperp,ea , ba, dea, ga,dga] = AnalysisFunctions.extractValues(drift_pos,edrifts,
+ efield, gradb, gdrifts, AnalyticalPrescriptions.e_descr3,AnalyticalPrescriptions.b_field_description1,
+ AnalyticalPrescriptions.gradb_description1, driftFunctions.ExBdrift, driftFunctions.gradBdrift)
 
 # Compute vectors and magnitude + errors for efield and drift
 [edrift_res_mag,edrift_an_mag,edrift_er_mag,edrift_er_vec,edrift_er_rel_mag,edrift_er_rel_vec] = AnalysisFunctions.extractVectorInfo(denum,dea)
@@ -57,25 +59,24 @@ 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']
+keys = ['edrifts', 'efield','edrifts_perp','efield_perp' ,'edrift_mag', 'edrift_an_mag','edrift_er_mag', 'edrift_vec_er','edrift_er_rel_mag', 'edrift_er_rel_vec',
+        'ef_mag', 'ef_an_mag', 'ef_mag_er', 'ef_vec_er','efield_er_rel_mag', 'efield_er_rel_vec', 
+        'gdrifts', 'gradb','gdrifts_perp','gradb_perp','gdrift_mag','gdrift_an_mag','gdrift_er_mag', 'gdrift_vec_er' ,'gdrift_er_rel_mag', 'gb_er_rel_vec',
+	'gb_mag','gb_an_mag', 'gb_mag_er', 'gb_vec_er','gb_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,
+allInfo = [denum, enum,denumperp,enumperp, 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,
+           dgnum,gnum, dgnumperp, gnumperp, 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]
 
-print(torSliceDict.get('ef_mag'))
-print(torSliceDict.get('ef_an_mag'))
-print(torSliceDict.get('efield_er_rel_mag'))
-print('max',max(torSliceDict.get('efield_er_rel_mag')))
+#print(torSliceDict.get('ef_mag'))
+#print(torSliceDict.get('ef_an_mag'))
+#print(torSliceDict.get('efield_er_rel_mag'))
+#print('max',max(torSliceDict.get('efield_er_rel_mag')))
 #print(len(torSliceDict.get('efield_er_rel_mag')))
 # Plot results
 
@@ -160,5 +161,5 @@ sc9 = plt.scatter(x, y, s=20, c= [x[2] for x in torSliceDict.get('efield_perp')]
 plt.colorbar(sc9)
 plt.savefig(path+title+'.png')
 
-#plt.show()
+plt.show()
 
diff --git a/Scripts_Matthieu/DoAnalysis.sh b/Scripts_Matthieu/DoAnalysis.sh
new file mode 100755
index 0000000..42d42b4
--- /dev/null
+++ b/Scripts_Matthieu/DoAnalysis.sh
@@ -0,0 +1,46 @@
+#!/bin/bash -l
+# Standard output and error:
+#SBATCH -o ./tjob.%j
+#SBATCH -e ./tjob.%j
+# Initial working directory:
+#SBATCH -D ./
+# Job Name:
+#SBATCH -J GeometriesforResolutionStudy
+# Queue (Partition):
+#SBATCH --partition=short
+# Don't specify 'SBATCH --nodes' !
+# Number of nodes and MPI tasks per node:
+#SBATCH --nodes=1
+#SBATCH --ntasks-per-node=1
+#SBATCH --mem=60G
+
+#SBATCH --mail-type=none
+#SBATCH --mail-user=user@mbj.mpg.de
+#
+# Wall clock limit:
+#SBATCH --time=00:10:00
+
+#python ResolutionStudyGeometries.py
+
+echo " Running Cases "
+
+cd ../../data-for-drift-computations/ResolutionStudies/
+
+for d in */ ; do
+    cd ./$d
+    echo "$d"
+    echo "$(pwd)"
+    declare -a DIRS=()
+    for dir in */ ; do
+	DIRS[${#DIRS[@]}]="$dir"
+    done
+    echo "${DIRS[*]}"
+    for case in 1 2 ; do
+        echo "Doing computations for description $case"
+	python ../../../python-scripts/Scripts_Matthieu/ErrorAnalysis.py $case ${DIRS[*]}
+    done
+    cd ..
+done
+
+echo "done all"
+    
diff --git a/Scripts_Matthieu/MeshFunctions.py b/Scripts_Matthieu/MeshFunctions.py
index 4f0ced4..fa6c54d 100644
--- a/Scripts_Matthieu/MeshFunctions.py
+++ b/Scripts_Matthieu/MeshFunctions.py
@@ -30,14 +30,13 @@ def create_toroidal_vertices(nb_tor,nb_rad,nb_pol,toroidal_range, R, Rmin):
     for i in range(nb_tor):
         for j in range(nb_pol):
             for k in range(nb_rad):
-                radial_positions[i,j*nb_rad+k] = R+r[k]*math.cos(j*2*math.pi/nb_pol)
-                horizontal_positions[i,j*nb_rad+k] = r[k]*math.sin(j*2*math.pi/nb_pol)
+                radial_positions[i,j*nb_rad+k] = R+r[k]*math.cos(j*2*math.pi/(nb_pol-1))
+                horizontal_positions[i,j*nb_rad+k] = r[k]*math.sin(j*2*math.pi/(nb_pol-1))
     return [toroidal_surfaces,radial_positions,horizontal_positions]
 
 
-def writeGEOMETRY3DDATA(tor_surf,rad_pos,hor_pos, ntheta, nr, nz, path):
+def writeGEOMETRY3DDATA(tor_surf,rad_pos,hor_pos, ntheta, nr, nz, filepath):
     # Open file to write the geometry information
-    filepath = path+"\\GEOMETRY_3D_DATA"
     gd = open(filepath, "w+")
     gd.write(str(nr)+"   ")
     gd.write(str(nz)+"   ")
@@ -138,7 +137,7 @@ def CreateGeneralInfo(filepath, nz,r,z,theta):
         gi.write("\n")
         gi.write("2")
         gi.write("\n")
-        gi.write("3  "+str(i)+"  1")
+        gi.write("0  "+str(i)+"  1")
         gi.write("\n")
         gi.write("0  "+str(z[i]-2)+"  0  "+str(theta[i]-2))
         gi.write("\n")
diff --git a/Scripts_Matthieu/PrepareGeometry.py b/Scripts_Matthieu/PrepareGeometry.py
index 392f68f..93a9c49 100644
--- a/Scripts_Matthieu/PrepareGeometry.py
+++ b/Scripts_Matthieu/PrepareGeometry.py
@@ -4,11 +4,11 @@ from MeshFunctions import *
 import os
 import shutil
 # The folder where the different cases should be stored is
-folder = "D:\Matthieu\Documents\\1.Master\Thesis\Data"
+folder = '/u/mbj/Drifts/data-for-drift-computations'
 # A specific folder for the case should be defined
 casefolder = "TestCase"
 # create the folder
-path = folder+"\\"+casefolder
+path = folder+"/"+casefolder
 # Check if the folder should be removed => Control to not loose simulations
 if os.path.isdir(path):
     Flag = int(input("Do you want to remove"+path+" Answer should be 1/0"))
@@ -17,46 +17,35 @@ if os.path.isdir(path):
     else:
         path = folder+"\\"+input("Give a new folder name")
 os.mkdir(path)
-path = path+"\\Inputs"
+path = path+"/Inputs"
 os.mkdir(path)
 # Choose the parameters for the case to test
 # Geometry
 # Number of zones considered
 nz = 1
 # Number of toroidal, radial, poloidal surfaces
-nt = 20
-nr = 10*10
-np = 10*5
+nt = 5
+nr = 10
+np = 10
 # toroidal range (degrees)
 trange = 36
 # Major and minor radius (considering a torus) (in cm!)
 Rmajor = 500
 Rminor = 100
 # Create the file with the general info
-gi_path = path+"\\input.geo"
+gi_path = path+"/input.geo"
 CreateGeneralInfo(gi_path, nz,[nr],[np],[nt])
 
-
+bfielddescription = AnalyticalPrescriptions.bfieldstrength1
 # Create the GEOMETRY_3D_DATA file
 [tor,rad,hor] = create_toroidal_vertices(nt,nr,np,trange,Rmajor,Rminor)
 print(tor)
-writeGEOMETRY3DDATA(tor,rad,hor,nt,nr,np,path)
+geopath = path+'/GEOMETRY_3D_DATA'
+writeGEOMETRY3DDATA(tor,rad,hor,nt,nr,np,geopath)
 
 # Create the magnetic field data
 # Create an analytical description
-def b_field_description(r,z,theta, Rmajor, Rminor):
-    # Radial B-field (linear and max at central surfaces of the torus
-    R = math.sqrt((r-Rmajor)**2 + z**2)+0.01+0*theta
-    bval = 0.01/R*100  # R in cm
-    return bval
-bfieldpath = path+"\\BFIELD_STRENGTH"
-imposeAnalytic(tor,rad,hor,b_field_description,Rmajor, Rminor, bfieldpath)
+bfieldpath = path+"/BFIELD_STRENGTH"
+imposeAnalytic(tor,rad,hor,bfielddescription,Rmajor, Rminor, bfieldpath)
+
 
-# Create the potential data
-# At the moment we use the same grid for magnetic and plasma cells
-def potential_description(r,z,theta,Rmajor,Rminor):
-    R = (r-Rmajor)
-    potval = math.exp(-R/100)*100
-    return potval
-potentialpath = path+"\\POTENTIAL"
-imposeAnalytic(tor,rad,hor,potential_description,Rmajor, Rminor,potentialpath)
diff --git a/Scripts_Matthieu/PrepareGeometryCompareNodal.py b/Scripts_Matthieu/PrepareGeometryCompareNodal.py
new file mode 100644
index 0000000..41d3637
--- /dev/null
+++ b/Scripts_Matthieu/PrepareGeometryCompareNodal.py
@@ -0,0 +1,51 @@
+from MeshFunctions import *
+import os
+import shutil
+import AnalyticalPrescriptions
+# The folder where the different cases should be stored is
+basepath = '/u/mbj/Drifts/data-for-drift-computations'
+folder = '/u/mbj/Drifts/data-for-drift-computations/CompareNodal/'
+os.mkdir(folder)
+# Choose the parameters for the case to test
+# Geometry
+# Number of zones considered
+nz = 1
+# Loop over different cases to prepare
+nts = [9, 18, 36]
+nrs = [15, 30, 60]
+nps = [100, 200, 400]
+#cases = [nrs,nps]
+caseNames = ['Case1','Case2','Case3']
+# Number of toroidal, radial, poloidal surfaces
+bfielddescription = AnalyticalPrescriptions.bfieldstrength1
+for i in range(len(nps)):
+    # A specific folder for the case should be defined
+    path = folder + "/" + caseNames[i]
+    os.mkdir(path)
+    path = path + "/Inputs"
+    os.mkdir(path)
+    nt = nts[i]
+    nr = nrs[i]
+    np = nps[i]
+    # toroidal range (degrees)
+    trange = 36
+    # Major and minor radius (considering a torus) (in cm!)
+    Rmajor = 500
+    Rminor = 100
+    # Create the file with the general info
+    gi_path = path+"/input.geo"
+    CreateGeneralInfo(gi_path, nz,[nr],[np],[nt])
+
+
+    # Create the GEOMETRY_3D_DATA file
+    [tor,rad,hor] = create_toroidal_vertices(nt,nr,np,trange,Rmajor,Rminor)
+    print(tor)
+    filepath = path+'/GEOMETRY_3D_DATA'
+    writeGEOMETRY3DDATA(tor,rad,hor,nt,nr,np,filepath)
+
+    # Create the magnetic field data
+    # Create an analytical description
+    bfieldpath = path+"/BFIELD_STRENGTH"
+    imposeAnalytic(tor,rad,hor,bfielddescription,Rmajor, Rminor, bfieldpath)
+
+  
diff --git a/Scripts_Matthieu/PrepareGeometryDescription1.py b/Scripts_Matthieu/PrepareGeometryDescription1.py
index 6e8a832..e275eb0 100644
--- a/Scripts_Matthieu/PrepareGeometryDescription1.py
+++ b/Scripts_Matthieu/PrepareGeometryDescription1.py
@@ -3,40 +3,42 @@ import os
 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\TorSurf"
+basepath = '/u/mbj/Drifts/data-for-drift-computations'
+folder = '/u/mbj/Drifts/data-for-drift-computations/StudyLSQ'
 if os.path.isdir(folder):
     Flag = int(input("Do you want to remove" + folder + " Answer should be 1/0"))
     if Flag:
         shutil.rmtree(folder)
     else:
-        folder = basepath + "\\" + input("Give a new folder name")
+        folder = basepath + "/" + input("Give a new folder name")
 os.mkdir(folder)
 # Choose the parameters for the case to test
 # Geometry
 # Number of zones considered
 nz = 1
 # Loop over different cases to prepare
-nts = [5,10,20,40,80,160]
-nr = 200
-np= nr
+nts = [9, 18, 36, 48, 72]
+nrs = [12, 25, 50, 75, 100]
+nps = [38, 75, 150, 225, 300]
 #cases = [nrs,nps]
-caseNames = ['Case5','Case10', 'Case20', 'Case40','Case80', 'Case160']
+caseNames = ['Case1', 'Case2', 'Case3','Case4', 'Case5']
 # Number of toroidal, radial, poloidal surfaces
 
 for i in range(len(nts)):
     # A specific folder for the case should be defined
-    path = folder + "\\" + caseNames[i]
+    path = folder + "/" + caseNames[i]
     if os.path.isdir(path):
         Flag = int(input("Do you want to remove" + path + " Answer should be 1/0"))
         if Flag:
             shutil.rmtree(path)
         else:
-            path = folder + "\\" + input("Give a new folder name")
+            path = folder + "/" + input("Give a new folder name")
     os.mkdir(path)
-    path = path + "\\Inputs"
+    path = path + "/Inputs"
     os.mkdir(path)
     nt = nts[i]
+    nr = nrs[i]
+    np = nps[i]
     potentialdescription = AnalyticalPrescriptions.potential_description1
     bfielddescription = AnalyticalPrescriptions.bfieldstrength1
     # toroidal range (degrees)
@@ -45,22 +47,23 @@ for i in range(len(nts)):
     Rmajor = 500
     Rminor = 100
     # Create the file with the general info
-    gi_path = path+"\\input.geo"
+    gi_path = path+"/input.geo"
     CreateGeneralInfo(gi_path, nz,[nr],[np],[nt])
 
 
     # Create the GEOMETRY_3D_DATA file
     [tor,rad,hor] = create_toroidal_vertices(nt,nr,np,trange,Rmajor,Rminor)
     print(tor)
-    writeGEOMETRY3DDATA(tor,rad,hor,nt,nr,np,path)
+    filepath = path+'/GEOMETRY_3D_DATA'
+    writeGEOMETRY3DDATA(tor,rad,hor,nt,nr,np,filepath)
 
     # Create the magnetic field data
     # Create an analytical description
-    bfieldpath = path+"\\BFIELD_STRENGTH"
+    bfieldpath = path+"/BFIELD_STRENGTH"
     imposeAnalytic(tor,rad,hor,bfielddescription,Rmajor, Rminor, bfieldpath)
 
     # Create the potential data
     # At the moment we use the same grid for magnetic and plasma cells
 
-    potentialpath = path+"\\POTENTIAL"
+    potentialpath = path+"/POTENTIAL"
     imposeAnalytic(tor,rad,hor,potentialdescription,Rmajor, Rminor,potentialpath)
diff --git a/Scripts_Matthieu/PrepareGeometryOptimalResolution.py b/Scripts_Matthieu/PrepareGeometryOptimalResolution.py
new file mode 100644
index 0000000..5066ef8
--- /dev/null
+++ b/Scripts_Matthieu/PrepareGeometryOptimalResolution.py
@@ -0,0 +1,51 @@
+from MeshFunctions import *
+import os
+import shutil
+import AnalyticalPrescriptions
+# The folder where the different cases should be stored is
+basepath = '/u/mbj/Drifts/data-for-drift-computations'
+folder = '/u/mbj/Drifts/data-for-drift-computations/FinalNoNoise/'
+os.mkdir(folder)
+# Choose the parameters for the case to test
+# Geometry
+# Number of zones considered
+nz = 1
+# Loop over different cases to prepare
+nts = [5, 9, 18, 36, 72]
+nrs = [8, 15, 30, 60, 120]
+nps = [50, 100, 200, 400, 800]
+#cases = [nrs,nps]
+caseNames = ['Case5_8_50','Case9_15_100','Case18_30_200','Case36_60_400','Case72_120_800']
+# Number of toroidal, radial, poloidal surfaces
+bfielddescription = AnalyticalPrescriptions.bfieldstrength1
+for i in range(len(nps)):
+    # A specific folder for the case should be defined
+    path = folder + "/" + caseNames[i]
+    os.mkdir(path)
+    path = path + "/Inputs"
+    os.mkdir(path)
+    nt = nts[i]
+    nr = nrs[i]
+    np = nps[i]
+    # toroidal range (degrees)
+    trange = 36
+    # Major and minor radius (considering a torus) (in cm!)
+    Rmajor = 500
+    Rminor = 100
+    # Create the file with the general info
+    gi_path = path+"/input.geo"
+    CreateGeneralInfo(gi_path, nz,[nr],[np],[nt])
+
+
+    # Create the GEOMETRY_3D_DATA file
+    [tor,rad,hor] = create_toroidal_vertices(nt,nr,np,trange,Rmajor,Rminor)
+    print(tor)
+    filepath = path+'/GEOMETRY_3D_DATA'
+    writeGEOMETRY3DDATA(tor,rad,hor,nt,nr,np,filepath)
+
+    # Create the magnetic field data
+    # Create an analytical description
+    bfieldpath = path+"/BFIELD_STRENGTH"
+    imposeAnalytic(tor,rad,hor,bfielddescription,Rmajor, Rminor, bfieldpath)
+
+  
diff --git a/Scripts_Matthieu/ResolutionStudyGeometries.py b/Scripts_Matthieu/ResolutionStudyGeometries.py
new file mode 100644
index 0000000..6641169
--- /dev/null
+++ b/Scripts_Matthieu/ResolutionStudyGeometries.py
@@ -0,0 +1,82 @@
+from MeshFunctions import *
+import os
+import shutil
+import AnalyticalPrescriptions
+import InputsGeometry
+# The folder where the different cases should be stored is
+basepath = '/u/mbj/Drifts/data-for-drift-computations'
+folder = '/u/mbj/Drifts/data-for-drift-computations/ResolutionStudies'
+os.mkdir(folder)
+# Choose the parameters for the case to test
+# Geometry
+# Number of zones considered
+nz = 1
+# Analytical Profiles
+potentialdescription = AnalyticalPrescriptions.potential_description1
+bfielddescription = AnalyticalPrescriptions.bfieldstrength1
+# General Info
+# toroidal range (degrees)
+trange = 36
+# Major and minor radius (considering a torus) (in cm!)
+Rmajor = 500
+Rminor = 100
+
+# Loop over different cases to prepare
+nts = [10, 20, 40, 80, 160]
+nrs = [12, 25, 50, 100, 200]
+nps = [40, 80, 160, 320, 640]
+
+#cases = [nrs,nps]
+caseNamesT = ['Case10', 'Case20', 'Case40','Case80', 'Case160']
+caseNamesP = ['Case40', 'Case80', 'Case160', 'Case320', 'Case640']
+caseNamesR = ['Case12', 'Case25', 'Case50', 'Case100', 'Case200']
+# Number of toroidal, radial, poloidal surfaces
+
+
+# First toroidal resolution
+nr = nrs[2]
+np = nps[2]
+path = folder + "/TorRes"
+os.mkdir(path)
+for i in range(len(nts)):
+    # A specific folder for the case should be defined
+    path = folder + "/TorRes/" + caseNamesT[i]
+    os.mkdir(path)
+    path = path + "/Inputs"
+    os.mkdir(path)
+    nt = nts[i]
+    
+    InputsGeometry.createInputFiles(nz, nr, np, nt,path, trange, Rmajor, Rminor,
+    bfielddescription, potentialdescription)
+
+
+nr = nrs[2]
+nt = nts[2]  
+path = folder + "/PolRes"
+os.mkdir(path)     
+for i in range(len(nps)):
+    # A specific folder for the case should be defined
+    path = folder + "/PolRes/" + caseNamesP[i]
+    os.mkdir(path)
+    path = path + "/Inputs"
+    os.mkdir(path)
+    np = nps[i]
+    InputsGeometry.createInputFiles(nz, nr, np, nt,path, trange, Rmajor, Rminor,
+    bfielddescription, potentialdescription)
+
+
+nt = nts[2]
+np = nps[2]   
+path = folder + "/RadRes"
+os.mkdir(path)
+for i in range(len(nrs)):
+    # A specific folder for the case should be defined
+    path = folder + "/RadRes/" + caseNamesR[i]
+    os.mkdir(path)
+    path = path + "/Inputs"
+    os.mkdir(path)
+    nr = nrs[i]
+    
+    InputsGeometry.createInputFiles(nz, nr, np, nt,path, trange, Rmajor, Rminor,
+    bfielddescription, potentialdescription)
+   
diff --git a/Scripts_Matthieu/driftFunctions.py b/Scripts_Matthieu/driftFunctions.py
index 9391799..56b209c 100644
--- a/Scripts_Matthieu/driftFunctions.py
+++ b/Scripts_Matthieu/driftFunctions.py
@@ -2,14 +2,12 @@ import numpy as np
 import math
 def ExBdrift(E,B):
     E = np.array(E)
- #   print('E', E)
     B = np.array(B)
-#    print('B',B)
     normB = np.linalg.norm(B)
     return np.cross(E,B)*1/normB**2
 
 def gradBdrift(B,gradB):
-    Tq = 150 # keV => no e for electroncharge (should be 000)
+    Tq = 1 # keV => no e for electroncharge (should be 000)
     return 2*Tq * np.cross(B,gradB)/np.linalg.norm(B)**3
 
 def cyltocart(V,p):
@@ -42,4 +40,4 @@ def carttocyl(V):
     C[0] = math.sqrt(V[0]**2+V[1]**2)
     C[1] = math.atan2(V[1],V[0])
     C[2] = V[2]
-    return C
\ No newline at end of file
+    return C
diff --git a/Scripts_Matthieu/test.py b/Scripts_Matthieu/test.py
new file mode 100644
index 0000000..99a93fd
--- /dev/null
+++ b/Scripts_Matthieu/test.py
@@ -0,0 +1,31 @@
+import math
+import driftFunctions
+import AnalyticalPrescriptions
+import Constants
+import Utilities
+import pandas as pd
+import matplotlib.pyplot as plt
+import PlottingFunctions
+import AnalysisFunctions
+import sys
+import os
+import seaborn as sns
+
+method = 'LSQ'
+anCase = '1'
+#Case = sys.argv[3]
+path = '' # Case+'/'
+allInfo = []
+
+plot = plt.figure(1)
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+limits = [1e3, 1e7, 1e-7, 1]
+nbs = [1e4, 1e5, 1e6]
+MaxErrors = [1e-1, 1e-3, 1e-5]
+plt.axis(limits)
+sc1 = plt.scatter(nbs,MaxErrors,label = 'Max Error ')
+vl = plt.vlines([2000, 5000, 10000],1e-7, 1,colors = ['g','r','b'])
+plt.show()
+
-- 
GitLab