diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..84fce74
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,9 @@
+*.csv
+*.mod
+*.o
+*.bkp
+binario
+result
+binario_gradiente
+binario_gradiente1
+pequeno
diff --git a/csv2vtk_particles.py b/csv2vtk_particles.py
index 8c67503..f3d8966 100755
--- a/csv2vtk_particles.py
+++ b/csv2vtk_particles.py
@@ -13,6 +13,14 @@
 import numpy as np
 import csv, os, shutil
 import configparser
+try:
+    import progressbar
+    pbar = True 
+except:
+    print('progressbar not available. Try one:')
+    print('conda install -c conda-forge progressbar2')
+    print('conda install -c conda-forge/label/gcc7 progressbar2')
+    pbar = False
 
 via = os.getcwd()
 
@@ -42,7 +50,7 @@
 
 config = configparser.ConfigParser()
 config.read('settings.ini')
-
+print('Reading settings.ini')
 N = int(config['global']['N'].split()[0])
 nimpre =  int(config['global']['nimpre'].split()[0])
 ntype = int(config['global']['Ntype'].split()[0])
@@ -95,6 +103,9 @@
 
 nID = np.linspace(1,N,N)
 
+print('Converting...')
+if pbar:
+    bar = progressbar.ProgressBar(max_value=nimpre)
 for fnum in range(0,nimpre+1):
     with open('temp/position.csv.'+str(fnum),encoding='utf-8') as file_locus:
         csv_lector = csv.reader(file_locus,delimiter = ',')
@@ -131,7 +142,8 @@
         pointsToVTK(via+'/'+folder +'/grupo'+ str(grupo) + '_' +str(fnum), xs, ys, zs, data = {"Vx" : vxs, "Vy" : vys, "Tipo" : tipos, "raio_solido" : rsols, "nID" : nIDs })       
         grupo += 1
 
-
+    if pbar:
+        bar.update(fnum)
 shutil.rmtree('temp')
         
         
@@ -144,6 +156,18 @@
 except:
     print('No settings file found!\n')
 
+include_folder = True 
+if os.path.isfile('.gitignore'):
+    with open('.gitignore','a+') as file:
+        for line in file:
+            if line == folder:
+                include_folder = False
+        if include_folder:
+            file.write(folder+'\n')
+else:
+    with open('.gitignore','a+') as file:
+        file.write(folder+'\n')
+
     #with open('cell.csv.'+str(fnum),encoding='utf-8') as file_velocitas:
         #csv_lector = csv.reader(file_velocitas,delimiter = ',')
         #i = 0
diff --git a/settings.ini b/settings.ini
index 7780d21..eac851b 100755
--- a/settings.ini
+++ b/settings.ini
@@ -1,6 +1,6 @@
 # PARAMETROS GLOBAIS
 [global]
-    nimpre  = 200 # quntidade de saidas 
+    nimpre  = 5000 # quntidade de saidas 
     N = 1089 # número de partículas 
     Ntype = 1 # número de tipos de partícula
     t_fim = 50  # tempo de execução 
diff --git a/viscosity.py b/viscosity.py
index ba1c9fa..b7fa108 100644
--- a/viscosity.py
+++ b/viscosity.py
@@ -5,7 +5,14 @@
 import matplotlib.pyplot as plt
 from glob import glob
 import time
-import curses # para progressbar
+try:
+    import progressbar
+    pbar = True 
+except:
+    print('progressbar not available. Try one:')
+    print('conda install -c conda-forge progressbar2')
+    print('conda install -c conda-forge/label/gcc7 progressbar2')
+    pbar = False
 
 def integral(F,i,j,dt):
     # sh = np.shape(F)
@@ -24,30 +31,37 @@ def integral(F,i,j,dt):
                 #     )
     return I
 
-def progressbar(step,total,message, stdscr):
-    # usar import curses
-    # inicializar uma strig vazia fora do loop para stdcr
-    if step == 0:  
-        stdscr = curses.initscr() 
-        stdscr.addstr(0, 0, message + "|" + " "*20 + "| {} %".format(0))
-        stdscr.refresh()
-        return stdscr
-    elif step < total-1:
-        a = int(20*step/total)
-        b = 20-a
-        stdscr.addstr(0, 0, message + "|" + "#"*a + " "*b + "| {} %".format(5*a))
-        stdscr.refresh()
-        return stdscr
-        #print("#",end='')
-    else:
-        a = int(20*step/total)
-        b = 20-a        
-        stdscr.addstr(0, 0, message + "|" + "#"*a + " "*b + "| {} %".format(5*a))
-        curses.endwin()
-        print("\n")    
-        return "0"
-    
+def diff1(F,h):
+    B = np.zeros((len(F),1))
+    for i in range(len(F)):
+        if i > 1 and i < len(F) - 2:
+            B[i] = ((1/12)*F[i-2] + (-2/3)*F[i-1] + (2/3)*F[i+1] + (-1/12)*F[i+2])/h   
+        elif i == 0:
+            B[i] = ((-11/6)*F[0] + 3*F[1] + (-3/2)*F[2] + (1/3)*F[3])/h 
+        elif i == 1:
+            B[i] = (-2*F[0]-3*F[1]+6*F[2]-1*F[3])/(6*h)
+        elif i == len(F)-1:
+            B[i] = (1*F[i-2]-4*F[i-1]+3*F[i+0])/(2*1.0*h**1)
+        elif i == len(F)-2:
+            B[i] = (1*F[i-2]-6*F[i-1]+3*F[i+0]+2*F[i+1])/(6*1.0*h**1)
+    return B
+
+def greenkubo(Tk,Tp,i,j,na,dt):
+    itv = int(Tp.shape[2]/na)
+    kk = np.zeros((itv))
+    kp = np.zeros((itv))
+    pp = np.zeros((itv))
+    for a in range(itv):
+        for b in range(na):
+            kk[a] += np.trapz(Tk[i,j,b*itv:b*itv+a],dx=dt)**2
+            kp[a] += np.trapz(Tp[i,j,b*itv:b*itv+a],dx=dt)*np.trapz(Tk[i,j,b*itv:b*itv+a],dx=dt)
+            pp[a] += np.trapz(Tp[i,j,b*itv:b*itv+a],dx=dt)**2
+    kk = kk/na
+    pp = pp/na
+    kp = kp/na 
+    return kk, kp, pp
 
+    
 dirname = os.getcwd() #os.path.dirname(os.path.abspath(__file__))
 dirlist = glob(dirname + "/*/")
 print("Choose a folder there the results are contained:\nNo | Folder")
@@ -64,7 +78,7 @@ def progressbar(step,total,message, stdscr):
 n_files = int(config['out_files']['out_files'].split()[0])
 ntype = int(config['global']['Ntype'].split()[0])
 t_fim = float(config['global']['t_fim'].split()[0])
-dt = float(config['global']['dt'].split()[0])
+#dt = float(config['global']['dt'].split()[0])
 F = np.array([0,0])
 quant = []
 sigma = []
@@ -110,12 +124,12 @@ def progressbar(step,total,message, stdscr):
 
 step = 0 
 n1,n2 = 0,0
-stdscr = 's' #para barra de progresso
 
+if pbar:
+    bar = progressbar.ProgressBar(max_value=nsteps)
 while step <= nsteps: 
     particle_map = [[[] for _ in range(mesh[1])] for _ in range(mesh[0])]
-    #stdscr = progressbar(step,nsteps,'Computing:',stdscr)
-    print(step)
+
     pos = pd.read_csv(res_dir+"position.csv."+str(step), header=None, names = ["x","y"])
     vel = pd.read_csv(res_dir+"velocity.csv."+str(step), header=None, names = ["v_x", "v_y"])
     n = [x for x in range(len(pos))]
@@ -165,6 +179,7 @@ def progressbar(step,total,message, stdscr):
                     tauxyp[i,j,step] += -F[1]*(x1[0] - x2[0])  # Falta multiplicar por 1/Volume
                 # na celula vizinha i+1, j 
                 if i+1 < region[1]:
+                    print('A')
                     for nn2 in range(1,len(particle_map[i+1][j])):
                         n2 = particle_map[i+1][j][nn2]
                         rs2 = rs[pos_vel.loc[n2,'tipo']]
@@ -184,6 +199,7 @@ def progressbar(step,total,message, stdscr):
                             time.sleep(2)
                 # na celula vizinha i+1, j+1 
                 if i+1 < region[1] and j+1 < region[3]:
+                    print('B')
                     for nn2 in range(1,len(particle_map[i+1][j+1])):
                         n2 = particle_map[i+1][j+1][nn2]
                         rs2 = rs[pos_vel.loc[n2,'tipo']]
@@ -203,6 +219,7 @@ def progressbar(step,total,message, stdscr):
                             time.sleep(2)
                 # na celula vizinha i, j+1 
                 if j+1 < region[3]:
+                    print('C')
                     for nn2 in range(1,len(particle_map[i][j+1])):
                         n2 = particle_map[i][j+1][nn2]
                         rs2 = rs[pos_vel.loc[n2,'tipo']]
@@ -222,6 +239,7 @@ def progressbar(step,total,message, stdscr):
                             time.sleep(2)
                 # na celula vizinha i-1, j+1 
                 if i-1 >= region[0] and j+1 < region[3]:
+                    print('D')
                     for nn2 in range(1,len(particle_map[i-1][j+1])):
                         n2 = particle_map[i-1][j+1][nn2]
                         rs2 = rs[pos_vel.loc[n2,'tipo']]
@@ -239,6 +257,8 @@ def progressbar(step,total,message, stdscr):
                         if r == 0:
                             print('nn {} {}'.format(nn,nn2))
                             time.sleep(2)
+    if pbar:
+        bar.update(step)
     step += 1
 
 
@@ -246,25 +266,29 @@ def progressbar(step,total,message, stdscr):
 tauxyp = tauxyp/Vol
 tauxyk = tauxyk/Vol
 
-Ik = np.zeros((mesh[0],mesh[1],nsteps))
-Ip = np.zeros((mesh[0],mesh[1],nsteps))
-etakk = np.zeros((mesh[0],mesh[1],nsteps))
-etakp = np.zeros((mesh[0],mesh[1],nsteps))
-etapp = np.zeros((mesh[0],mesh[1],nsteps))
-
 # vamos fazer um gráfico 
-for step in range(nsteps):
-    for i in range(region[0], region[1]):
-        for j in range(region[2], region[3]):
-            Ik[i,j,step] = integral(tauxyk,i,j,dt)
-            Ip[i,j,step] = integral(tauxyp,i,j,dt)
-            
-            # tem que dividir por 2kT
-            etakk[i,j,step] = Vol*Ik[i,j,step]**2/density_map[i,j,step]
-            etakp[i,j,step] = Vol*Ik[i,j,step]*Ip[i,j,step]/density_map[i,j,step]
-            etapp[i,j,step] = Vol*Ip[i,j,step]**2/density_map[i,j,step]
 
-plt.plot(etakk[0,0,:]+etakp[0,0,:]+etapp[0,0,:])
+na = int(input('Enter de number of assembles (na) to average.\nThe number of steps for correlation will be nsteps/na: '))
+etakk = np.zeros((mesh[0],mesh[1],int(nsteps/na)))
+etakp = np.zeros((mesh[0],mesh[1],int(nsteps/na)))
+etapp = np.zeros((mesh[0],mesh[1],int(nsteps/na)))
+for i in range(region[0], region[1]):
+    for j in range(region[2], region[3]):
+        etakk[i,j,:],etakp[i,j,:],etapp[i,j,:] = greenkubo(tauxyk,tauxyp,i,j,na,dt)
+
+ij = input('i j pro grafico: ').split()
+i,j = int(ij[0]),int(ij[1])
+t = np.linspace(0,t_fim/na,int(nsteps/na))
+plt.figure(1)
+plt.plot(t,etapp[i,j,:],label='etapp*t')
+plt.plot(t,etakp[i,j,:],label='etakp*t')
+plt.plot(t,etakk[i,j,:],label='etakk*t')
+etatot = etakk[i,j,:]+etakp[i,j,:]+etapp[i,j,:]
+plt.plot(t,etatot,'.k',label='sum')
+m,b = np.polyfit(t, etatot, 1)
+y = m*t + b 
+plt.plot(t,y,'-k',label='Eta*t fit')
+plt.legend()
 plt.show()
 
 
diff --git a/viscosity_allen11.py b/viscosity_allen11.py
new file mode 100644
index 0000000..80bc622
--- /dev/null
+++ b/viscosity_allen11.py
@@ -0,0 +1,187 @@
+#viscosidade calculada com o livro do Allen usando as partes do Q lá
+
+import numpy as np 
+import pandas as pd 
+import os
+import configparser
+import matplotlib.pyplot as plt
+from glob import glob
+import time
+try:
+    import progressbar
+    pbar = True 
+except:
+    print('progressbar not available. Try one:')
+    print('conda install -c conda-forge progressbar2')
+    print('conda install -c conda-forge/label/gcc7 progressbar2')
+    pbar = False
+
+def diff1(F,h):
+    B = np.zeros((len(F),1))
+    for i in range(len(F)):
+        if i > 1 and i < len(F) - 2:
+            B[i] = ((1/12)*F[i-2] + (-2/3)*F[i-1] + (2/3)*F[i+1] + (-1/12)*F[i+2])/h   
+        elif i == 0:
+            B[i] = ((-11/6)*F[0] + 3*F[1] + (-3/2)*F[2] + (1/3)*F[3])/h 
+        elif i == 1:
+            B[i] = (-2*F[0]-3*F[1]+6*F[2]-1*F[3])/(6*h)
+        elif i == len(F)-1:
+            B[i] = (1*F[i-2]-4*F[i-1]+3*F[i+0])/(2*1.0*h**1)
+        elif i == len(F)-2:
+            B[i] = (1*F[i-2]-6*F[i-1]+3*F[i+0]+2*F[i+1])/(6*1.0*h**1)
+    return B
+
+def einstein_relation(A,na):
+    itv = int(len(A)/na)
+    B = np.zeros((itv,1))
+    for j in range(itv):
+        for i in range(na):
+            B[j] += (A[i*itv+j] - A[i*itv])**2
+    B = B/na 
+    return B
+
+dirname = os.getcwd() #os.path.dirname(os.path.abspath(__file__))
+dirlist = glob(dirname + "/*/")
+print("Choose a folder there the results are contained:\nNo | Folder")
+for a in range(len(dirlist)):
+    print("{} | {}\n".format(a,dirlist[a]))
+a = int(input("Enter the number of the folder\n"))
+res_dir = dirlist[a]
+
+config = configparser.ConfigParser()
+config.read(res_dir + 'settings.txt')
+N = int(config['global']['N'].split()[0])
+dimx = float(config['global']['dimX'].split()[0])
+dimy = float(config['global']['dimY'].split()[0])
+n_files = int(config['out_files']['out_files'].split()[0])
+
+ntype = int(config['global']['Ntype'].split()[0])
+t_fim = float(config['global']['t_fim'].split()[0])
+dt = float(config['global']['dt'].split()[0])
+F = np.array([0,0])
+quant = []
+sigma = []
+epsilon = []
+rs = [] # raio sólido
+mass = []
+tipo = [0]*N
+
+dt = t_fim/n_files
+
+a = input("Enter the subdomains mesh dimensions.\n")
+a = a.split()
+mesh = np.array([int(a[0]), int(a[1])])
+a = input("Enter a location. xmin xmax ymin ymax\n")
+if a == '':
+    region = [0, mesh[0], 0, mesh[1]]
+else:
+    region = [int(x) for x in a.split()]
+
+Vol = (dimx/mesh[0])*(dimy/mesh[1]) # volume dos elementos da malha
+for i in range(ntype):
+    quant.append(int(config['par_'+str(i)]['quantidade'].split()[0]))
+    rs.append(float(config['par_'+str(i)]['rs'].split()[0]))
+    sigma.append(float(config['par_'+str(i)]['sigma'].split()[0]))
+    epsilon.append(float(config['par_'+str(i)]['epsilon'].split()[0]))
+    mass.append(float(config['par_'+str(i)]['m'].split()[0]))
+j,k = 0,0
+for i in range(len(quant)):
+    for j in range(quant[i]):
+        tipo[j+k] = i
+        
+    k = sum(quant[0:i+1])
+tipo = pd.DataFrame(tipo, columns=["tipo"]) # numero id da partícula
+
+hx = dimx/mesh[0]
+hy = dimy/mesh[1]
+
+nsteps = n_files # int(input('Enter the number of steps:\n'))
+
+density_map = np.zeros((mesh[0],mesh[1], nsteps+1))
+r = np.zeros((mesh[0],mesh[1],nsteps+1))
+dQm = np.zeros(nsteps)
+eta = np.zeros(nsteps)
+
+
+step = 0 
+n1,n2 = 0,0
+stdscr = 's' #para barra de progresso
+sample_list = []
+r0 = []
+p0 = []
+Qxy = np.zeros((nsteps+1,1))
+Qyx = np.zeros((nsteps+1,1))
+#eta = np.zeros((nsteps,1))
+
+if pbar:
+    bar = progressbar.ProgressBar(max_value=nsteps)
+while step <= nsteps: 
+    particle_map = [[[] for _ in range(mesh[1])] for _ in range(mesh[0])]
+    pos = pd.read_csv(res_dir+"position.csv."+str(step), header=None, names = ["x","y"])
+    vel = pd.read_csv(res_dir+"velocity.csv."+str(step), header=None, names = ["v_x", "v_y"])
+    n = [x for x in range(len(pos))]
+    n = pd.DataFrame(n, columns=["n"]) # numero id da partícula
+    pos_vel = pd.concat([n,pos,vel,tipo],axis=1)
+   
+    for nn in range(len(pos_vel)):
+        xp = int(pos_vel.loc[nn,'x']//hx)
+        yp = int(pos_vel.loc[nn,'y']//hy)
+        if xp == mesh[0]:
+            xp = xp - 1
+        if yp == mesh[1]:
+            yp = yp - 1
+        particle_map[xp][yp].append( pos_vel.loc[nn,'n'] )
+        density_map[xp,yp,step] += 1
+
+    if (step == 0):
+        for i in range(region[0],region[1]):
+            for j in range(region[2],region[3]):
+                for nn in range(len(particle_map[i][j])):
+                    sample_list.append(particle_map[i][j][nn])
+                    n1 = particle_map[i][j][nn]
+                    r0.append([pos_vel.loc[n1,'x'], pos_vel.loc[n1,'y']])
+                    m = mass[pos_vel.loc[n1,'tipo']]
+                    p0.append([pos_vel.loc[n1,'v_x']*m, pos_vel.loc[n1,'v_y']*m])
+        r0 = np.array(r0)
+        r1 = np.zeros( (len(sample_list),2) )
+        p0 = np.array(p0)
+        p1 = np.zeros( (len(sample_list),2) )
+        Qxy[0] = np.sum(r0[:,0]*p0[:,1])/Vol
+        Qyx[0] = np.sum(r0[:,1]*p0[:,0])/Vol
+    else:
+        for nn in range(len(sample_list)):
+            n1 = sample_list[nn]
+            r1[nn,:] =  [pos_vel.loc[n1,'x'], pos_vel.loc[n1,'y']]
+            m = mass[pos_vel.loc[n1,'tipo']]
+            p1[nn,:] = [pos_vel.loc[n1,'v_x']*m, pos_vel.loc[n1,'v_y']*m]
+
+        Qxy[step] = np.sum(r1[:,0]*p1[:,1], axis=0)/Vol
+        Qyx[step] = np.sum(r1[:,1]*p1[:,0], axis=0)/Vol
+
+    if pbar:
+        bar.update(step)
+    step += 1
+
+
+na = int(input('Enter de number of assembles (na) to average.\nThe number of steps for correlation will be nsteps/na: '))
+
+etaxy2t = einstein_relation(Qxy,na)
+etayx2t = einstein_relation(Qyx,na)
+
+plt.figure(1)
+plt.plot(etaxy2t, label='2t*etaxy')
+plt.plot(etayx2t, label='2t*etayx')
+plt.plot((etayx2t+etaxy2t)/2,'--k', label='2*t*eta mean')
+plt.legend()
+
+
+h = t_fim/nsteps 
+etaxy = diff1(etaxy2t,h)/2
+etayx = diff1(etayx2t,h)/2
+plt.figure(2)
+plt.plot(etaxy[1:len(etaxy)], label='etaxy*kT')
+plt.plot(etayx[1:len(etaxy)], label='etayx*kT')
+plt.plot((etayx[1:len(etaxy)]+etaxy[1:len(etaxy)])/2,'--k', label='eta*kT mean')
+plt.legend()
+
+plt.show()
\ No newline at end of file
diff --git a/viscosity_allen21.py b/viscosity_allen21.py
new file mode 100644
index 0000000..ecf9226
--- /dev/null
+++ b/viscosity_allen21.py
@@ -0,0 +1,182 @@
+import numpy as np 
+import pandas as pd 
+import os
+import configparser
+import matplotlib.pyplot as plt
+from glob import glob
+import time
+try:
+    import progressbar
+    pbar = True 
+except:
+    print('progressbar not available. Try one:')
+    print('conda install -c conda-forge progressbar2')
+    print('conda install -c conda-forge/label/gcc7 progressbar2')
+    pbar = False
+
+
+def diff1(F,h):
+    B = np.zeros((len(F),1))
+    for i in range(len(F)):
+        if i > 1 and i < len(F) - 2:
+            B[i] = ((1/12)*F[i-2] + (-2/3)*F[i-1] + (2/3)*F[i+1] + (-1/12)*F[i+2])/h   
+        elif i == 0:
+            B[i] = ((-11/6)*F[0] + 3*F[1] + (-3/2)*F[2] + (1/3)*F[3])/h 
+        elif i == 1:
+            B[i] = (-2*F[0]-3*F[1]+6*F[2]-1*F[3])/(6*h)
+        elif i == len(F)-1:
+            B[i] = (1*F[i-2]-4*F[i-1]+3*F[i+0])/(2*1.0*h**1)
+        elif i == len(F)-2:
+            B[i] = (1*F[i-2]-6*F[i-1]+3*F[i+0]+2*F[i+1])/(6*1.0*h**1)
+    return B
+
+def einstein_relation(A,na):
+    itv = int(len(A)/na)
+    B = np.zeros((itv,1))
+    for j in range(itv):
+        for i in range(na):
+            B[j] += (A[i*itv+j] - A[i*itv])**2
+    B = B/na 
+    return B
+    
+    
+
+dirname = os.getcwd() #os.path.dirname(os.path.abspath(__file__))
+dirlist = glob(dirname + "/*/")
+print("Choose a folder there the results are contained:\nNo | Folder")
+for a in range(len(dirlist)):
+    print("{} | {}\n".format(a,dirlist[a]))
+a = int(input("Enter the number of the folder\n"))
+res_dir = dirlist[a]
+
+config = configparser.ConfigParser()
+config.read(res_dir + 'settings.txt')
+N = int(config['global']['N'].split()[0])
+dimx = float(config['global']['dimX'].split()[0])
+dimy = float(config['global']['dimY'].split()[0])
+n_files = int(config['out_files']['out_files'].split()[0])
+
+ntype = int(config['global']['Ntype'].split()[0])
+t_fim = float(config['global']['t_fim'].split()[0])
+dt = float(config['global']['dt'].split()[0])
+F = np.array([0,0])
+quant = []
+sigma = []
+epsilon = []
+rs = [] # raio sólido
+mass = []
+tipo = [0]*N
+
+dt = t_fim/n_files
+
+a = input("Enter the subdomains mesh dimensions.\n")
+a = a.split()
+mesh = np.array([int(a[0]), int(a[1])])
+a = input("Enter a location. xmin xmax ymin ymax\n")
+if a == '':
+    region = [0, mesh[0], 0, mesh[1]]
+else:
+    region = [int(x) for x in a.split()]
+
+Vol = (dimx/mesh[0])*(dimy/mesh[1]) # volume dos elementos da malha
+for i in range(ntype):
+    quant.append(int(config['par_'+str(i)]['quantidade'].split()[0]))
+    rs.append(float(config['par_'+str(i)]['rs'].split()[0]))
+    sigma.append(float(config['par_'+str(i)]['sigma'].split()[0]))
+    epsilon.append(float(config['par_'+str(i)]['epsilon'].split()[0]))
+    mass.append(float(config['par_'+str(i)]['m'].split()[0]))
+j,k = 0,0
+for i in range(len(quant)):
+    for j in range(quant[i]):
+        tipo[j+k] = i
+        
+    k = sum(quant[0:i+1])
+tipo = pd.DataFrame(tipo, columns=["tipo"]) # numero id da partícula
+
+hx = dimx/mesh[0]
+hy = dimy/mesh[1]
+
+nsteps = n_files # int(input('Enter the number of steps:\n'))
+
+density_map = np.zeros((mesh[0],mesh[1], nsteps+1))
+r = np.zeros((mesh[0],mesh[1],nsteps+1))
+dQm = np.zeros(nsteps)
+eta = np.zeros(nsteps)
+
+
+step = 0 
+n1,n2 = 0,0
+stdscr = 's' #para barra de progresso
+sample_list = []
+Qxy = np.zeros((nsteps+1,1))
+Qyx = np.zeros((nsteps+1,1))
+Npart = np.zeros((nsteps+1,1))
+#eta = np.zeros((nsteps,1))
+
+if pbar:
+    bar = progressbar.ProgressBar(max_value=nsteps)
+while step <= nsteps: 
+    r0 = []
+    p0 = []
+
+    particle_map = [[[] for _ in range(mesh[1])] for _ in range(mesh[0])]
+    pos = pd.read_csv(res_dir+"position.csv."+str(step), header=None, names = ["x","y"])
+    vel = pd.read_csv(res_dir+"velocity.csv."+str(step), header=None, names = ["v_x", "v_y"])
+    n = [x for x in range(len(pos))]
+    n = pd.DataFrame(n, columns=["n"]) # numero id da partícula
+    pos_vel = pd.concat([n,pos,vel,tipo],axis=1)
+   
+    for nn in range(len(pos_vel)):
+        xp = int(pos_vel.loc[nn,'x']//hx)
+        yp = int(pos_vel.loc[nn,'y']//hy)
+        if xp == mesh[0]:
+            xp = xp - 1
+        if yp == mesh[1]:
+            yp = yp - 1
+        particle_map[xp][yp].append( pos_vel.loc[nn,'n'] )
+        density_map[xp,yp,step] += 1
+
+
+    for i in range(region[0],region[1]):
+        for j in range(region[2],region[3]):
+            for nn in range(len(particle_map[i][j])):
+                Npart[step] += 1 #numero de partículas na região neste passo 
+                n1 = particle_map[i][j][nn]
+                r0.append([pos_vel.loc[n1,'x'], pos_vel.loc[n1,'y']])
+                m = mass[pos_vel.loc[n1,'tipo']]
+                p0.append([pos_vel.loc[n1,'v_x']*m, pos_vel.loc[n1,'v_y']*m])
+    r0 = np.array(r0)
+    p0 = np.array(p0)
+    Qxy[step] += np.sum(r0[:,0]*p0[:,1], axis=0)/Vol 
+    Qyx[step] += np.sum(r0[:,1]*p0[:,0], axis=0)/Vol
+    if pbar:
+        bar.update(step)
+    step += 1
+
+na = int(input('Enter de number of assembles (na) to average.\nThe number of steps for correlation will be nsteps/na: '))
+
+etaxy2t = einstein_relation(Qxy,na)
+etayx2t = einstein_relation(Qyx,na)
+
+t = np.linspace(0,t_fim/na,len(etaxy2t))
+plt.figure(1)
+plt.plot(t,etaxy2t,'.r', label='etaxy')
+m,b = np.polyfit(t, etaxy2t, 1)
+y = m*t + b 
+plt.plot(t,y,'-r', label='etaxy fit')
+plt.plot(t,etayx2t,'.b', label='etayx')
+m,b = np.polyfit(t, etayx2t, 1)
+y = m*t + b 
+plt.plot(t,y,'-b', label='etaxy fit')
+plt.legend()
+
+h = t_fim/nsteps 
+etaxy = diff1(etaxy2t,h)/2
+etayx = diff1(etayx2t,h)/2
+plt.figure(2)
+plt.plot(etaxy[1:len(etaxy)], label='etaxy*kT')
+plt.plot(etayx[1:len(etaxy)], label='etayx*kT')
+plt.plot((etayx[1:len(etaxy)]+etaxy[1:len(etaxy)])/2,'--k', label='eta*kT mean')
+plt.legend()
+
+plt.show()
\ No newline at end of file