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

from __future__ import division
from visual import *
import math

spacing = 1.0E-10
k = 5.0E-6 # Tension compression constant
kt = 5.0E-6 # 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 = 0.0*math.pi/180.0
bondangle = 105.*math.pi/180.0
tilt = 20.0*math.pi/180.0

hy1 = sphere()
hy1.pos = rotate(vector(spacing*math.cos(theta),spacing*math.sin(theta),0.),
angle = tilt, axis=(1,0,0))
hy1.mass = 1.0*u
hy1.charge = ioncharge
hy1.p = vector(0.,0.,0.)
hy1.color=color.green
hy1.radius = 0.4*rad

hy2 = sphere()
hy2.pos = rotate(vector(spacing*math.cos(theta+bondangle),
spacing*math.sin(theta+bondangle),0.),
angle=tilt, axis=(1,0,0))
hy2.mass = 1.0*u
hy2.charge = ioncharge
hy2.p = vector(0.,0.,0.)
hy2.color=color.green
hy2.radius = 0.4*rad

oxy1 = sphere()
oxy1.pos = vector(0.,0.0*spacing,0.)
oxy1.mass = 16.0*u
oxy1.charge = -2.0*ioncharge
oxy1.p = vector(0.,0.,0.)
oxy1.color=color.red
oxy1.radius = 1.2*rad

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

spring2 = helix(pos=oxy1.pos,axis=hy2.pos-oxy1.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 = 120.0 #100

frequency = 0.25 # units of 10^10 Hz, 2.8 is good

frequency *= (2.0*math.pi)*1.000001E10


timestep = 5.0E-13

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

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

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


v1 = hy1.pos - oxy1.pos
v2 = hy2.pos - oxy1.pos

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

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

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

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

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

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

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

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

spring2.pos = oxy1.pos
spring2.axis = v2