Making a Ring Galaxy

OK, in this assignment you are going to make a 2D simulation of the formation of a ring galaxy. One way to form a ring galaxy is to drop a small satellite galaxy down the rotation axis of a large disk galaxy. We are going to mock this up in 2D by simulating a disk galaxy where we slowly turn on and then turn off a satellite galaxy potential at the center of the disk galaxy.

Galaxy Potentials:

For the disk galaxy potential, use the logarithmic potential: phi = 0.5*v_c²*ln(r²). If we are simulating the Milky Way, what should the value of v_c be?
For the satellite galaxy potential, use the Hernquist potential, with mass M_tot=10¹⁰ M_sun and scale radius a=300 pc.

Step 1: Pushing a particle

Ignore the satellite galaxy stuff for now. Set a particle at the position x=4000 y=0. Give it the proper velocity vx and vy to make it go on a circular orbit.
Choose a reasonable timestep, and integrate the orbit using the midpoint method. Demonstrate that the particle travels on a circular orbit.
Tweak the velocities a bit so it no longer travels on a circular orbit, and demonstrate that the new orbit is reasonable.

Step 2: Initializing the disk

Now use the initialization code given in class to initialize a few thousand particles in an exponential disk with scale length 2.5 kpc and max radius 10 kpc
Assign the particles the proper velocities to move on circular orbits -- think carefully about the velocity vectors here!
Integrate the model forward for one rotation period at the solar radius (8.5 kpc) and demonstrate that the model is stable.

Step 3: Adding a satellite

Add the satellite galaxy potential to the galaxy model, but have the mass be a gaussian function of time: M(t) = M_tot * exp(-(t-t₀)²/(2*sig²)) where t₀=50 Myr and sig=10Myr. This mimics a satellite galaxy passing through the center of the galaxy.
Watch and describe what happens to the disk galaxy over ~ 500 Myr period.
Compare to pictures of ring galaxies and discuss.

Tips:

work in units where length is in parsecs, masses are in solar masses, velocities are km/s, and time is in Myr. In these units G=4.43x10^-3.
in your integration loop, use a timestep counter and ask for plots every Nth timestep (i.e., every 25th) like this:

do i=0,nstep {

define time $($i*$dtime)

(step model forward)

if (($i % 25)==0) {

(plot model)
echo Time=$time
define next ?

}