PY525 hw 3

import numpy

boxsize = 10.
boxdiv = (10-10*0.1)/2
num = 16
snum = numpy.int(numpy.sqrt(num))
x = numpy.linspace(-boxdiv,boxdiv,snum)
y = numpy.linspace(-boxdiv,boxdiv,snum)
dt = 0.01
time = 4
ttime = numpy.arange(0,time+0.00001,dt)

p = numpy.zeros((num,ttime.size,2))

for i in range(0,snum,1):
    for j in range(0,snum,1):
        p[(i+j*snum),0,0,] = x[i]
        p[(i+j*snum),0,1] = y[j]

#Choose random velocities
v = numpy.zeros((num,2))

for i in range(0,num,1):
    flipx = numpy.random.rand()
    if flipx > .5: v[i,0] = flipx*3.5
    else: v[i,0] = flipx*-3.5
    flipy = numpy.random.rand()
    if flipy > .5: v[i,1] = flipy*3.5
    else: v[i,1] = flipx*-3.5

for i in range(1,ttime.size,1):
    for j in range(0,num,1):
        #Assign positions based on velocity
        position = p[j,i-1,:]

        for k in range(0,num,1):
            r = numpy.sqrt((position[0]-p[k,i-1,0])**2+(position[1]-p[k,i-1,1])**2)
            if (r <= 3) and (j != k):
                v[j,0] = -24*(2*(1/(r**13))+1/(r**7))*dt
                v[j,1] = -24*(2*(1/(r**13))+1/(r**7))*dt
        position[0] += v[j,0]*dt
        position[1] += v[j,1]*dt

        #Flip position if on a border
        if (position[0] >= 5 or position[0] <= -5):
            position[0] = position[0]*-1
        if (position[1] >= 5 or position[1] <= -5):
            position[1] = position[1]*-1

        p[j,i,:] = position