Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.mso.anu.edu.au/pfrancis/phys1101/Lectures/L09/threebody.py Дата изменения: Thu Mar 10 08:11:14 2011 Дата индексирования: Tue Oct 2 14:20:24 2012 Кодировка:

from __future__ import division
from visual import *
R=4e11
s1 = sphere(pos=(0,R,0), radius=1e10, color=color.magenta)
s2 = sphere(pos=(0,-R,0), radius=1e10, color=color.cyan)
s3 = sphere(pos=(R,0.05*R,0), radius=1e10, color=color.green)# Change to 0.04
s1.m = 1e30
s2.m = 1e30
s3.m = 1e30
G = 6.7e-11
#v = sqrt(G*s1.m/(4*R))
v = 10000
s1.p = s1.m*vector(-0.2*v,-0.3*v,0)
s2.p = s2.m*vector(0,0.2*v,0)
s3.p = s3.m*vector(-0.15*v,0,0)
t1=curve(color=s1.color)
t2=curve(color=s2.color)
t3=curve(color=s3.color)
dt = 0.6*1000
t = 0
while t < 100000*2000:
#rate(1000)
r21 = s2.pos - s1.pos
F21 = -norm(r21)*G*s1.m*s2.m/mag(r21)**2
r32 = s3.pos - s2.pos
F32 = -norm(r32)*G*s3.m*s2.m/mag(r32)**2
r31= s3.pos - s1.pos
F31 = -norm(r31)*G*s3.m*s1.m/mag(r31)**2
s3.p = s3.p + (F32+F31)*dt
s2.p = s2.p + (F21-F32)*dt
s1.p = s1.p + (-F21-F31)*dt
s1.pos = s1.pos + (s1.p/s1.m)*dt
s2.pos = s2.pos + (s2.p/s2.m)*dt
s3.pos = s3.pos + (s3.p/s3.m)*dt
t1.append(pos=s1.pos)
t2.append(pos=s2.pos)
t3.append(pos=s3.pos)
t = t+dt
print s1.pos, s2.pos, s3.pos