Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.mso.anu.edu.au/pfrancis/simulations/raytrace.py Дата изменения: Thu Oct 17 03:26:24 2013 Дата индексирования: Sun Apr 10 04:24:24 2016 Кодировка: Поисковые слова: http www.astronomy.ru forum index.php topic 4644.0.html

import numpy
import pylab

# Set up incident rays
irange = 2.5
nrays = 1000000
x1 = numpy.random.uniform(-irange,irange,nrays)
y1 = numpy.random.uniform(-irange,irange,nrays)

#Set up first lens
mass1 = 20.0
xmass1 = 0.0
ymass1 = 0.0

#Set up second lens
mass2 = 0.0
xmass2 = -1.2
xvel = 0.02
ymass2 = -0.15
yvel = 0.0

nplots = 1
for i in range(0,nplots):
print "Generating plot ", i
imname = "movie/im" + '%(#)04d' % {"#": i} + '.png'

# Deflect by first mass
delx = x1 - xmass1
dely = y1 - ymass1
r = numpy.sqrt(delx*delx+dely*dely)
dtheta = mass1*0.01/r
dx = delx*dtheta/r
dy = dely*dtheta/r

# Deflect by second mass
xmass2 += xvel
ymass2 += yvel
delx = x1 - xmass2
dely = y1 - ymass2
r = numpy.sqrt(delx*delx+dely*dely)
dtheta = mass2*0.01/r
dx += delx*dtheta/r
dy += dely*dtheta/r

x2 = x1 - dx
y2 = y1 - dy

#pylab.plot(x1,y1,"+k")
pylab.clf()
pylab.plot(x2,y2,".b",markersize=0.1)
pylab.axis("equal")
pylab.xlim([-1,1])
pylab.ylim([-1,1])

pylab.savefig(imname)