# Particle rotation in galaxy - runge kutta 2nd order help (1 Viewer)

### Users Who Are Viewing This Thread (Users: 0, Guests: 1)

#### InuYasha

1. The problem statement, all variables and given/known data
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.

2. Relevant 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

3. 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
1. The problem statement, all variables and given/known data

2. Relevant equations

3. The attempt at a solution

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

### The Physics Forums Way

We Value Quality
• Topics based on mainstream science
• Proper English grammar and spelling
We Value Civility
• Positive and compassionate attitudes
• Patience while debating
We Value Productivity
• Disciplined to remain on-topic
• Recognition of own weaknesses
• Solo and co-op problem solving