From 0637ee43e9ad608c3a5c882b7bb82fad84a1faab Mon Sep 17 00:00:00 2001
From: Matthieu Jacobs <matthieu.jacobs@student.kuleuven.be>
Date: Wed, 11 May 2022 11:13:39 +0200
Subject: [PATCH] Modified sripts

---
 GridQuality/grid_analyze.py                   |  27 +-
 Scripts_Matthieu/AnalysisFunctions.py         |  33 +-
 Scripts_Matthieu/AnalyticalPrescriptions.py   |  10 +
 Scripts_Matthieu/Constants.py                 |   2 +-
 Scripts_Matthieu/ErrorAnalysis.py             | 520 ++++++++++++++++--
 Scripts_Matthieu/ExtractNoiseInfo.py          |  23 +-
 Scripts_Matthieu/PlotTrajectorySingle.py      |  32 +-
 .../PrepareGeometryCompareNodal.py            |   8 +-
 Scripts_Matthieu/TrajectoryDeviation.py       |  13 +-
 Scripts_Matthieu/writeInfo.py                 |   6 +-
 10 files changed, 572 insertions(+), 102 deletions(-)

diff --git a/GridQuality/grid_analyze.py b/GridQuality/grid_analyze.py
index 31fc5fa..5be17c1 100644
--- a/GridQuality/grid_analyze.py
+++ b/GridQuality/grid_analyze.py
@@ -1,6 +1,7 @@
 import numpy as np
 from matplotlib import pyplot as plt
-from combine.quick_g import quick_g_plot
+from quick_g import quick_g_plot
+import sys
 
 # not whole loop, only tenth -> phi doesnt loop    
 def yield_grid_cells(g, r_lahead=1, tht_lahead=1, phi_lahead=0):
@@ -184,7 +185,8 @@ def nonconvex(g):
 
 if __name__ == "__main__":
     try:
-        g = np.load("g_possible.npy")
+        name = sys.argv[1]
+        g = np.load(name)
         # g = np.load("g_kaputt.npy")
     except:
         print("Need to generate grid array 'g.npy' with another program")
@@ -193,11 +195,16 @@ if __name__ == "__main__":
     
     #nonc = nonconvex(g) # TODO  1 == good, 0 == nonconvex, -1 == weird
     
-    # volumes = volume(g)
-    # angles_r, angles_tht = non_orthogonality(g)
-    # uneven_r, uneven_tht = unevenness(g)
-    # skewness_r, skewness_tht = skewness(g)
-    # print(f"{np.sum(volumes):.3E} <- Whole volume in m^3")
-    # print(f"{np.min(volumes):.3E}, {np.max(volumes):.3E} <- Min, Max cell volume in m^3")
-    # print(f"{np.mean(volumes):.3E}, {np.std(volumes):.3E} <- Mean, Std of cell volume in m^3")
-    
\ No newline at end of file
+    volumes = volume(g)
+    angles_r, angles_tht = non_orthogonality(g)
+    uneven_r, uneven_tht = unevenness(g)
+    skewness_r, skewness_tht = skewness(g)
+    print(f"{np.sum(volumes):.3E} <- Whole volume in m^3")
+    print(f"{np.min(volumes):.3E}, {np.max(volumes):.3E} <- Min, Max cell volume in m^3")
+    print(f"{np.mean(volumes):.3E}, {np.std(volumes):.3E} <- Mean, Std of cell volume in m^3")
+    print(f"{np.mean(uneven_r)} <- Uneven r")
+    print(f"{np.mean(uneven_tht)} <- Uneven theta")
+    print(f"{np.mean(skewness_r)} <- Skewness r")
+    print(f"{np.mean(skewness_tht)} <- Skewness theta")
+    print(f"{np.mean(angles_r)} <- Unorthogonality r")
+    print(f"{np.mean(angles_tht)} <- Unorthogonality theta")
diff --git a/Scripts_Matthieu/AnalysisFunctions.py b/Scripts_Matthieu/AnalysisFunctions.py
index af91043..8c8e67f 100644
--- a/Scripts_Matthieu/AnalysisFunctions.py
+++ b/Scripts_Matthieu/AnalysisFunctions.py
@@ -7,7 +7,7 @@ import pandas as pd
 import numpy as np
 import ast
 
-def extractValues(positions,gbpositions, drifts, efield, gradb, drgradb, bcells,bnodes, edescr, bdescr, gradbdescr, ExB, gradBdrift, ggcorr):
+def extractValues(positions,gbpositions, drifts, efield, gradb, drgradb, bcells,bnodes, edescr, bdescr, gradbdescr, ExB, gradBdrift, ggcorr, Rmajor):
     pos = []
     gbpos = []
     enum = []
@@ -40,24 +40,24 @@ def extractValues(positions,gbpositions, drifts, efield, gradb, drgradb, bcells,
 
         # filter our boundary cells => These have zero field and drift and are not useful for interpretation!
 	# quick fix with nan
-        if  bnd != 1 and not math.isnan(en[0]) and not math.isnan(den[0]): #and (abs(max(en))) <500 and abs(min(en))<500: 
+        if  bnd != 1 and not math.isnan(en[0]) and not math.isnan(den[0]) and (abs(max(en))) <300 and abs(min(en))<300: 
             denum.append(den)
             pos.append(dpn)
             enum.append(en)
             dgnum.append(dbn)
             # Compute analytical efield and drifts at these positions
-            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)
+            efa = edescr(dpn[0], dpn[2], dpn[1], Rmajor)
+            bfa = bdescr(dpn[0], dpn[2], dpn[1],Rmajor, Constants.Rminor)
+            gfa = gradbdescr(dpn[0], dpn[2], dpn[1], Rmajor, Constants.Rminor)
             dgfa = gradBdrift(bfa,gfa)
             dga.append(dgfa.tolist())
             ea.append(efa)
             ba.append(bfa)
             ggcorrnum.append(gg)
             dea.append(ExB(efa, bfa).tolist())
-            denumperp.append(transformRadPerptoCyl(denum[-1],pos[-1]))
-            enumperp.append(transformRadPerptoCyl(enum[-1],pos[-1]))
-            dgnumperp.append(transformRadPerptoCyl(dgnum[-1],pos[-1]))
+            denumperp.append(transformRadPerptoCyl(denum[-1],pos[-1], Rmajor))
+            enumperp.append(transformRadPerptoCyl(enum[-1],pos[-1], Rmajor))
+            dgnumperp.append(transformRadPerptoCyl(dgnum[-1],pos[-1], Rmajor))
 	    
     for i in range(len(gradb)):
         gbn = [float(x) for x in gradb[i].split()]
@@ -65,10 +65,10 @@ def extractValues(positions,gbpositions, drifts, efield, gradb, drgradb, bcells,
         dpn = [float(x) for x in gbpositions[i].split()]
         if  bnd != 1:
             gnum.append(gbn)
-            gfa = gradbdescr(dpn[0], dpn[2], dpn[1], Constants.Rmajor, Constants.Rminor)
+            gfa = gradbdescr(dpn[0], dpn[2], dpn[1], Rmajor, Constants.Rminor)
             ga.append(gfa)
             gbpos.append(dpn)
-            gnumperp.append(transformRadPerptoCyl(gnum[-1],gbpos[-1]))
+            gnumperp.append(transformRadPerptoCyl(gnum[-1],gbpos[-1],Rmajor))
     return [pos,gbpos, enum, denum,gnum,dgnum, enumperp,denumperp, gnumperp ,dgnumperp,ea , ba, dea, ga,dga, ggcorrnum]
 
 def extractVectorInfo(numeric,analytical):
@@ -115,9 +115,9 @@ def getAllSlices(positions,allInfo,direction,coordinate):
     result.append(pos_slice)
     return result
 
-def getAllInfo(positions,gbpositions,drifts,efield, gradb, drgradb,bcells,bnodes,edescr,bdescr, gradbdescr, ExB, gradBdrift,GGcorr, name):
+def getAllInfo(positions,gbpositions,drifts,efield, gradb, drgradb,bcells,bnodes,edescr,bdescr, gradbdescr, ExB, gradBdrift,GGcorr,Rmajor, name):
     [pos, gbpos, enum, denum,gnum,dgnum, enumperp,denumperp, gnumperp ,dgnumperp,ea , ba, dea, ga,dga, GGcorr] = extractValues(positions, gbpositions,
-                        drifts, efield, gradb, drgradb, bcells,bnodes, edescr, bdescr,gradbdescr, ExB, gradBdrift, GGcorr)
+                        drifts, efield, gradb, drgradb, bcells,bnodes, edescr, bdescr,gradbdescr, ExB, gradBdrift, GGcorr, Rmajor)
 
     [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)
@@ -155,9 +155,9 @@ def readText(filepath):
     return result
     
     
-def transformCyltoRadPerp(v,pos):
+def transformCyltoRadPerp(v,pos, Rmajor):
     result = []
-    phi = math.atan2(pos[2],pos[0]-500)
+    phi = math.atan2(pos[2],pos[0]-Rmajor)
     c = math.cos(phi)
     s = math.sin(phi)
     result.append(v[0]*c+v[2]*s)
@@ -167,6 +167,7 @@ def transformCyltoRadPerp(v,pos):
     
 def readAllInfo(method,descr,noise, folder):
     result = {}
+    print(folder+'AllInfo'+method+'Descr'+descr+'Noise'+noise+'.csv')
     df = pd.read_csv(folder+'AllInfo'+method+'Descr'+descr+'Noise'+noise+'.csv')
     keys = list(df.columns)
     for i in range(1,len(keys)):
@@ -201,9 +202,9 @@ def gradientLength(gradientmags,functionvalues):
         result.append(functionvalues[i]/(gradientmags[i]+1e-6))
     return result
 
-def transformRadPerptoCyl(v,pos):
+def transformRadPerptoCyl(v,pos,Rmajor):
     result = []
-    phi = math.atan2(pos[2],pos[0]-500)
+    phi = math.atan2(pos[2],pos[0]-Rmajor)
     c = math.cos(phi)
     s = math.sin(phi)
     result.append(v[0]*c-v[2]*s)
diff --git a/Scripts_Matthieu/AnalyticalPrescriptions.py b/Scripts_Matthieu/AnalyticalPrescriptions.py
index 793841f..da6f92c 100644
--- a/Scripts_Matthieu/AnalyticalPrescriptions.py
+++ b/Scripts_Matthieu/AnalyticalPrescriptions.py
@@ -92,3 +92,13 @@ def e_descr2(r,z,theta, Rmajor):
     return [Er,Ephi,Ez]
 #    return Ez
 
+def e_descr6(r,z,theta, Rmajor):
+    r = r/100
+    Rmajor = 590/100
+    z = z/100
+    D = math.sqrt((r - Rmajor) ** 2 + z ** 2)
+    Er = -500*math.exp(-D)*1/D*(r-Rmajor)
+    Ephi = 0
+    Ez = -500*math.exp(-D)*1/D*z
+    return [Er,Ephi, Ez]
+
diff --git a/Scripts_Matthieu/Constants.py b/Scripts_Matthieu/Constants.py
index 02dbe68..ff0c3c7 100644
--- a/Scripts_Matthieu/Constants.py
+++ b/Scripts_Matthieu/Constants.py
@@ -1,4 +1,4 @@
 # This file contains constant values, to be used in other files
 Rmajor = 500
 Rminor = 100
-ns = 100
\ No newline at end of file
+ns = 100
diff --git a/Scripts_Matthieu/ErrorAnalysis.py b/Scripts_Matthieu/ErrorAnalysis.py
index a1b8cd3..4c2ca94 100644
--- a/Scripts_Matthieu/ErrorAnalysis.py
+++ b/Scripts_Matthieu/ErrorAnalysis.py
@@ -21,8 +21,8 @@ Cases = sys.argv[2:]
 path = os.getcwd()
 print(path)
 print(Cases)
-methods = ['1', '2', '3', '4','5'] 
-
+methodnames = ['GG Cell-centered', 'LS Cell-centered', 'GG Nodal', 'LS Nodal','Mapping'] 
+methods = ['1','2','3','4','5']
 # Define structures to be computed    
 efield_error_key = 'efield_er_rel_mag'
 L2Errors = []
@@ -33,11 +33,9 @@ GGErrors = []
 lens = []
 
 nbs = []
-siner = 1000
-dmin = 1000000000
-r = 550
+r = 530
 phi = 0.15
-z = 30
+z = 45
 gg = 10
 for method in methods:
 #    errorsDF = pd.DataFrame()
@@ -48,7 +46,16 @@ for method in methods:
     charLength = []
     maxErrors = []
     errorsGG = []
+    distances = []
+    distancesnontor = []
+    singledistances = []
+    singledistancesnontor = []
+    
     for case in Cases:
+        siner = 1000
+        dmin = 1000000000
+        distmax = 0
+        distnontormax = 0
         print(case)
         [caseInfo, caseInfoGB] = AnalysisFunctions.readAllInfo(method, analyticalCase, '0', case)
         caseErrors = caseInfo.get(efield_error_key)
@@ -63,17 +70,41 @@ for method in methods:
         vol = float(vol[0])
 #        print(vol)
         charLength.append((vol/10**6)**(1/3))
-        print(charLength)
         maxErrors.append(max(caseErrors))
+# Read geometrical data
+        geoData = AnalysisFunctions.readText(path+'/'+case+'Output'+method+'Descr'+analyticalCase+'Noise0/NBS')	
 # Get single error
         for i in range(len(caseErrors)):
-            ipos = positions[i]
-            d = (z-ipos[2])**2+(r-ipos[0])**2+(phi-ipos[0])**2
-            if d < dmin :
-                siner = caseErrors[i] 
-                gg = ggErrors[i]
+           ipos = positions[i]
+           d = (z-ipos[2])**2+(r-ipos[0])**2+(phi-ipos[0])**2
+           if d<dmin:
+               siner = caseErrors[i] 
+               gg = ggErrors[i]
+               dmin = d
+               singlepos = ipos
+        for i in range(math.floor(len(geoData)/5)):
+            ipos = [float(x) for x in geoData[5*i+2].split()]
+            dis = [float(x) for x in geoData[5*i+3].split()] + [float(x) for x in geoData[5*i+4].split()]
+            bound = int(geoData[5*i+1])
+            dist = max(dis)
+#            print(bound, dist,distmax)
+#            print(bound == 0)
+#            print(type(bound))
+            distnontor = max(dis[0:4])
+            if dist > distmax and bound == 0:
+                distmax = dist
+            if distnontor > distnontormax and bound == 0:
+                distnontormax = distnontor
+            if np.array_equal(ipos,singlepos):
+                sindist = dist
+                sindistnontor = distnontor
+                
         errorsSingle.append(siner)
         errorsGG.append(gg)	
+        distances.append(distmax)
+        distancesnontor.append(distnontormax)
+        singledistances.append(sindist)
+        singledistancesnontor.append(sindistnontor)
 #        print(errors)
 #        print(maxErrors)
 #        print(nbCells)
@@ -94,47 +125,60 @@ print('Single', SingleErrors)
 print('Max',MaxErrors)
 print('GG',GGErrors)
 print(nbs) 
+print('distances', distances)
+print('distancesnontor', distancesnontor)
+print(singledistances)
+print(singledistancesnontor)
+print('char length', lens)
 # Compute measure for gradient lengths and typical grid size
 typsize = 50*36*150
 #maxGL = 
 #avGL =  
  
 
-upperX = 1e0
-lowerX = 1e-5
+upperX = 1e-1
+lowerX = 5*1e-3
 upperY = 1e-1
-lowerY = 1e-9   
+lowerY = 1e-6   
+
+##########################################################
+# Characteristic Length
+
+
 xlabel = 'Characteristic Cell Length [m]'
 ylabel = 'Relative Magnitude Error'
+xarray = lens[0]
 
 
 plot1 = plt.figure(1)
 v1 = plt.vlines(typsize,lowerY,upperY, color = 'g')
-sc1 = plt.scatter(lens[0],MeanErrors[0],label = methods[0])
-sc2 = plt.scatter(lens[1],MeanErrors[1],label = methods[1])
-sc3 = plt.scatter(lens[1],MeanErrors[2],label = methods[2]) 
-sc4 = plt.scatter(lens[1],MeanErrors[3],label = methods[3]) 
-sc5 = plt.scatter(lens[1],MeanErrors[4],label = methods[4]) 
+sc1 = plt.scatter(xarray,MeanErrors[0],label = methodnames[0])
+sc2 = plt.scatter(xarray,MeanErrors[1],label = methodnames[1])
+sc3 = plt.scatter(xarray,MeanErrors[2],label = methodnames[2]) 
+sc4 = plt.scatter(xarray,MeanErrors[3],label = methodnames[3]) 
+sc5 = plt.scatter(xarray,MeanErrors[4],label = methodnames[4]) 
 ax = plt.gca()
 ax.set_yscale('log')
 ax.set_xscale('log')
 limits = [lowerX, upperX, lowerY, upperY]
 ax.legend()
 x = np.linspace(0.000001,1,100)
-y = x**(1)/100
-plt.loglog(x,y)
+y = x**(1)/10
+y2 = x**2/10
+plt.loglog(x,y,label = 'first order')
+plt.loglog(x,y2, label = 'second order')
 plt.axis(limits)
-plt.title('Mean Errors Description'+analyticalCase)
+plt.title('Mean Errors Description '+analyticalCase)
 plt.xlabel(xlabel)
 plt.ylabel(ylabel)
-plt.savefig(path+'/Descr'+analyticalCase+'MeanError.png')
+plt.savefig(path+'/Descr'+analyticalCase+'MeanError'+xlabel+'.png')
 
 plot2 = plt.figure(2)
-sc1 = plt.scatter(lens[0],L2Errors[0],label = methods[0])
-sc2 = plt.scatter(lens[1],L2Errors[1],label = methods[1])
-sc3 = plt.scatter(lens[0],L2Errors[2],label = methods[2])
-sc4 = plt.scatter(lens[1],L2Errors[3],label = methods[3])
-sc5 = plt.scatter(lens[0],L2Errors[4],label = methods[4])
+sc1 = plt.scatter(xarray,L2Errors[0],label = methodnames[0])
+sc2 = plt.scatter(xarray,L2Errors[1],label = methodnames[1])
+sc3 = plt.scatter(xarray,L2Errors[2],label = methodnames[2])
+sc4 = plt.scatter(xarray,L2Errors[3],label = methodnames[3])
+sc5 = plt.scatter(xarray,L2Errors[4],label = methodnames[4])
 
 v2 = plt.vlines(typsize,lowerY,upperY, color = 'g')
 ax = plt.gca()
@@ -144,20 +188,103 @@ limits = [lowerX, upperX, lowerY, upperY]
 ax.legend()
 x = np.linspace(0.000001,1,100)
 y = x**(1)/100
-plt.loglog(x,y)
+y2 = x**2/100
+plt.loglog(x,y,label = 'first order')
+plt.loglog(x,y2, label = 'second order')
 plt.axis(limits)
 plt.title('L2 Errors Description'+analyticalCase)
 plt.xlabel(xlabel)
 plt.ylabel(ylabel)
-plt.savefig(path+'/Descr'+analyticalCase+'L2Error.png')
+plt.savefig(path+'/Descr'+analyticalCase+'L2Error'+xlabel+'.png')
 
 
 plot3 = plt.figure(3)
-sc1 = plt.scatter(nbs[0],MaxErrors[0],label = 'Max Error '+methods[0])
-sc2 = plt.scatter(lens[1],MaxErrors[1],label = 'Max Error '+methods[1])
-sc3 = plt.scatter(lens[0],MaxErrors[2],label = 'Max Error '+methods[2])
-sc4 = plt.scatter(lens[1],MaxErrors[3],label = 'Max Error '+methods[3])
-sc5 = plt.scatter(lens[0],MaxErrors[4],label = 'Max Error '+methods[4])
+sc1 = plt.scatter(xarray,MaxErrors[0],label = 'Max Error '+methodnames[0])
+sc2 = plt.scatter(xarray,MaxErrors[1],label = 'Max Error '+methodnames[1])
+sc3 = plt.scatter(xarray,MaxErrors[2],label = 'Max Error '+methodnames[2])
+sc4 = plt.scatter(xarray,MaxErrors[3],label = 'Max Error '+methodnames[3])
+sc5 = plt.scatter(xarray,MaxErrors[4],label = 'Max Error '+methodnames[4])
+
+v2 = plt.vlines(typsize,lowerY,upperY, color = 'g')
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+limits = [lowerX, upperX, lowerY, upperY]
+ax.legend()
+x = np.linspace(0.000001,1,100)
+y = x**(1)/10
+y2 = x**2/10
+plt.loglog(x,y,label = 'first order')
+plt.loglog(x,y2, label = 'second order')
+plt.axis(limits)
+plt.title('Max Errors Description'+analyticalCase)
+plt.xlabel(xlabel)
+plt.ylabel(ylabel)
+plt.savefig(path+'/Descr'+analyticalCase+'MaxError'+xlabel+'.png')
+
+plot4 = plt.figure(4)
+sc1 = plt.scatter(xarray,GGErrors[0],label = methodnames[0])
+sc2 = plt.scatter(xarray,GGErrors[1],label = methodnames[1])
+sc3 = plt.scatter(xarray,GGErrors[2],label = methodnames[2])
+sc4 = plt.scatter(xarray,GGErrors[3],label = methodnames[3])
+sc5 = plt.scatter(xarray,GGErrors[4],label = methodnames[4])
+
+v2 = plt.vlines(typsize,lowerY,upperY, color = 'g')
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+limits = [lowerX, upperX, lowerY, upperY]
+ax.legend()
+x = np.linspace(0.000001,1,100)
+y = x**(1)/10
+y2 = x**2/10
+or1 = plt.loglog(x,y,label = 'first order')
+or2 = plt.loglog(x,y2, label = 'second order')
+plt.axis(limits)
+plt.title('GG Correction Error'+analyticalCase)
+plt.xlabel(xlabel)
+plt.ylabel(ylabel)
+plt.savefig(path+'/Descr'+analyticalCase+'GG correction Error'+xlabel+'.png')
+
+
+xlabel = 'Max Cell Length [m]'
+ylabel = 'Relative Magnitude Error'
+xarray = distances
+upperX = 1e0
+lowerX = 5*1e-2
+
+######################################################################
+# Max distance
+
+plot5 = plt.figure(5)
+v1 = plt.vlines(typsize,lowerY,upperY, color = 'g')
+sc1 = plt.scatter(xarray,MeanErrors[0],label = methodnames[0])
+sc2 = plt.scatter(xarray,MeanErrors[1],label = methodnames[1])
+sc3 = plt.scatter(xarray,MeanErrors[2],label = methodnames[2]) 
+sc4 = plt.scatter(xarray,MeanErrors[3],label = methodnames[3]) 
+sc5 = plt.scatter(xarray,MeanErrors[4],label = methodnames[4]) 
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+limits = [lowerX, upperX, lowerY, upperY]
+ax.legend()
+x = np.linspace(0.000001,1,100)
+y = x**(1)/400
+y2 = x**2/400
+plt.loglog(x,y,label = 'first order')
+plt.loglog(x,y2, label = 'second order')
+plt.axis(limits)
+plt.title('Mean Errors Description'+analyticalCase)
+plt.xlabel(xlabel)
+plt.ylabel(ylabel)
+plt.savefig(path+'/Descr'+analyticalCase+'MeanError'+xlabel+'.png')
+
+plot6 = plt.figure(6)
+sc1 = plt.scatter(xarray,L2Errors[0],label = methodnames[0])
+sc2 = plt.scatter(xarray,L2Errors[1],label = methodnames[1])
+sc3 = plt.scatter(xarray,L2Errors[2],label = methodnames[2])
+sc4 = plt.scatter(xarray,L2Errors[3],label = methodnames[3])
+sc5 = plt.scatter(xarray,L2Errors[4],label = methodnames[4])
 
 v2 = plt.vlines(typsize,lowerY,upperY, color = 'g')
 ax = plt.gca()
@@ -167,20 +294,102 @@ limits = [lowerX, upperX, lowerY, upperY]
 ax.legend()
 x = np.linspace(0.000001,1,100)
 y = x**(1)/100
-plt.loglog(x,y)
+y2 = x**2/100
+plt.loglog(x,y,label = 'first order')
+plt.loglog(x,y2, label = 'second order')
+plt.axis(limits)
+plt.title('L2 Errors Description'+analyticalCase)
+plt.xlabel(xlabel)
+plt.ylabel(ylabel)
+plt.savefig(path+'/Descr'+analyticalCase+'L2Error'+xlabel+'.png')
+
+
+plot7 = plt.figure(7)
+sc1 = plt.scatter(xarray,MaxErrors[0],label = 'Max Error '+methodnames[0])
+sc2 = plt.scatter(xarray,MaxErrors[1],label = 'Max Error '+methodnames[1])
+sc3 = plt.scatter(xarray,MaxErrors[2],label = 'Max Error '+methodnames[2])
+sc4 = plt.scatter(xarray,MaxErrors[3],label = 'Max Error '+methodnames[3])
+sc5 = plt.scatter(xarray,MaxErrors[4],label = 'Max Error '+methodnames[4])
+
+v2 = plt.vlines(typsize,lowerY,upperY, color = 'g')
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+limits = [lowerX, upperX, lowerY, upperY]
+ax.legend()
+x = np.linspace(0.000001,1,100)
+y = x**(1)/200
+y2 = x**2/200
+plt.loglog(x,y,label = 'first order')
+plt.loglog(x,y2, label = 'second order')
 plt.axis(limits)
 plt.title('Max Errors Description'+analyticalCase)
 plt.xlabel(xlabel)
 plt.ylabel(ylabel)
-plt.savefig(path+'/Descr'+analyticalCase+'MaxError.png')
+plt.savefig(path+'/Descr'+analyticalCase+'MaxError'+xlabel+'.png')
 
+plot8 = plt.figure(8)
+sc1 = plt.scatter(xarray,GGErrors[0],label = methodnames[0])
+sc2 = plt.scatter(xarray,GGErrors[1],label = methodnames[1])
+sc3 = plt.scatter(xarray,GGErrors[2],label = methodnames[2])
+sc4 = plt.scatter(xarray,GGErrors[3],label = methodnames[3])
+sc5 = plt.scatter(xarray,GGErrors[4],label = methodnames[4])
 
-plot4 = plt.figure(4)
-sc1 = plt.scatter(lens[0],SingleErrors[0],label = methods[0])
-sc2 = plt.scatter(lens[1],SingleErrors[1],label = methods[1])
-sc3 = plt.scatter(lens[0],SingleErrors[2],label = methods[2])
-sc4 = plt.scatter(lens[1],SingleErrors[3],label = methods[3])
-sc5 = plt.scatter(lens[0],SingleErrors[4],label = methods[4])
+v2 = plt.vlines(typsize,lowerY,upperY, color = 'g')
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+limits = [lowerX, upperX, lowerY, upperY]
+ax.legend()
+x = np.linspace(0.000001,1,100)
+y = x**(1)/100
+y2 = x**2/100
+plt.loglog(x,y,label = 'first order')
+plt.loglog(x,y2, label = 'second order')
+plt.axis(limits)
+plt.title('GG Correction Error'+analyticalCase+xlabel)
+plt.xlabel(xlabel)
+plt.ylabel(ylabel)
+plt.savefig(path+'/Descr'+analyticalCase+'GG correction Error'+xlabel+'.png')
+
+############################################################################
+# Max length, non toroidal
+xlabel = 'Max Length (excluding toroidal direction) [m]'
+ylabel = 'Relative Magnitude Error'
+xarray = distancesnontor
+
+upperX = 1e-1
+lowerX = 5*1e-3 
+
+plot9 = plt.figure(9)
+v1 = plt.vlines(typsize,lowerY,upperY, color = 'g')
+sc1 = plt.scatter(xarray,MeanErrors[0],label = methodnames[0])
+sc2 = plt.scatter(xarray,MeanErrors[1],label = methodnames[1])
+sc3 = plt.scatter(xarray,MeanErrors[2],label = methodnames[2]) 
+sc4 = plt.scatter(xarray,MeanErrors[3],label = methodnames[3]) 
+sc5 = plt.scatter(xarray,MeanErrors[4],label = methodnames[4]) 
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+limits = [lowerX, upperX, lowerY, upperY]
+ax.legend()
+x = np.linspace(0.000001,1,100)
+y = x**(1)/10
+y2 = x**2/10
+plt.loglog(x,y,label = 'first order')
+plt.loglog(x,y2, label = 'second order')
+plt.axis(limits)
+plt.title('Mean Errors Description'+analyticalCase)
+plt.xlabel(xlabel)
+plt.ylabel(ylabel)
+plt.savefig(path+'/Descr'+analyticalCase+'MeanError'+xlabel+'.png')
+
+plot10 = plt.figure(10)
+sc1 = plt.scatter(xarray,L2Errors[0],label = methodnames[0])
+sc2 = plt.scatter(xarray,L2Errors[1],label = methodnames[1])
+sc3 = plt.scatter(xarray,L2Errors[2],label = methodnames[2])
+sc4 = plt.scatter(xarray,L2Errors[3],label = methodnames[3])
+sc5 = plt.scatter(xarray,L2Errors[4],label = methodnames[4])
 
 v2 = plt.vlines(typsize,lowerY,upperY, color = 'g')
 ax = plt.gca()
@@ -190,19 +399,22 @@ limits = [lowerX, upperX, lowerY, upperY]
 ax.legend()
 x = np.linspace(0.000001,1,100)
 y = x**(1)/100
-plt.loglog(x,y)
+y2 = x**2/100
+plt.loglog(x,y,label = 'first order')
+plt.loglog(x,y2, label = 'second order')
 plt.axis(limits)
-plt.title('Single Cell Error'+analyticalCase)
+plt.title('L2 Errors Description'+analyticalCase)
 plt.xlabel(xlabel)
 plt.ylabel(ylabel)
-plt.savefig(path+'/Descr'+analyticalCase+'SingleError.png')
+plt.savefig(path+'/Descr'+analyticalCase+'L2Error'+xlabel+'.png')
 
-plot5 = plt.figure(5)
-sc1 = plt.scatter(nbs[0],GGErrors[0],label = methods[0])
-sc2 = plt.scatter(lens[1],GGErrors[1],label = methods[1])
-sc3 = plt.scatter(lens[0],GGErrors[2],label = methods[2])
-sc4 = plt.scatter(lens[1],GGErrors[3],label = methods[3])
-sc5 = plt.scatter(lens[0],GGErrors[4],label = methods[4])
+
+plot11 = plt.figure(11)
+sc1 = plt.scatter(xarray,MaxErrors[0],label = 'Max Error '+methodnames[0])
+sc2 = plt.scatter(xarray,MaxErrors[1],label = 'Max Error '+methodnames[1])
+sc3 = plt.scatter(xarray,MaxErrors[2],label = 'Max Error '+methodnames[2])
+sc4 = plt.scatter(xarray,MaxErrors[3],label = 'Max Error '+methodnames[3])
+sc5 = plt.scatter(xarray,MaxErrors[4],label = 'Max Error '+methodnames[4])
 
 v2 = plt.vlines(typsize,lowerY,upperY, color = 'g')
 ax = plt.gca()
@@ -211,12 +423,212 @@ ax.set_xscale('log')
 limits = [lowerX, upperX, lowerY, upperY]
 ax.legend()
 x = np.linspace(0.000001,1,100)
+y = x**(1)/10
+y2 = x**2/10
+plt.loglog(x,y,label = 'first order')
+plt.loglog(x,y2, label = 'second order')
+plt.axis(limits)
+plt.title('Max Errors Description'+analyticalCase)
+plt.xlabel(xlabel)
+plt.ylabel(ylabel)
+plt.savefig(path+'/Descr'+analyticalCase+'MaxError'+xlabel+'.png')
+
+plot12 = plt.figure(12)
+sc1 = plt.scatter(xarray,GGErrors[0],label = methodnames[0])
+sc2 = plt.scatter(xarray,GGErrors[1],label = methodnames[1])
+sc3 = plt.scatter(xarray,GGErrors[2],label = methodnames[2])
+sc4 = plt.scatter(xarray,GGErrors[3],label = methodnames[3])
+sc5 = plt.scatter(xarray,GGErrors[4],label = methodnames[4])
+
+v2 = plt.vlines(typsize,lowerY,upperY, color = 'g')
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+limits = [lowerX, upperX, lowerY, upperY]
+ax.legend()
+x = np.linspace(0.000001,1,100)
+y = x**(1)/10
+y2 = x**2/10
+plt.loglog(x,y,label = 'first order')
+plt.loglog(x,y2, label = 'second order')
+plt.axis(limits)
+plt.title('GG Correction Error'+analyticalCase)
+plt.xlabel(xlabel)
+plt.ylabel(ylabel)
+plt.savefig(path+'/Descr'+analyticalCase+'GG correction Error'+xlabel+'.png')
+
+###########################################################################"
+# Nb of cells
+xlabel = 'Number of cells'
+ylabel = 'Relative Magnitude Error'
+xarray = nbs[0]
+upperX = 1e7
+lowerX = 1000
+upperY = 1e-1
+lowerY = 1e-6 
+
+
+plot13 = plt.figure(13)
+v1 = plt.vlines(typsize,lowerY,upperY, color = 'g')
+sc1 = plt.scatter(xarray,MeanErrors[0],label = methodnames[0])
+sc2 = plt.scatter(xarray,MeanErrors[1],label = methodnames[1])
+sc3 = plt.scatter(xarray,MeanErrors[2],label = methodnames[2]) 
+sc4 = plt.scatter(xarray,MeanErrors[3],label = methodnames[3]) 
+sc5 = plt.scatter(xarray,MeanErrors[4],label = methodnames[4]) 
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+limits = [lowerX, upperX, lowerY, upperY]
+ax.legend()
+x = np.linspace(0.000001,1,100)
 y = x**(1)/100
+y2 = x**2/100
 plt.loglog(x,y)
+plt.loglog(x,y2)
+plt.axis(limits)
+plt.title('Mean Errors Description'+analyticalCase+xlabel)
+plt.xlabel(xlabel)
+plt.ylabel(ylabel)
+plt.savefig(path+'/Descr'+analyticalCase+'MeanError'+xlabel+'.png')
+
+plot14 = plt.figure(14)
+sc1 = plt.scatter(xarray,L2Errors[0],label = methodnames[0])
+sc2 = plt.scatter(xarray,L2Errors[1],label = methodnames[1])
+sc3 = plt.scatter(xarray,L2Errors[2],label = methodnames[2])
+sc4 = plt.scatter(xarray,L2Errors[3],label = methodnames[3])
+sc5 = plt.scatter(xarray,L2Errors[4],label = methodnames[4])
+
+v2 = plt.vlines(typsize,lowerY,upperY, color = 'g')
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+limits = [lowerX, upperX, lowerY, upperY]
+ax.legend()
+x = np.linspace(0.000001,1,100)
+y = x**(1)/100
+y2 = x**2/100
+plt.loglog(x,y,label = 'first order')
+plt.loglog(x,y2, label = 'second order')
+plt.axis(limits)
+plt.title('L2 Errors Description'+analyticalCase)
+plt.xlabel(xlabel)
+plt.ylabel(ylabel)
+plt.savefig(path+'/Descr'+analyticalCase+'L2Error'+xlabel+'.png')
+
+
+plot15 = plt.figure(15)
+sc1 = plt.scatter(xarray,MaxErrors[0],label = 'Max Error '+methodnames[0])
+sc2 = plt.scatter(xarray,MaxErrors[1],label = 'Max Error '+methodnames[1])
+sc3 = plt.scatter(xarray,MaxErrors[2],label = 'Max Error '+methodnames[2])
+sc4 = plt.scatter(xarray,MaxErrors[3],label = 'Max Error '+methodnames[3])
+sc5 = plt.scatter(xarray,MaxErrors[4],label = 'Max Error '+methodnames[4])
+
+v2 = plt.vlines(typsize,lowerY,upperY, color = 'g')
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+limits = [lowerX, upperX, lowerY, upperY]
+ax.legend()
+x = np.linspace(0.000001,1,100)
+y = x**(1)/100
+y2 = x**2/100
+plt.loglog(x,y,label = 'first order')
+plt.loglog(x,y2, label = 'second order')
+plt.axis(limits)
+plt.title('Max Errors Description'+analyticalCase)
+plt.xlabel(xlabel)
+plt.ylabel(ylabel)
+plt.savefig(path+'/Descr'+analyticalCase+'MaxError'+xlabel+'.png')
+
+plot16 = plt.figure(16)
+sc1 = plt.scatter(xarray,GGErrors[0],label = methodnames[0])
+sc2 = plt.scatter(xarray,GGErrors[1],label = methodnames[1])
+sc3 = plt.scatter(xarray,GGErrors[2],label = methodnames[2])
+sc4 = plt.scatter(xarray,GGErrors[3],label = methodnames[3])
+sc5 = plt.scatter(xarray,GGErrors[4],label = methodnames[4])
+
+v2 = plt.vlines(typsize,lowerY,upperY, color = 'g')
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+limits = [lowerX, upperX, lowerY, upperY]
+ax.legend()
+x = np.linspace(0.000001,1,100)
+y = x**(1)/100
+y2 = x**2/100
+plt.loglog(x,y,label = 'first order')
+plt.loglog(x,y2, label = 'second order')
 plt.axis(limits)
 plt.title('GG Correction Error'+analyticalCase)
 plt.xlabel(xlabel)
 plt.ylabel(ylabel)
-plt.savefig(path+'/Descr'+analyticalCase+'GG correction Error.png')
+plt.savefig(path+'/Descr'+analyticalCase+'GG correction Error'+xlabel+'.png')
+
+
+
+############################################################################
+# Single Cell
+upperX = 1e0
+lowerX = 5*1e-2
+upperY = 1e-1
+lowerY = 1e-6 
+xlabel = 'Max Connection length'
+ylabel = 'Relative Magnitude Error'
+
+plot = plt.figure(17)
+sc1 = plt.scatter(singledistances,SingleErrors[0],label = methodnames[0])
+sc2 = plt.scatter(singledistances,SingleErrors[1],label = methodnames[1])
+sc3 = plt.scatter(singledistances,SingleErrors[2],label = methodnames[2])
+sc4 = plt.scatter(singledistances,SingleErrors[3],label = methodnames[3])
+sc5 = plt.scatter(singledistances,SingleErrors[4],label = methodnames[4])
+
+v2 = plt.vlines(typsize,lowerY,upperY, color = 'g')
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+limits = [lowerX, upperX, lowerY, upperY]
+ax.legend()
+x = np.linspace(0.000001,1,100)
+y = x**(1)/10
+y2 = x**2/10
+plt.loglog(x,y,label = 'first order')
+plt.loglog(x,y2, label = 'second order')
+plt.axis(limits)
+plt.title('Single Cell Error'+analyticalCase)
+plt.xlabel(xlabel)
+plt.ylabel(ylabel)
+plt.savefig(path+'/Descr'+analyticalCase+'SingleError.png')
+
+
+upperX = 1e-1
+lowerX = 5*1e-3
+xlabel = 'Max Connection length excluding toroidal direction'
+plot = plt.figure(18)
+sc1 = plt.scatter(singledistancesnontor,SingleErrors[0],label = methodnames[0])
+sc2 = plt.scatter(singledistancesnontor,SingleErrors[1],label = methodnames[1])
+sc3 = plt.scatter(singledistancesnontor,SingleErrors[2],label = methodnames[2])
+sc4 = plt.scatter(singledistancesnontor,SingleErrors[3],label = methodnames[3])
+sc5 = plt.scatter(singledistancesnontor,SingleErrors[4],label = methodnames[4])
+
+v2 = plt.vlines(typsize,lowerY,upperY, color = 'g')
+ax = plt.gca()
+ax.set_yscale('log')
+ax.set_xscale('log')
+limits = [lowerX, upperX, lowerY, upperY]
+ax.legend()
+x = np.linspace(0.000001,1,100)
+y = x**(1)/10
+y2 = x**2/10
+plt.loglog(x,y, label='first order')
+plt.loglog(x,y2, label='second order')
+plt.axis(limits)
+plt.title('Single Cell Error Non Tor '+analyticalCase)
+plt.xlabel(xlabel)
+plt.ylabel(ylabel)
+plt.savefig(path+'/Descr'+analyticalCase+'Single Cell Error.png')
+
+
+
+
 
 plt.show()
diff --git a/Scripts_Matthieu/ExtractNoiseInfo.py b/Scripts_Matthieu/ExtractNoiseInfo.py
index adf9b82..9ec01ca 100644
--- a/Scripts_Matthieu/ExtractNoiseInfo.py
+++ b/Scripts_Matthieu/ExtractNoiseInfo.py
@@ -5,17 +5,24 @@ import matplotlib.pyplot as plt
 
 path = os.getcwd()
 
-xnames = ['_1','_2','_3','_4','_5','_7','_8','_9','_10']
+xnames = ['_1','_2','_3','_4','_5','_6','_7','_8','_9']#,'_10']
 Tdata = []
 
 for i in range(len(xnames)):
-    Temp = np.loadtxt('/u/mbj/Drifts/data-for-drift-computations/EIM/N03.00_P05.00_D05_C02.50-TEST/x-out'+xnames[i]+'/x-out/TE_TI', max_rows = 121554 )
-    Tlength = len(Temp)
-#    TempFiltered = Temp[Temp != 0]
-#    TempFiltered = np.nonzero(Temp)
-    TempFiltered = np.reshape(Temp,Tlength*6)
-    print('max',np.amax(TempFiltered))
-    Tdata.append(TempFiltered)
+    file = open('/u/mbj/Drifts/data-for-drift-computations/Reference-FullRes/X-TEST-HR/x-out'+xnames[i]+'/x-out/TE_TI','r')
+    Temptemp = file.readlines()
+    Temp = []
+    for i in range(len(Temptemp)):
+        Temp = Temp+[float(x) for x in Temptemp[i].split()]
+#    print(Temp)
+    print(len(Temp))
+    Temp = np.array(Temp)
+#    Temp = np.loadtxt('/u/mbj/Drifts/data-for-drift-computations/Reference-FullRes/X-TEST-HR/x-out'+xnames[i]+'/x-out/TE_TI', max_rows = 25851)
+#    Temp =np.concatenate(Temp.flatten(),np.loadtxt('/u/mbj/Drifts/data-for-drift-computations/Reference-FullRes/X-TEST-HR/x-out'+xnames[i]+'/x-out/TE_TI',skiprows=25851,max_rows = 1))
+#    Tlength = len(Temp)
+#    print(Tlength)
+#    Temp.flatten()
+    Tdata.append(Temp)
 
 
 avTemp = Tdata[0]    
diff --git a/Scripts_Matthieu/PlotTrajectorySingle.py b/Scripts_Matthieu/PlotTrajectorySingle.py
index 43579f7..26df7bf 100644
--- a/Scripts_Matthieu/PlotTrajectorySingle.py
+++ b/Scripts_Matthieu/PlotTrajectorySingle.py
@@ -4,6 +4,7 @@ import sys
 import AnalysisFunctions
 import math
 import numpy as np
+import plotly.graph_objects as go
 
 path = os.getcwd()
 method = sys.argv[1]
@@ -20,7 +21,7 @@ plt.xlabel('Radial Position [cm]')
 plt.ylabel('Vertical Position [cm]')
 title = 'Trajectory'+ name
 plt.title(title)
-positions = AnalysisFunctions.readText(path+ '/Output' + method+'Descr'+analyticalP+'Noise'+noise+'/TRAJECTORY1')
+positions = AnalysisFunctions.readText(path+ '/Output' + method+'Descr'+analyticalP+'Noise'+noise+'/TRAJECTORY4')
 for i in range(len(positions)):
     if i%2 == 0:
         p = [float(x) for x in positions[i].split()]
@@ -35,13 +36,38 @@ y = [r[2] for r in posNum]
 R = math.sqrt((x[0]-500)**2+y[0]**2)
 print(R)
 
-sc = plt.scatter(x, y, s=10, c = posLen, cmap = 'plasma')
+sc = plt.scatter(x, y, s=50, c = posLen, cmap = 'plasma')
 #circle = plt.Circle((500,0),R, color = 'green')
 angle = np.linspace(0,2*np.pi,150)
 x = (500+R*np.cos(angle))
 y = R*np.sin(angle)
 
 axes.plot(x,y)
-plt.show()
+
 plt.savefig(path+'/Trajectory'+method+analyticalP+noise+name+'.png')
 print(path+method+analyticalP+noise+name+'.png')
+
+
+x = [r[0]*math.cos(r[1]) for r in posNum]
+y = [r[0]*math.sin(r[1]) for r in posNum]
+z = [r[2] for r in posNum]
+
+x = x[0:250]
+y = y[0:250]
+z = z[0:250]
+
+#fig = plt.figure()
+#ax = fig.add_subplot(projection='3d')
+#sc = plt.scatter(x, y, z, cmap = 'plasma')
+
+
+fig = go.Figure(data=[go.Scatter3d(x=x, y=y, z=z,
+                                   mode='markers',
+                                   marker=dict(
+        size=12,
+        color=posLen[0:250],                # set color to an array/list of desired values
+        colorscale='Viridis',   # choose a colorscale
+        opacity=0.8
+    ))])
+fig.show()
+plt.show()
diff --git a/Scripts_Matthieu/PrepareGeometryCompareNodal.py b/Scripts_Matthieu/PrepareGeometryCompareNodal.py
index 80235a1..3a2b94f 100644
--- a/Scripts_Matthieu/PrepareGeometryCompareNodal.py
+++ b/Scripts_Matthieu/PrepareGeometryCompareNodal.py
@@ -4,16 +4,16 @@ 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/CompareNodal75'
+folder = '/u/mbj/Drifts/data-for-drift-computations/CompareNodal105'
 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 = [6, 9,12,18,24]
-nrs = [15,23, 30,45, 60]
-nps = [75,110, 150,225, 300]
+nts = [6, 12,24, 36, 48]
+nrs = [10, 20, 40, 60,80]
+nps = [50, 100, 200, 300,400]
 #cases = [nrs,nps]
 caseNames = ['Case1','Case2','Case3', 'Case4', 'Case5']
 # Number of toroidal, radial, poloidal surfaces
diff --git a/Scripts_Matthieu/TrajectoryDeviation.py b/Scripts_Matthieu/TrajectoryDeviation.py
index 55a025c..9aa7488 100644
--- a/Scripts_Matthieu/TrajectoryDeviation.py
+++ b/Scripts_Matthieu/TrajectoryDeviation.py
@@ -5,7 +5,7 @@ import AnalysisFunctions
 
 analyticalCase = sys.argv[1]
 path = os.getcwd()
-methods = ['6','8']
+methods = ['6','7','8','9', '10']
 noise = '0'
 
 pathlengths = []
@@ -23,7 +23,7 @@ for method in methods:
         t = [float(x) for x in trajInfo[i].split()]
         print(t)
         pls.append(t[0])
-        devs.append(t[1])
+        devs.append(abs(t[1]))
         ts.append(t[2])
         rds.append(devs[i]/pls[i])
     pathlengths.append(pls)
@@ -33,8 +33,8 @@ for method in methods:
 
 upperX = 1e-4
 lowerX = 1e-9
-upperY = 1e-3
-lowerY = 1e-8   
+upperY = 100
+lowerY = 1e-12  
 xlabel = 'Timestep [s]'
 ylabel = 'Relative Deviation (after 1s)'
 
@@ -42,6 +42,9 @@ ylabel = 'Relative Deviation (after 1s)'
 plot1 = plt.figure(1)
 sc1 = plt.scatter(timesteps[0],reldev[0],label = methods[0])
 sc2 = plt.scatter(timesteps[1],reldev[1],label = methods[1])
+sc3 = plt.scatter(timesteps[0],reldev[2],label = methods[2])
+sc4 = plt.scatter(timesteps[1],reldev[3],label = methods[3])
+sc5 = plt.scatter(timesteps[0],reldev[4],label = methods[4])
 
 ax = plt.gca()
 ax.set_yscale('log')
@@ -51,7 +54,7 @@ ax.legend()
 #x = np.linspace(0.000001,1,100)
 #y = x**(1)/100
 #plt.loglog(x,y)
-#plt.axis(limits)
+plt.axis(limits)
 plt.title('Relative Deviation'+analyticalCase)
 plt.xlabel(xlabel)
 plt.ylabel(ylabel)
diff --git a/Scripts_Matthieu/writeInfo.py b/Scripts_Matthieu/writeInfo.py
index bb4fa86..36f50c8 100644
--- a/Scripts_Matthieu/writeInfo.py
+++ b/Scripts_Matthieu/writeInfo.py
@@ -26,10 +26,14 @@ boundarynodes = AnalysisFunctions.readText(path+ '/Output' + method+'Descr'+anal
 gbpositions =  AnalysisFunctions.readText(path+ '/Output' + method+'Descr'+analyticalP+'Noise'+noise+'/GB_POS')
 ggcorr = AnalysisFunctions.readText(path+ '/Output' + method+'Descr'+analyticalP+'Noise'+noise+'/GGCORR')
 
+Rmajor = 500
 if aP ==1:
     edescr = AnalyticalPrescriptions.e_descr1
 elif aP == 5:
     edescr = AnalyticalPrescriptions.e_descr5
+elif aP==6:
+    edescr = AnalyticalPrescriptions.e_descr6
+    Rmajor = 590
 elif aP >= 2:
     edescr = AnalyticalPrescriptions.e_descr2
     
@@ -38,6 +42,6 @@ gdescr = AnalyticalPrescriptions.gradb_description1
 exb = driftFunctions.ExBdrift
 gB = driftFunctions.gradBdrift
 
-allInfo = AnalysisFunctions.getAllInfo(positions,gbpositions,edrifts,efield,gradb,gdrifts,boundarycells,boundarynodes,edescr,bdescr,gdescr,exb, gB, ggcorr, method+'Descr'+analyticalP+'Noise'+noise)
+allInfo = AnalysisFunctions.getAllInfo(positions,gbpositions,edrifts,efield,gradb,gdrifts,boundarycells,boundarynodes,edescr,bdescr,gdescr,exb, gB, ggcorr,Rmajor, method+'Descr'+analyticalP+'Noise'+noise)
 
 print('CSV file written')
-- 
GitLab