Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.mso.anu.edu.au/pfrancis/simulations/raytrace_caustics.py Дата изменения: Tue May 5 10:17:11 2015 Дата индексирования: Sun Apr 10 05:30:46 2016 Кодировка: Поисковые слова: http astrokuban.info astrokuban

import numpy
import pylab

# Set up incident rays
irange = 1.0
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 = 1.0
xmass2 = -0.4
xvel = 0.0
ymass2 = 0.0
yvel = 0.0

nplots = 1000
npoints = nrays/nplots
zeros = numpy.zeros(npoints)
x1 = numpy.array([])
y1 = numpy.array([])
for i in range(0,nplots):
print "Generating plot ", i
imname = "plots2/im" + '%(#)04d' % {"#": i} + '.png'

x = -1.0*irange + 2*i*irange/nplots
xtemp = zeros+x
ytemp = numpy.arange(-irange,irange,2*irange/npoints)

# Deflect by first mass
delx = xtemp - xmass1
dely = ytemp - 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

delx = xtemp - xmass2
dely = ytemp - ymass2
r = numpy.sqrt(delx*delx+dely*dely)
dtheta = mass2*0.01/r
dx += delx*dtheta/r
dy += dely*dtheta/r

x2 = xtemp - dx
y2 = ytemp - dy
x1 = numpy.append(x1,x2)
y1 = numpy.append(y1,y2)

#pylab.plot(x1,y1,"+k")
pylab.clf()
pylab.plot(x2,y2,".k",markersize=0.4) #0.03
pylab.plot(x1,y1,".g",markersize=0.03) #0.03

pylab.axis("equal")
pylab.xlim([-0.4,0.4])
pylab.ylim([-0.4,0.4])

pylab.savefig(imname)