Take Home Exam 3

INTRODUCTION

In this THE, you will be doing some physics and have the chance to visualize your solutions graphically. Your task is to implement a simulation of a gravitational many-body problem (see the gure on the next page).

In this problem, you will be dealing with particles that have mass (denoted by m), which is an invariant property. In addition to mass, the particles have dynamic properties, characterized by two vectors: namely the position ~r, and the velocity ~v. In this THE, for the sake of simplicity, all vectors (position and velocity) will be two-dimensional.

To predict (i.e., compute) a particle system’s behavior in time, we have to know the masses and the initial values of the dynamic properties (i.e., positions and velocities) of the particles. Then, the laws of physics allow us to predict the future of the dynamics1.

Our system will consist of n particles. We will label our particles with integer subscripts. For example, m4 will denote the mass of the 4th particle. The values in the following table will be provided as input.

Particle
Mass
Initial position
Initial velocity




p1
m1
(x01; y01)
(v0x1; v0y1)
p2
m2
(x02; y02)
(v0x2; v0y2)
.
.
.
.
.
.
.
.
.
.
.
.
pn
mn
(x0n; y0n)
(v0xn; v0yn)

In the  gure on the next page, you see an example of a many-body particle system with 3 particles.


1We con ne ourself to classical mechanics, which is quite correct at the macro scale. If we go down to atomic sizes, then this approach will no more predict correct results and you must pick a quantum mechanical approach.
Updating the system

Consider a single particle: As you know, velocity means \change in position in unit time". Normally, we express this fact as a derivative, considering \unit-time" as an in nitely small number:

~v = d~rdt


If we discretize this (i.e., approximate the in nitely small dt by some ‘small’ number), then we write the above derivative as:




~v
=
~r







t






Y




M1: 17.5



R
: (-2,3.5)





O 1




R
: (2,1)
V
: (0.5,-1.5)

O2
O1




V
: (-1,0)






O2







M2: 7.25







X



M3: 42.4




R

: (-1.5,-1)



O 3





V

:
(1,0)




O







3









The di erence now is that t is just a ‘small’ number. Then, at any time t, if a particle is at position ~rt, then after the unit-time update, the new position ~rt+ t will be:

~rt+  t = ~rt +  t ~v
(1)

Equation 1 always holds, even if the velocity is changing.

If the particle is under any force, then its velocity will change. So, we should also update the velocity ~v after unit-time has passed. The change in the velocity is called acceleration and is denoted by the vector ~a. Newton’s second law provides a relation that gives us the acceleration at any time the particle will attain, if we know the force acting on the particle at time t as well as its mass:
~
~a = mF


By de nition, ~a is:

~a = d~vdt


Again, if we discretize the time and consider a ‘small’ unit-time, a velocity ~v will be updated according to:


~

~vt+  t = ~vt +  t
Ft
(2)

m



In order to calculate all the dynamics of a single particle, we have to know the forces acting
~
~
acting on a
on it. Force is superposable, which means that if there are two forces F1
and F2

~
particle, you can assume that these forces can be substituted by a single force Fsum, which equals

~
~
. Hence, generally speaking, if a particle is under k number of forces, the resulting force
F1
+ F2

is the vector sum of all of them:
k






~
~


Fsum =
Xi



Fi



=1

Now, what kind of force is exerted by a particle p0 on the particle p? In this THE, we will consider a gravitational force. As you now, the gravitational force is inverse quadratic proportional to the distance between the two particles (p and p0):
F~= G

mp

mp0
2
u^rr0


~rp0

~rp
k












k






where G is the so-called universal gravitational constant (a real number). Furthermore, u^rr0  is the

unit vector in the direction from particle p towards p0. In other words, u^rr0  = (~rp0    ~rpk

So, we could have written the force as:

F~= G


mp

mp0
3
(~rp0
~rp)



~rp0

~rp
k















k







If we wish to calculate the resulting force acting on a particle pi by all other particles making use of the superposition principle, we can write:


m


mi  mj



F~    =
G
k



(r~   r~ )
isum
Xj


r~j    r~i
k
3

ji

=1

















j6=i

Taking all constants out and rearranging the order of subtraction :


m









~
Xj



mj





Fisum = G mi


k
r~i
k
3
(r~i
r~j )
(3)

=1


r~j


















j6=i

With the equations (1), (2) and (3), it is possible to calculate new positions and velocities for all particles after a unit-time t has passed (Provided that the initial positions and velocities are given).


PROBLEM

You will be given the G constant, the unit-time t, and the mass and the initial (position and velocity) values for a set of particles as a list. This list will be provided to you by a function call of get_data(), which will be provided by the evaluators (i.e., us) at the time of grading.

You are expected to write a function that you will exactly name as new_move(), which will take no arguments and (at each call) return the ~r ( ~r = t ~v) values for each particle (for the t time interval that has passed). This returned list is simply the information of how much each particle will move in the x- and y-directions. The evaluators will be able to observe the whole event simulation by successive calls to new_move() .

What you are supposed to submit is the the3.py le holding the function new_move() (and its helpers). However, we provide you with a visual tool, so that you can visualize how your new_move() performs. This Python code, when loaded, attempts to load your Python code that lives in the3.py, makes the necessary initializations, opens up an 800 800 screen, places circles of sizes proportional to the masses of the particles and displays the action by internally calling your new_move() successively.
HOW-TO-USE draw.py

The use of draw.py is optional. It serves only visualization purposes. Your code the3.py will be evaluated using black-box batch testing and not visually.

Place the3.py, draw.py, evaluator.py in your working directory and in Python do a from draw import *

We provide sample the3.py and evaluator.py les with dummy new_move() and get_data() function de nitions for your convenience. The sample new_move() does no gravitational force calculation at all, and always returns a constant result; it is only for you to see how draw.py works.

draw.py uses a coordinate system where the origin is at the top-left corner. x-axis is in the left-to-right direction and y-axis is in the top-to-bottom direction.

There are two global variables used by draw.py. If you want to alter their default values, do so by editing their assignment lines in draw.py.

DELAY: The delay (time in milliseconds) between two successive new_move() calls. The default is 200 (ms).

SCALE: A oating point constant to determine which screen coordinate corresponds to which physical coordinate. The relation is simple:

hScreen coordinatei = hPhysical coordinatei    SCALE

If you want to have a zoom-out view, you can set SCALE < 1:0. The default is set to

be 1:0.

To avoid name clashes, do not re-de ne something with the names of the global variables and functions of draw.py.


SPECIFICATIONS

There will be at least 1 particle and at most 20 particles.

To get your input data, you shall execute a call to the function get_data() in the3.py. This call will return a list like:

[G,  t,[m1,x01,y01,v0x1,v0y1],[m2,x02,y02,v0x2,v0y2],,[mn,x0n,y0n,v0xn,v0yn]]


G is the gravitational constant. t is the ‘small’ unit time. The following lists provide the initial values for the particles (namely the mass, initial x and y coordinates, followed by the initial velocity’s x- and y-components). All numbers will be oating point values. Multiple calls to get_data() in a single run will return the same values.

In your the3.py, you are allowed to import only from the math module. You cannot import form any other le=library.

Evaluating the correctness of the updates calculated by your solution will require performing


oating point comparison (by our scripts). We will consider two values f1 and f2 the same if jf1 f2j < 0:001 is true.
To test your own code, you will need a get_data() function. For that purpose, you should construct your own get_data() function. We suggest that you write the code that de nes your get_data() into a le named evaluator.py, which will be automatically loaded by the draw.py package.

{ Do not de ne this get_data() (that you will be using for your own testing) in the the3.py le.

{ Do not submit your evaluator.py le or any other le except the3.py. As evaluators, we are not interested in them. We will implement our own get_data().

Make sure that new_move() is implemented in the3.py, and conforms to the speci cations. You are free to implement helper functions as much as you like. Moreover, you can (and probably will) use global variables. Do not execute a call of new_move() in the3.py since this is the job of the evaluators.

Each call to new_move() should return a list of the form

[[  x1,  y1],[  x2,  y2],,[  xn,  yn]]

where the xi, yi values are the coordinate increments for the ith particle’s position in a t time (certainly these values can be negative as well as positive).

Do not perform any error checks. The input will be error free.

You can de ne other functions or variables in the3.py as long as there are no name clashes with other functions and variables in the other les.

You can only use the concepts covered in the lectures until the deadline of THE3.


GRADING

Comply with the speci cations. Since your returned results will be evaluated automatically, non-compliant results will be considered as incorrect by our evaluation system.

Your program will be tested with multiple data (a distinct run for each data).

Any program that performs only 30% and below will enter a glass-box test (eye inspection by the grader TA). The TA will judge an overall THE2 grade in the range of [0,30]. The glass-box test grade is not open to negotiation nor explanation.

A program based on randomness will be graded zero.