Misc notes:

1. On the initial condition generator for the disk galaxy model,
don't bother adding velocity dispersion to the particles (unless
you are really keen to...). Just make sure they are all on circular
orbits. All the relevant physics of ring galaxy making will still
be captured...

2. to use the runge kutta integrator, first make a "deriv" subroutine
that calculates wdot given w, ie (in fortran)

      subroutine deriv(t,w,wdot)

      real t,w(6),wdot(6)

      wdot[1]=w[4];
      wdot[2]=w[5];
      wdot[3]=w[6];
      wdot[4]=(calculate x acceleration here);
      wdot[5]=(calculate y acceleration here);
      wdot[6]=(calculate z acceleration here);

      return
      end

then to integrate a particle orbit, make an phase space array of 
x,y,z,vx,vy,vz (i.e. say w(1)=x, w(2)=y, w(3)=z, w(4)=vx, w(5)=vy,
w(6)=vz), then call the RK4 integrator to integrate w. 

in fortran, i would say

      RK4(w,wdot,6,tnow,dtime,wnew,deriv)

and do this for each particle at each time step.

3. to do more than on particle, store x,y,z,vx,vy,vz as arrays. then
loop over time and particle to integrate the model.

conceptually:

   do istep=1,nstep
     do ipart=1,npart
       (load w, wdot)
       RK4(w,wdot,6,tnow,dtime,wnew,deriv)
       (load wnew back into x,y,z, etc)
     enddo
   enddo


4. i have a fortran code on the tools page that calculates surface density
and velocity as a function of radius for a distribution of particles (ie your
disk model). It assumes your data is stored, one particle per line, as:

x1 y1 z1 vx1 vy1 vz1
x2 y2 z2 vx2 vy2 vz2
...
xN yN zN vxN vyN vzN