修订版 | 07ecc6cc5344d8d718d0fbfda3e83526854440d0 (tree) |
---|---|
时间 | 2008-02-26 02:42:18 |
作者 | iselllo |
Commiter | iselllo |
I added videcoq.py which calculates the DLVO potential used by Videcoq.
@@ -0,0 +1,67 @@ | ||
1 | +#! /usr/bin/env python | |
2 | + | |
3 | +import scipy as s | |
4 | +import numpy as n | |
5 | +import pylab as p | |
6 | + | |
7 | + | |
8 | + | |
9 | + | |
10 | +def videcoq_pot(A,B,d,k,r): | |
11 | + V_vdw=s.zeros(len(r)) | |
12 | + V_el=s.zeros(len(r)) | |
13 | + | |
14 | + #V[:]=s.where(r<=r_core,slope*r+intercept,V[:]) | |
15 | + V_vdw=-A*(d**2./(r**2.-d**2.)+(d/r)**2.+2*s.log((r**2.-d**2.)/r**2.)) | |
16 | + V_el=B*s.exp(-k*(r-d)) | |
17 | + | |
18 | + V_tot=V_vdw+V_el | |
19 | + | |
20 | + return V_tot | |
21 | + | |
22 | + | |
23 | +#Now it is time to calculate some parameters: | |
24 | + | |
25 | + | |
26 | +d=500e-9 # particle diameter | |
27 | +A=4.76e-20/12. | |
28 | +print "A is, ", A | |
29 | + | |
30 | + | |
31 | +# | |
32 | +T=300. #temperature in kelvin | |
33 | + | |
34 | +k_B=1.38e-23 | |
35 | +z=1. | |
36 | +e=1.6e-19 #electron charge in Coulombs | |
37 | +epsi_0=8.85e-12 #vacuum permittivity | |
38 | +epsi_r=81. | |
39 | +psi=32e-3 #in Volts | |
40 | +k=3.35e8 | |
41 | + | |
42 | + | |
43 | +B=s.pi*epsi_0*epsi_r*(4.*k_B*T/(z*e)*s.tanh(z*e*psi/(4.*k_B*T)))**2. | |
44 | +B=B*d | |
45 | + | |
46 | +r=s.linspace(1.01,6.,1000) | |
47 | +r=r*d | |
48 | + | |
49 | + | |
50 | + | |
51 | +potential=videcoq_pot(A,B,d,k,r) | |
52 | + | |
53 | +p.plot(r/d,potential/(k_B*T),linewidth=2.) | |
54 | +p.ylim(min(potential/(k_B*T))-0.2,max(0.1,max(potential/(k_B*T)))) | |
55 | +p.xlim(0.9,max(r/d)) | |
56 | +p.xlabel('r') | |
57 | +p.ylabel('My Radial Potential') | |
58 | +#pylab.legend(('prey population','predator population')) | |
59 | +p.title('Radial Potential') | |
60 | +p.grid(True) | |
61 | +p.savefig('videcoq_dimensional_potential.pdf') | |
62 | +p.hold(False) | |
63 | +p.clf() | |
64 | + | |
65 | + | |
66 | +print "So far so good" | |
67 | + |