Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.mso.anu.edu.au/pfrancis/simulations/co2.py Дата изменения: Fri Oct 25 13:34:33 2013 Дата индексирования: Sun Apr 10 03:25:45 2016 Кодировка: Поисковые слова: annular

from __future__ import division
from visual import *
import math

spacing = 1.0E-10
k = 8.0E-4 # Tension compression constant
kt = 3.0E-5 # Torsion constant
u = 1.660E-27
ioncharge = 1.0*1.602E-19
rad = 3.0E-11

scale = 2.0*spacing
scene = display(range=(scale,scale,scale),width=1200,height=800, background=color.black)

theta = 30.*math.pi/180.0

oxy1 = sphere()
oxy1.pos = vector(spacing*math.cos(theta),spacing*math.sin(theta),0.)
oxy1.mass = 16.0*u
oxy1.charge = ioncharge
oxy1.p = vector(0.,0.,0.)
oxy1.color=color.green
oxy1.radius = rad
oxy1pos0 = oxy1.pos[1]

oxy2 = sphere()
oxy2.pos = -1.0*oxy1.pos
oxy2.mass = 16.0*u
oxy2.charge = ioncharge
oxy2.p = vector(0.,0.,0.)
oxy2.color=color.green
oxy2.radius = rad

carb1 = sphere()
carb1.pos = vector(0.,0.0*spacing,0.)
carb1.mass = 12.0*u
carb1.charge = -2.0*ioncharge
carb1.p = vector(0.,0.,0.)
carb1.color=color.red
carb1.radius = 0.8*rad


spring1 = helix(pos=carb1.pos,axis=oxy1.pos-carb1.pos, radius = 0.3*rad,
thickness = 0.1*rad, color=color.yellow, coils = 10)

spring2 = helix(pos=carb1.pos,axis=oxy2.pos-carb1.pos, radius = 0.3*rad,
thickness = 0.1*rad, color=color.yellow, coils = 10)

evector = arrow(pos=(1.7*spacing,0.,0.), axis=(0.,0.,0.),
shaftwidth=0.1*spacing, color=color.blue)


E0 = 200.0 #200

frequency = 8.0 # 5.291 for linear, 1.5 for bending, 2.76 for non active sym
frequency *= (2.0*math.pi)*1.00001E10


timestep = 1.0E-13

time = 0
flag = 1
old = 0.
while 1:
rate(500)

time += timestep
E = E0*math.cos(time*frequency)*vector(0.,1.,0.)
evector.axis=0.5*spacing*E/E0

if (flag == 1) & (carb1.pos[1] < 0):
new = time
flag = 0
#print 1.0e-10/(new-old)
elif (flag == 0) & (carb1.pos[1] > 0):
flag = 1
old = new


v1 = oxy1.pos - carb1.pos
v2 = oxy2.pos - carb1.pos

angle = math.acos(dot(v1, -1.0*v2)/(v1.mag*v2.mag))

# Oxy1 forces
coulombforce = oxy1.charge*E
stretchforce = -1.0*k*(mag(v1)-spacing)*norm(v1)
torqueforce = -1.0*kt*spacing*angle*norm(v1+v2)

oxy1.p += timestep*(coulombforce+stretchforce+torqueforce)
oxy1.pos += timestep*oxy1.p/oxy1.mass

# Oxy2 forces
coulombforce = oxy2.charge*E
stretchforce = -1.0*k*(mag(v2)-spacing)*norm(v2)
torqueforce = -1.0*kt*spacing*angle*norm(v1+v2)

oxy2.p += timestep*(coulombforce+stretchforce+torqueforce)
oxy2.pos += timestep*oxy2.p/oxy2.mass

# Carb1 forces
coulombforce = carb1.charge*E
stretchforce = k*(mag(v1)-spacing)*norm(v1) + k*(mag(v2)-spacing)*norm(v2)
torqueforce = 2.0*kt*spacing*angle*norm(v1+v2)

carb1.p += timestep*(coulombforce+stretchforce+torqueforce)
carb1.pos += timestep*carb1.p/carb1.mass

# Update springs
spring1.pos = carb1.pos
spring1.axis = v1

spring2.pos = carb1.pos
spring2.axis = v2