views:

85

answers:

1

Hi,

I have a programming question, as follows, for which my solution does not produce the desired output

This particle simulator operates in a universe with different laws of physics to ours. Each particle has a position (x, y), velocity (vx, vy) and an acceleration (ax, ay). Every particle exerts an attractive force on every other particle. This force is the same regardless of how far away the particle is.
The acceleration of a particle in the x direction is given by

(ax=number of particles to right of x−number of particles to left of x) / 10.0

The particle will then move left or right with velocity vx + ax.
Similarly, the acceleration of a particle in the y direction is given by

(ay=number of particles above y−number of particles below y) / 10.0

The particle will then move up or down with velocity vy + ay.
The particles are bound in a chamber with dimensions −300 < x < 300 and −200 < y < 200. If a particle hits the wall of the chamber it should bounce. Bouncing involves setting the x or y coordinate to the boundary, and reversing the direction of the velocity. For example, if a particle ends up with position x=305 then you should set x=300 and vx = -vx. Note that x must be set to the integer value 300 in order to get the same output values as our test cases.

Write a program to read in a file named particles.txt containing the initial positions, velocities and accelerations of a number of particles. The first line of the file contains the number of iterations to run the simulation for (5 in this example). Every other row contains data about one particle, in the format x y vx vy ax ay like this:

5
0 -30 3 0 0 0
100 50 0 1 0 0
20 10 0 3 0 0
-80 15 2 -2 0 0

Your program should create a Particle object to store the data for each particle. Then for every iteration of the simulation you should

  • calculate the acceleration on each particle (using the equations above)
  • then calculate the new velocity for each particle (vx = vx + ax)
  • then calculate the new position for each particle (x = x + vx)
  • The output of your program should be a list of positions for each particle at every step of the simulation, in CSV format x1,y1,x2,y2,x3,y3,x4,y4 for the 4-particle example shown above:
    3.1,-29.7,99.7,50.7,19.9,13.1,-77.7,12.9
    6.3,-29.1,99.1,51.1,19.7,16.1,-75.1,10.9
    9.6,-28.2,98.2,51.2,19.4,19.0,-72.2,9.0
    13.0,-27.0,97.0,51.0,19.0,21.8,-69.0,7.2
    16.5,-25.5,95.5,50.5,18.5,24.5,-65.5,5.5

    To produce these numbers, you should call str on the x and y coordinates of each of the particles.

    My code is as follows:

    class Particle(object):
        def __init__(self, (x, y, vx, vy, ax, ay)):
            # Set initial values
            self.x, self.y = float(x), float(y)
            self.vx, self.vy = float(vx), float(vy)
            self.ax, self.ay = float(ax), float(ay)
    
        def __str__(self):
            return '(' + str(self.x) + ', ' + str(self.y) + ')'
    
        # Calculate new acceleration
        def calc_acc(self, part_list):
            left, right = 0, 0
            up, down = 0, 0
            for particle in part_list:
                # Count particles on left & right
                if particle.x < self.x:
                    left += 1
                elif particle.x > self.x:
                    right += 1
                # Count particles above & below
                if particle.y > self.y:
                    up += 1
                elif particle.y < self.y:
                    down += 1
            # Calculate x-acceleration
            self.ax = (right - left) / 10.0
            # Calculate y-acceleration
            self.ay = (up - down) / 10.0
    
        # Calculate new velocity
        def calc_vel(self):
            self.vx = self.vx + self.ax
            self.vy = self.vy + self.ay
    
        # Move the particle
        def move(self, p_list):
            # Calculate new acceleration & velocity
            self.calc_acc(p_list)
            self.calc_vel()
            # Make move
            self.x += self.vx
            self.y += self.vy
            # Check for bounce, and bounce
            # X-axis
            if self.x > 300.0:
                self.x = 300
                self.vx = -(self.vx)
            elif self.x < -300.0:
                self.x = -300
                self.vx = -(self.vx)
            # Y-axis
            if self.y > 200.0:
                self.y = 200
                self.vy = -(self.vy)
            elif self.y < -200.0:
                self.y = -200
                self.vy = -(self.vy)
            # Return resulting position
            return self.x, self.y
    
    # Take file input
    input_file = []
    for line in open('particles2.txt', 'rU'):
        input_file.append(line)
    
    # Take number of iterations, and leave particle info only
    times = int(input_file.pop(0))
    
    # Remove '\n' from particle info, and convert particle info to a list of floats
    for line in input_file:
        input_file[input_file.index(line)] = line.strip('\n').split()
    
    # Create list of Particle objects
    particles = []
    for line in input_file:
        particles.append(Particle(line))
    
    # Create result position array
    results = []
    for i in range(times):
        results.append([None for x in range(len(particles))])
    
    # Run simulation for 'times' iterations
    for iteration in range(times):
        i = 0
        for particle in particles:
            results[iteration][i] = particle.move(particles)
            i += 1
    
    # Create csv formatted string for output
    result_output = ''
    for iteration in results:
        for particle in iteration:
            result_output += str(particle[0]) + ',' + str(particle[1]) + ','
        result_output += '\n'
    
    result_output = result_output.replace(',\n', '\n')
    
    print result_output
    

    output is:

    21.9,2.0,-18.9,10.1
    23.7,4.1,-17.7,20.1
    25.4,6.3,-16.4,30.0
    27.0,8.6,-15.0,39.8
    28.5,11.0,-13.5,49.5
    29.9,13.5,-11.9,59.1
    31.2,16.1,-10.2,68.6
    32.4,18.8,-8.4,78.0
    33.5,21.6,-6.5,87.3
    34.5,24.5,-4.5,96.5
    35.4,27.5,-2.4,105.6
    36.2,30.6,-0.2,114.6
    36.9,33.8,2.1,123.5
    37.5,37.1,4.5,132.3
    38.0,40.5,7.0,141.0
    38.4,44.0,9.6,149.6
    38.7,47.6,12.3,158.1
    38.9,51.3,15.1,166.5
    39.0,55.1,18.0,174.8
    39.0,59.0,21.0,183.0
    38.9,63.0,24.1,191.1
    38.7,67.1,27.3,199.1
    38.4,71.3,30.6,200
    38.0,75.6,34.0,192.0
    37.5,80.0,37.5,183.9
    37.1,84.5,40.9,175.7
    36.8,89.1,44.2,167.4
    36.6,93.8,47.4,159.0
    36.5,98.6,50.5,150.5
    36.5,103.5,53.5,141.9
    36.6,108.5,56.4,133.2
    36.8,113.6,59.2,124.4
    37.1,118.8,61.9,115.5
    37.5,123.9,64.5,106.7
    38.0,128.9,67.0,98.0
    38.6,133.8,69.4,89.4
    39.3,138.6,71.7,80.9
    40.1,143.3,73.9,72.5
    41.0,147.9,76.0,64.2
    42.0,152.4,78.0,56.0
    43.1,156.8,79.9,47.9
    44.3,161.1,81.7,39.9
    45.6,165.3,83.4,32.0
    47.0,169.4,85.0,24.2
    48.5,173.4,86.5,16.5
    50.1,177.3,87.9,8.9
    51.8,181.1,89.2,1.4
    53.6,184.8,90.4,-6.0
    55.5,188.4,91.5,-13.3
    57.5,191.9,92.5,-20.5
    

    when it should be:

    21.9,2.0,-18.9,10.0
    23.7,4.1,-17.7,19.9
    25.4,6.3,-16.4,29.7
    27.0,8.6,-15.0,39.4
    28.5,11.0,-13.5,49.0
    29.9,13.5,-11.9,58.5
    31.2,16.1,-10.2,67.9
    32.4,18.8,-8.4,77.2
    33.5,21.6,-6.5,86.4
    34.5,24.5,-4.5,95.5
    35.4,27.5,-2.4,104.5
    36.2,30.6,-0.2,113.4
    36.9,33.8,2.1,122.2
    37.5,37.1,4.5,130.9
    38.0,40.5,7.0,139.5
    38.4,44.0,9.6,148.0
    38.7,47.6,12.3,156.4
    38.9,51.3,15.1,164.7
    39.0,55.1,18.0,172.9
    39.0,59.0,21.0,181.0
    38.9,63.0,24.1,189.0
    38.7,67.1,27.3,196.9
    38.4,71.3,30.6,200
    38.0,75.6,34.0,192.1
    37.5,80.0,37.5,184.1
    37.1,84.5,40.9,176.0
    36.8,89.1,44.2,167.8
    36.6,93.8,47.4,159.5
    36.5,98.6,50.5,151.1
    36.5,103.5,53.5,142.6
    36.6,108.5,56.4,134.0
    36.8,113.6,59.2,125.3
    37.1,118.8,61.9,116.5
    37.5,123.9,64.5,107.8
    38.0,128.9,67.0,99.2
    38.6,133.8,69.4,90.7
    39.3,138.6,71.7,82.3
    40.1,143.3,73.9,74.0
    41.0,147.9,76.0,65.8
    42.0,152.4,78.0,57.7
    43.1,156.8,79.9,49.7
    44.3,161.1,81.7,41.8
    45.6,165.3,83.4,34.0
    47.0,169.4,85.0,26.3
    48.5,173.4,86.5,18.7
    50.1,177.3,87.9,11.2
    51.8,181.1,89.2,3.8
    53.6,184.8,90.4,-3.5
    55.5,188.4,91.5,-10.7
    57.5,191.9,92.5,-17.8
    

    For some reason, which I am unable to locate, the floats in the last column are not quite right. They differ from the desired output by anything between 2 and 0.5 or so.

    I have no idea why this is the case!

    Thanks for any help!

    +6  A: 

    You need to make a snapshot of the state and perform your calculations on the snapshot. If you move the particles around during the calculations as you are currently doing, you will get inconsistent results.

    Something like this may work

    from copy import deepcopy
    for iteration in range(times):
        i = 0
        particles_snapshot = deepcopy(particles)
        for particle in particles:
            results[iteration][i] = particle.move(particles_snapshot)
            i += 1
    
    gnibbler
    AWESOME!It works!Thanks very much!
    Jasper