0

I wrote the following code,

import numpy as np
from random import gauss
from random import seed
from pandas import Series
import matplotlib.pyplot as plt 
import math
###variable declaration
R=0.000001 #radiaus of particle
beta=0.23 # shape factor
Ad=9.2#characteristic nanoscale defect area
gamma=4*10**-2 #surface tension
tau=.0001 #line tension

phi=-(math.pi/2)#spatial perturbation
lamda=((Ad)/(2*3.14*R)) #averge contcat line position

mu=0.001#viscosity of liquid
lamda_m=10**-9# characteristic size of adsorption site
KbT=(1.38**-24)*293 # boltzman constant with tempartaure
nu=0.001#moleculer volume in liquid phase 1
khi=3 #scaling factor
#deltaF=(beta*gamma*Ad)#surface energy perturbation
deltaF=19*KbT
# seed random number generator
seed(0)
# create white noise series
series = [gauss(0.0, 1.0) for i in range(1)]
series = Series(series)
#########################################
Z=0.0000001 #particle position
time=1
dt=1
for time in np.arange(1, 100, dt):
#####simulation loop#######
    theta=np.arccos(-Z/R) #contact angle
    theta_e=((math.pi*110)/180) #equilibrium contact angle
    Z_e=-R*np.cos(theta_e)#equilibrium position of particle
    C=3.14*gamma*(R-Z_e) #additive constant
    Fsz= (gamma*math.pi*(Z-Z_e)**2)+(tau*2*math.pi*math.sqrt(R**2-Z**2))+C
    Fz=Fsz+(0.5*deltaF*np.sin((2*math.pi/lamda)*(Z-Z_e)-phi))#surface force
    #dFz=(((gamma*Ad)/2)*np.sin(2*math.pi/lamda))+((Z-Z_e)*(2*gamma*math.pi))-((tau*2*math.pi*Z)/(math.sqrt(R**2-Z**2)))
    dFz=(deltaF*np.sin(2*math.pi/lamda))+((Z-Z_e)*(2*gamma*math.pi))-((tau*2*math.pi*Z)/(math.sqrt(R**2-Z**2)))
    w_a=gamma*lamda_m**2*(1-np.cos(theta_e)) #work of adhesion
    epsilon_z=2*math.pi*R*np.sin(theta)*mu*(nu/(lamda_m**3))*np.exp(w_a/KbT)#transitional drag
    epsilon_s=khi*mu*((4*math.pi**2*R**2)/math.sqrt(Ad))*(1-(Z/R)**2)
    epsilon=epsilon_z+epsilon_s
    Ft=math.sqrt(2*KbT*epsilon)*series #thermal force
    v=(dFz+Ft)/epsilon ##new velocity
    Z=Z+v*dt #new position    
    print('z=',Z)
    print('v=',v)
    print('Fz=',Fz)
    print('dFz',dFz)
    print('time',time)    
plt.plot(Z,time)
plt.show()

According to my code I suppose to have 99 values for everything (Fz, Z, v , time). While I print , I can see all the values but while I was trying to plot them with different parameters with respect to each other for analyzing, I never get any graph. Can anyone tell me, what is missing from my code with explanation?

2
  • I think that either your indentation is wrong -- right now you only plot the last Z and time value of your for loop. Anyway, without really trying to dig into your math, I think you can do the entire thing without for loop using numpy. What is the variable series for? I don't see it being used anywhere. Commented Jan 19, 2018 at 8:35
  • @ThomasKühn the variable series is used to generate random noise. It is used to calculate Ft Commented Jan 19, 2018 at 19:30

2 Answers 2

1

@AnttiA's answer is basically correct, but can be easily misunderstood, as can be seen from the OP's comment. Therefore here the complete code altered such that a plot is actually produced. Instead of making Z a list, define another variable as list, say Z_all = [], and then append the updated Z-values to that list. The same can be done for the time variable, i.e. time_all = np.arange(1,100,dt). Finally, take the plot command out of the loop and plot the entire data series at once.

Note that in your example you don't really have a series of random numbers, you pull one fixed number for one fixed seed, thus the plot is not really meaningful (it appears to be producing a straight line). Trying to interpret your intentions correctly, you probably want a series of random numbers that is as long as your time series. This is most easily done using np.random.normal

There are a lot of other ways that your code could be optimised. For instance, all mathematical functions from the math module are also found in numpy, so you could just not import math at all. The same goes for pandas. Also, you define some constant values inside the for-loop, which could be computed once before the loop. Lastly, @AnttiA is probably right, that you want time on the x axis and Z on the y axis. I therefore generate two plots -- on the left time over Z and on the right Z over time. Now finally the altered code:

import numpy as np
#from random import gauss
#from random import seed
#from pandas import Series
import matplotlib.pyplot as plt 
#import math

###variable declaration
R=0.000001 #radiaus of particle
beta=0.23 # shape factor
Ad=9.2#characteristic nanoscale defect area
gamma=4*10**-2 #surface tension
tau=.0001 #line tension

phi=-(np.pi/2)#spatial perturbation
lamda=((Ad)/(2*3.14*R)) #averge contcat line position

mu=0.001#viscosity of liquid
lamda_m=10**-9# characteristic size of adsorption site
KbT=(1.38**-24)*293 # boltzman constant with tempartaure
nu=0.001#moleculer volume in liquid phase 1
khi=3 #scaling factor
#deltaF=(beta*gamma*Ad)#surface energy perturbation
deltaF=19*KbT


##quantities moved out of the for-loop:
theta_e=((np.pi*110)/180) #equilibrium contact angle
Z_e=-R*np.cos(theta_e)#equilibrium position of particle
C=3.14*gamma*(R-Z_e) #additive constant
w_a=gamma*lamda_m**2*(1-np.cos(theta_e)) #work of adhesion

#########################################
Z=0.0000001 #particle position
##time=1
dt=1

Z_all = []
time_all = np.arange(1, 100, dt)
# seed random number generator
# seed(0)
np.random.seed(0)
# create white noise series
##series = [gauss(0.0, 1.0) for i in range(1)]
##series = Series(series)
series = np.random.normal(0.0, 1.0, len(time_all))

for time, S in zip(time_all,series):
#####simulation loop#######
    Z_all.append(Z)
    theta=np.arccos(-Z/R) #contact angle
    Fsz= (gamma*np.pi*(Z-Z_e)**2)+(tau*2*np.pi*np.sqrt(R**2-Z**2))+C
    Fz=Fsz+(0.5*deltaF*np.sin((2*np.pi/lamda)*(Z-Z_e)-phi))#surface force
    #dFz=(((gamma*Ad)/2)*np.sin(2*np.pi/lamda))+((Z-Z_e)*(2*gamma*np.pi))-((tau*2*np.pi*Z)/(np.sqrt(R**2-Z**2)))
    dFz=(deltaF*np.sin(2*np.pi/lamda))+((Z-Z_e)*(2*gamma*np.pi))-((tau*2*np.pi*Z)/(np.sqrt(R**2-Z**2)))
    epsilon_z=2*np.pi*R*np.sin(theta)*mu*(nu/(lamda_m**3))*np.exp(w_a/KbT)#transitional drag
    epsilon_s=khi*mu*((4*np.pi**2*R**2)/np.sqrt(Ad))*(1-(Z/R)**2)
    epsilon=epsilon_z+epsilon_s
    Ft=np.sqrt(2*KbT*epsilon)*S #series #thermal force
    v=(dFz+Ft)/epsilon ##new velocity
    Z=Z+v*dt #new position    
    print('z=',Z)
    print('v=',v)
    print('Fz=',Fz)
    print('dFz',dFz)
    print('time',time)    


fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(8,4))

axes[0].plot(Z_all,time_all)
axes[0].set_xlabel('Z')
axes[0].set_ylabel('t')

axes[1].plot(time_all, Z_all)
axes[1].set_xlabel('t')
axes[1].set_ylabel('Z')

fig.tight_layout()

plt.show()

The result looks like this:

result of above code

Sign up to request clarification or add additional context in comments.

1 Comment

@Thomas...this is the exact thing i wanted..actually I find a solution as well but problem remain with straight line that you have explained very clearly. Also regarding programming your explanation is very clear to me. good job sir.
1

I suppose, you will get plot anyway, y-values are probably from 94 to 104. Now you are plotting line with one point. Its length is zero, that's why you cannot see it, try: plt.plot(Z,time,'*'). Now you should get graph with an asterix in the middle.
As Thomas suggested, you should use arrays instead of using last calculated value. If you prefer loops (sometimes they are easier to modify), modify the lines...
Before loop:

Z = [0.0000001] # Initialize Z for time 0
time_vec = np.arange(1, 100, dt)

Inside loop:

Z.append(Z[-1] + v*dt) # new position

After loop:

plt.plot(Z[1:], time_vec)

Have no time to test it, hopefully works...

Note that first argument in plot command is x-axis values and second y-axis, I'd prefer time in x-axis.

1 Comment

This solution give me '' theta=np.arccos(-Z/R) #contact angle TypeError: bad operand type for unary -: 'list' '' this error.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.