Lesson 1: Simulating orbits
x1(0), vx2(0), vy1(0), vy2(0) the solution of these equations x1(t), x2(t), y1(t), y2(t) and vx1(t), vx2(t), vy1(t), vy2(t) describes the time-evolution of the system in the position and velocity space respectively.The solution to the Kepler problemJohannes Kepler proved that in the limit when one body is much more massive that the other (e.g., the Sun and a planet) the orbits are ellipses with the more massive body at one focus. This is known as Kepler's first law. This ellipse shown in the diagram above has a semi-major axis a and eccentricity e and is described by the following equation:having an orbital period of: Solving for Orbital motion on a Computer: The Leapfrog methodIn the leapfrog method we discretize time as t_(n+1/2)=t_n+ dt/2. This means that we use the velocity at time t_(n+1/2) to find the position at time t_(n+1) and the velocity at time t_(n+3/2) as is shown in the following diagram. If we know the position at time t_n, x1(t_n) we compute the position at time t_(n+1), x1(t_(n+1)) using the following equation : However, we will use a simplified version of it which is (hereafter equation (1)): x2(t_(n+1)), y1(t_(n+1)) and y2(t_(n+1)). If we know the velocity at time t_(n+1/2), vx1(t_(n+1/2)) we compute the velocity at time t_(n+3/2), vx1(t_(n+3/2)) sing the following equation: Similarly we do with vx2(t_(n+3/2)), vy1(t_(n+3/2)) and vy2(t_(n+3/2)) as follows:In a simplified version this equation becomes (hereafter equation (2)):Numerical modeling in Python The leapfrog method is implemented in the function leapfrog_step() defined in hw1.py. Homework Assignment1) The initial positions and velocities are stored in an 8-element array called rv. Organize this array such that rv[0:2] contains the x and y position of body 1, rv[2:4] the x and y components of the velocity of body 1 and similarly for body 2 (i.e., correct the first four #FIX ME! statements) 2) Update the positions to their values at time t1 (next two #FIX ME! statements) using equation (1) (see figure below) but the velocity at time t0. the correct way would be to use the velocity at time n+1/2 but we ignore this subtlety in the current homework. 3) Similarly Compute velocities at time t1 using equation (2) (see figure below) (next two #FIX ME! statements) 4) Store the new values back into the array rv (last four #FIX ME! statements) Solving for Orbital motion on a Computer: Choosing Initial conditionswhere m is the total mass of the system . This is implemented in the the function r0v0 (). A convenient choice of initial positions and velocities is starting the system while at closest approach, i.e., at periapsis. The initial positions and velocities at periapsis relative to center of mass of the system are then simplified as: where q is the mass ratio defined as q=m1/m2. This is implemented in the main function of our program hw1() that we explain below. Solving for Orbital motion on a Computer: Running the program The main function of our program hw1() : takes as arguments the initial conditions of our problem, i.e., the masses of the bodies m1 and m2, the semi-major axis a, the eccentricity e as well as the maximum time until which the system is evolved. The time step for the function is set to a default value tstep=1.0e-3 so you can use tmax=nsteps*tstep where nsteps is the number of steps you want to evolve the system. Using different values for nsteps and trying different initial conditions run the program clicking on the green play button and calling in the terminal the function hw1(). For example, typing in the terminal sets m1=1,m2=3,a=1,e=0.5 and runs the program for tmax=1=1000*1.0e-3. i.e., for 1000 time steps which is equivalent to one orbit. The program shows in the plot produced the positions of bodies 1 and 2 in blue and green respectively as well as their relative position in red. The relative motion of the two bodies is an ellipse with the semi-major axis and eccentricity defined in the initial conditions verifying Kepler's first law. |