Skip to content

Commit

Permalink
algumas mudanças + distribution
Browse files Browse the repository at this point in the history
  • Loading branch information
g7fernandes committed Sep 8, 2019
1 parent e2b75d5 commit eb87dca
Show file tree
Hide file tree
Showing 3 changed files with 137 additions and 29 deletions.
93 changes: 93 additions & 0 deletions distribution_particle.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,93 @@
import numpy as np
import vtk
from vtk.util.numpy_support import vtk_to_numpy
import configparser

def trapz(x,y):
res = 0
for i in range(1,len(x)):
res += (y[i-1]+y[i])*(x[i-1]+x[i])/2
return res

# LÊ A CONFIGURAÇÃO

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])
nimpre_init = int(config['global']['nimpre_init'].split()[0])
dimx = int(config['global']['dimX'].split()[0])
dimy = int(config['global']['dimY'].split()[0])

quant = []


for i in range(ntype):
quant.append(int(config['par_'+str(i)]['quantidade'].split()[0]))

print("Initial step: {}. Final step: {}. Number of particle groups (ntype): {}\n".format(nimpre_init,nimpre,ntype))

pgroup = input("Enter the particle group (>= 0, < ntype): ")
stepini = input("Enter the initial step: ")
stepfim = input("Enter the final step: ")

x = np.zeros((stepfim - stepini, quant[int(pgroup)], 2)

# LÊ OS VTK

grupo = "./result/grupo" + pgroup + "_"

reader = vtk.vtkXMLUnstructuredGridReader()
nfiles = 0
for file in range(stepini,stepfim)
reader.SetFileName(grupo+str(file))
reader.Update()
data = reader.GetOutput()

points = data.GetPoints()
x[nfiles,:,:] = vtk_to_numpy(points.GetData())[:,0:2]# coluna 1 = x, coluna 2 = y
nfiles += 1
#u = vtk_to_numpy(data.GetPointData().GetArray(0)) # colunas: Vx, Vy, Vz

# Todas as posições carregadas

aux1 = input("Enter the number of subdivisions and the direction (x or y).\nIf direction not given, will be x.\n")).split()
direc = 0
if len(aux1) > 1:
if aux[0] == 'y':
direc = 1


ndensiy = np.zeros((nfiles,aux1))
dhist = dimx/int(aux1[0]) # distancia no histograma

for file in range(nfiles):
for i in range(quant[int(pgroup)]):
ndensiy[file, (x[file,i,direc]//dhist) ] += 1


n2density = (np.sum(ndensity,axis=0)/nfiles)/(quant[int(pgroup)])
xx = np.linspace(0,dimx, int(aux1[0]))
plt.plot(n2density)

du = quant[int(pgroup)]/dimx #distribuição uniforme


DS = trapz(x,n2density*np.log(n2density/du))





#N = reader.GetNumberOfPoints() # numero de pontos
#npa = reader.GetNumberOfPointArrays() # numero de vetores de ponto

# n_arrays = reader.GetNumberOfPointArrays()
# print("Available data:\n")
# for i in range(n_arrays):
# print(reader.GetPointArrayName(i))



35 changes: 19 additions & 16 deletions lennard.f90
Original file line number Diff line number Diff line change
Expand Up @@ -845,7 +845,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,
ptr%p%F = aux2 + ptr%p%F +fR
ptrn%p%F = -aux2 + ptrn%p%F - fR
else if (rs1 > 0 .and. rs2 == 0) then
gamma = theta(-pa +ptr%p%n) + atan(coss/sine)
gamma = theta(-pa +ptr%p%n) + atan2(coss,sine)

aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * &
( (sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * &
Expand All @@ -869,7 +869,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,
! partlst(ptr%p%n - pa) = 1

else if (rs2 > 0 .and. rs1 == 0) then
gamma = theta(-pa +ptrn%p%n) + atan(coss/sine)
gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine)
aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ((sigma*(1+B*sin(beta*gamma))/r)**6 * &
(1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * &
[(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
Expand Down Expand Up @@ -956,7 +956,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,
ptrn%p%F = -aux2 + ptrn%p%F - fR
else if (rs1 > 0 .and. rs2 == 0) then

gamma = theta(-pa +ptr%p%n) + atan(coss/sine)
gamma = theta(-pa +ptr%p%n) + atan2(coss,sine)
aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 * &
(1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * &
[(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
Expand All @@ -971,7 +971,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,
ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss]
! partlst(ptr%p%n - pa) = 1
else if (rs2 > 0 .and. rs1 == 0) then
gamma = theta(-pa +ptrn%p%n) + atan(coss/sine)
gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine)
aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 &
* (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * &
[(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
Expand Down Expand Up @@ -1050,7 +1050,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,
ptrn%p%F = -aux2 + ptrn%p%F - fR
else if (rs1 > 0 .and. rs2 == 0) then

gamma = theta(-pa +ptr%p%n) + atan(coss/sine)
gamma = theta(-pa +ptr%p%n) + atan2(coss,sine)
aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ((sigma*(1+B*sin(beta*gamma))/r)**6 &
* (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * &
[(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
Expand All @@ -1066,7 +1066,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,
ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss]
! partlst(ptr%p%n - pa) = 1
else if (rs2 > 0 .and. rs1 == 0) then
gamma = theta(-pa +ptrn%p%n) + atan(coss/sine)
gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine)
aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * &
( (sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * &
[(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
Expand Down Expand Up @@ -1149,7 +1149,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,
ptrn%p%F = -aux2 + ptrn%p%F - fR
else if (rs1 > 0 .and. rs2 == 0) then

gamma = theta(-pa +ptr%p%n) + atan(coss/sine)
gamma = theta(-pa +ptr%p%n) + atan2(coss,sine)
aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * &
( (sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) &
* [(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
Expand All @@ -1163,7 +1163,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,
ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss]
! partlst(ptr%p%n - pa) = 1
else if (rs2 > 0 .and. rs1 == 0) then
gamma = theta(-pa +ptrn%p%n) + atan(coss/sine)
gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine)
aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 * &
(1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * &
[(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
Expand Down Expand Up @@ -1240,7 +1240,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,

else if (rs1 > 0 .and. rs2 == 0) then

gamma = theta(-pa +ptr%p%n) + atan(coss/sine)
gamma = theta(-pa +ptr%p%n) + atan2(coss,sine)
aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * &
((sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * &
[(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
Expand All @@ -1254,7 +1254,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,
ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss]
! partlst(ptr%p%n - pa) = 1
else if (rs2 > 0 .and. rs1 == 0) then
gamma = theta(-pa +ptrn%p%n) + atan(coss/sine)
gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine)
aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 * &
(1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * &
[(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
Expand Down Expand Up @@ -1331,7 +1331,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,
ptrn%p%F = -aux2 + ptrn%p%F - fR
else if (rs1 > 0 .and. rs2 == 0) then

gamma = theta(-pa +ptr%p%n) + atan(coss/sine)
gamma = theta(-pa +ptr%p%n) + atan2(coss,sine)
aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 * &
(1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * &
[(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
Expand All @@ -1345,7 +1345,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,
ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss]
! partlst(ptr%p%n - pa) = 1
else if (rs2 > 0 .and. rs1 == 0) then
gamma = theta(-pa +ptrn%p%n) + atan(coss/sine)
gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine)
aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * &
( (sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * &
[(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
Expand Down Expand Up @@ -1423,7 +1423,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,
ptrn%p%F = -aux2 + ptrn%p%F - fR
else if (rs1 > 0 .and. rs2 == 0) then

gamma = theta(-pa +ptr%p%n) + atan(coss/sine)
gamma = theta(-pa +ptr%p%n) + atan2(coss,sine)
aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ( (sigma*(1+B*sin(beta*gamma))/r)**6 &
* (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * &
[(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
Expand All @@ -1437,7 +1437,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,
ptrn%p%F = -aux2 + ptrn%p%F - aux3*[sine,-coss]
! partlst(ptr%p%n - pa) = 1
else if (rs2 > 0 .and. rs1 == 0) then
gamma = theta(-pa +ptrn%p%n) + atan(coss/sine)
gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine)
aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ((sigma*(1+B*sin(beta*gamma))/r)**6 * &
(1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * &
[(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
Expand Down Expand Up @@ -1516,7 +1516,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,
ptr%p%F = aux2 + ptr%p%F +fR
ptrn%p%F = -aux2 + ptrn%p%F - fR
else if (rs1 > 0 .and. rs2 == 0) then
gamma = theta(-pa +ptr%p%n) + atan(coss/sine)
gamma = theta(-pa +ptr%p%n) + atan2(coss,sine)

aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * &
( (sigma*(1+B*sin(beta*gamma))/r)**6 * (1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * &
Expand All @@ -1537,7 +1537,7 @@ subroutine comp_FT(GField,Hfield,theta,Tor,pr,mesh,malha,propriedade,r_cut,domx,
! partlst(ptr%p%n - pa) = 1

else if (rs2 > 0 .and. rs1 == 0) then
gamma = theta(-pa +ptrn%p%n) + atan(coss/sine)
gamma = theta(-pa +ptrn%p%n) + atan2(coss,sine)
aux2 = -24*epsil*(1+A*sin(alpha*gamma+ph))*(1/r**2) * ((sigma*(1+B*sin(beta*gamma))/r)**6 * &
(1-2* (sigma*(1+B*sin(beta*gamma))/r)**6 )) * &
[(x1(1)-x2(1)) - (rs1+rs2)*coss, (x1(2)-x2(2)) - (rs1+rs2)*sine]
Expand Down Expand Up @@ -6516,6 +6516,7 @@ program main
end if
!! CASE OF RUGGED LENNARD JONES THAT CAN MAKE THE PARTICLE ROTATE
if (abs(pr(4)+pr(3)) > 0) then ! PARTICLE ROTATES
print*, "CASE: Particle Rotates"
if (id == 0) then
deallocate(x,v)
allocate(x(N,3),v(N,3))
Expand Down Expand Up @@ -6826,6 +6827,7 @@ program main
end do

else if (termostato_nose_hoover) then
print*, "CASE: NOSE HOOVER"
! TERMOSTATO NOSE HOOVER
do while (t_fim > t)
!print*, "L 1990"
Expand Down Expand Up @@ -7010,6 +7012,7 @@ program main

else
! TERMOSTATO SCALING OU SEM TEMOSTATO
print*, "CASE: SCALING OR NO TERMOSTAT"
do while (t_fim > t)
! COMPF
call comp_F(GField, mesh,malha,propriedade,rcut,domx,domy,ids,id,t) !altera Força
Expand Down
38 changes: 25 additions & 13 deletions rot_par.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
import numpy as np
from numpy import sin, cos
from numpy import sin, cos, heaviside
import matplotlib.pyplot as plt


Expand All @@ -9,24 +9,29 @@ def abssin(x):
maxF = 40
sigma_ref = 0.341

delta = 0.005

delta = 0.01
sig0 = 1
ep0 = 4
ep0 = 1
rs = 5 # raio solido
rm = .245 / sigma_ref #raio do atomo
# .245 é a distância entra atomos em estrutura hexagonal rombica

x = y = np.arange(0, 3*sig0+rs, delta)
# caso com polimero
Np = 6
deltad = 2

x = y = np.arange(-(3*sig0+rs), 3*sig0+rs, delta)
X, Y = np.meshgrid(x, y)

r = np.sqrt(X**2 + Y**2)
np.place(r, r < rs, np.nan)
gamma = np.arctan2(Y,X)


A = 0.0001 #ep0/2 #(0.071/sigma_ref)/4 # no epsilon
B = 0.0001 #(0.071/sigma_ref)/8 # no sigma
alpha = beta = (2*np.pi*rs/rm)
A = 0 #ep0/8 #(0.071/sigma_ref)/4 # no epsilon
B = 0 #(0.071/sigma_ref)/8 # no sigma
alpha = beta = 6 # (2*np.pi*rs/rm)/2
ph = (2*np.pi/alpha)
print("A {}, B {}, alpha {}, fase {}, periodo {}\n".format(A, B, alpha ,ph,2*np.pi/alpha))

Expand All @@ -35,17 +40,20 @@ def abssin(x):
gmolec = np.linspace(0,2*np.pi/4,round((2*np.pi*rs/rm)/4))
#beta = 6

V = 4*ep0*(1+A*np.sin(alpha*gamma+ph))* ( (sig0*(1+B*np.sin(beta*gamma))/(r-rs) )**12 - (sig0*(1+B*np.sin(beta*gamma))/(r-rs))**6 )
#V = 4*ep0*(1+A*np.sin(alpha*gamma+ph))* ( (sig0*(1+B*abs(np.sin(beta*gamma)))/(r-rs) )**12 - (sig0*(1+B*abs(np.sin(beta*gamma)))/(r-rs))**6 )
# LJ rugoso
#V = 4*ep0*(1+A*np.sin(alpha*gamma+ph))* ( (sig0*(1+B*np.sin(beta*gamma))/(r-rs) )**12 - (sig0*(1+B*np.sin(beta*gamma))/(r-rs))**6 )
# com polimero
V = 4*ep0*(1+A*np.sin(alpha*gamma+ph))* ( (sig0*(1+B*np.sin(beta*gamma))/(r-rs) )**12 - (sig0*(1+B*np.sin(beta*gamma))/(r-rs))**6 ) + (Np/(2*deltad))*((deltad*2 + sig0) + (r-rs) *(np.log((deltad*2 + sig0)/(r-rs) ) -1))*np.sin(beta*gamma)**2*heaviside(r-(rs+sig0) , 1)*(1-heaviside(r-(deltad*2+rs+sig0) , 1))

#V = 4*ep0*(1+A*np.sin(alpha*gamma+ph))* ( (sig0*(1+B*abs(np.sin(beta*gamma)))/(r-rs) )**12 - (sig0*(1+B*abs(np.sin(beta*gamma)))/(r-rs))**6 )
#V = 4*ep0*(1+A*np.sin(alpha*gamma+ph))* ( (sig0/(r-rs*(1+B*np.sin(beta*gamma))) )**12 - (sig0/(r-rs*(1+B*np.sin(beta*gamma))))**6 )

#Fr = 24*ep0*(1+A*np.sin(alpha*gamma+ph))*(1/(r-rs)**2) * ( (sig0*(1+B*np.sin(beta*gamma))/(r-rs))**6 * (1-2* (sig0*(1+B*np.sin(beta*gamma))/(r-rs))**6 )) * (r-rs)
# Fg = (1/(r-rs))* ( (4*ep0*(A*np.sin(alpha*gamma+ph) + 1 ))* (sig0/(r-rs))**6*(B*(np.sin(beta*gamma))+1)**5 * \
# ( (sig0/(r-rs))**6 * 12*beta*B*np.cos(beta*gamma)*(B*(np.sin(beta*gamma))+1)**6 - 6*beta*B*np.cos(beta*gamma)) + \
# 4*alpha*A*ep0*np.cos(alpha*gamma+ph)*(sig0/(r-rs))**6*(B*(np.sin(beta*gamma))+1)**6 * \
# ((sig0/(r-rs))**6 * (B*(np.sin(beta*gamma))+1)**6 - 1) )

#rugoso
Fr = 24*ep0*(1+A*sin(alpha*gamma+ph))*(1/(r-rs)**2) * ((sig0*(1+B*sin(beta*gamma))/(r-rs))**6 * \
(1-2* (sig0*(1+B*sin(beta*gamma))/(r-rs))**6 )) * (r-rs)

Expand All @@ -54,6 +62,10 @@ def abssin(x):
4*alpha*A*ep0*np.cos(alpha*gamma+ph)*(sig0/(r-rs))**6*(B*np.sin(beta*gamma)+1)**6 * \
((sig0/(r-rs))**6 * (B*np.sin(beta*gamma)+1)**6 - 1) )

# polimero
Fr = 24*ep0*(1/(r-rs)**2) * ((sig0/(r-rs))**6 *(1-2* (sig0/(r-rs))**6 )) * (r-rs) + (Np/(2*deltad))*((np.log((deltad*2 + sig0)/(r-rs) ) -1) + (deltad*2+sig0))*np.sin(beta*gamma)**2*heaviside(r-(rs+sig0) , 1)*(1-heaviside(r-(deltad*2+rs+sig0) , 1))
Fg = 2*beta*cos(beta*gamma)*sin(beta*gamma)*(Np/(2*deltad))*((deltad*2 + sig0) + (r-rs) *(np.log((deltad*2 + sig0)/(r-rs) ) -1))*np.sin(beta*gamma)**2*heaviside(r-(rs+sig0) , 1)*(1-heaviside(r-(deltad*2+rs+sig0) , 1))
###
# Fr = 24*ep0*(1+A*np.sin(alpha*gamma+ph))*(1/(r-rs)**2) * ( (sig0*(1+B*abs(np.sin(beta*gamma)))/(r-rs))**6 * (1-2* (sig0*(1+B*abs(np.sin(beta*gamma)))/(r-rs))**6 )) * (r-rs)

# Fg = (1/(r-rs))* ( (4*ep0*(A*np.sin(alpha*gamma+ph) + 1 ))* (sig0/(r-rs))**6*(B*abs(np.sin(beta*gamma))+1)**5 * \
Expand All @@ -66,9 +78,9 @@ def abssin(x):
np.place(Fr,abs(Fr) > maxF, np.nan)
np.place(Fg,abs(Fg) > maxF, np.nan)

plt.figure()
plt.plot(X[int(len(X)/2),:] ,V[int(len(X)/2),:])
plt.title("potential")
# plt.figure()
# plt.plot(X[int(len(X)/2),:] ,V[int(len(X)/2),:])
# plt.title("potential")


fig1, ax1 = plt.subplots()
Expand Down

0 comments on commit eb87dca

Please sign in to comment.