# Particle rotation in galaxy - runge kutta 2nd order help

## 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

Code:
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
Code:
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