Particle rotation in galaxy - runge kutta 2nd order help

[InuYasha]

Homework Statement

Hi everyone,
I'm a bit stuck and have been staring at my fortran code all day. For my project this year I'm writing a program to model a spiral galaxy in which gas clouds circulate, collide and produce star formation. I've written code for gravitational acceleration due to the disk, halo, bulge and spiral arms. So I have gaccel(x,y,z). I was then able to produce a plot showing the correct galaxy rotation curve using the code below. p_TOT = total gas clouds, factor was an adjustment to get km/s

I now want to make the clouds move using a runge kutta 2nd order method like below. I wrote an RK2 program to solve an ODE as practice and it worked ok (though the method is slightly different from the one I've been asked to implement). I've got a mental block on how to implement RK2 to make the clouds go round the galaxy. Am I being a bit dumb and missing something obvious :-/
Any help would be really appreciated.
Thanks.

Homework Equations

a(n) = a(r(n))
r(n+1/2) = r(n) + v(n)*tstep/2
v(n+1/2) = v(n) + a*tstep/2
a(n+1/2) = a(r(n+1/2)
r(n+1) = r(n) + v(n+1/2)*tstep
v(n+1) = v(n) +a(n+1/2)*tstep
t(n+1) = t(n) +tstep

The Attempt at a Solution

DO p=1,p_TOT                                      ! loop over particles
    r_MAG=SQRT(r(1,p)**2+r(2,p)**2+r(3,p)**2)     ! compute distance from centre
    a_MAG=SQRT(a(1,p)**2+a(2,p)**2+a(3,p)**2)     ! compute magnitude of acceleration
    v(1,p)=+a(2,p)                                ! compute x-component of velocity
    v(2,p)=-a(1,p)                                ! compute y-component of velocity
    v(3,p)=0.                                     ! set z-component of velocity to zero
    v(1:3,p)=v(1:3,p)*SQRT(r_MAG/a_MAG)           ! scale velocity to centrifugal balance
    v_MAG=SQRT(v(1,p)**2+v(2,p)**2+v(3,p)**2)     ! compute magnitude of velocity
    x(p)=r_MAG                                    ! compute abscissa of plotted point
    y(p)=factor*v_MAG                             ! compute ordinate of plotted point
ENDDO                                             ! end loop over particles

PROGRAM RK2
IMPLICIT NONE

REAL(kind=8) :: t=0.0,h=0.1
INTEGER :: i, tmax=30
REAL(kind=8),DIMENSION(30) :: y
y(1)=5
DO i=1, tmax
    CALL RUNGEKUTTA(t,y,h,i)
    t=i*h
    WRITE (*,*) "y=",y(i+1)
ENDDO

END PROGRAM RK2

SUBROUTINE RUNGEKUTTA(t,y,h,i)
IMPLICIT NONE
INTEGER :: i
REAL(kind=8) :: t,k1,k2,h
REAL(kind=8),DIMENSION(30) :: y
REAL :: EQN
k1=EQN(t,y,i)
k2=EQN(t+h,y+k1*h,i)
y(i+1)=y(i)+0.5*(k1+k2)*h
END SUBROUTINE RUNGEKUTTA

REAL FUNCTION EQN(t,y,i)
IMPLICIT NONE
INTEGER :: i
REAL(kind=8) :: t
REAL(kind=8),DIMENSION(30) :: y
EQN=3*exp(-t)-0.4*y(i)
RETURN
END FUNCTION EQN

 

[InuYasha]

I think what I'm asking is what would the likely function be that I'm integrating?
I'm pretty sure the code regarding scaling the velocity to centrifugal force only gets called once. Would I need to call the subroutine calculating the gravitational accel (I haven't posted this as it's pretty long. It just gives an acceleration value for x,y,z at the end though). That isn't a differential equation though.
It doesn't help that we haven't been exposed to much programming apart from Mathcad and a little Matlab and now in 3rd year physics we're expected to code everything in languages we don't know. Ay carumba.