C N-body calculation using the modified Euler method,
C which is globally 2nd-order.
C Uses a ring organization to calculate accelerations.
C
C Data arrays (position, velocity, acceleration, mass)
C are organized:
C
C 1 Particle number (i) --> N
C +--------------------------------------------+
C | x | |
C (c) y | |
C v z | |
C +--------------------------------------------+
C
C The parameter N and common blocks /NBODYFE/, /NBODYPE/ with arrays
C
C X0, X1 - Current, next positions
C V0, V1 - velocities
C M, GM - mass, G*mass (redundant)
C c, i - array coordinates
C and
C NB - integer <= N, actual no. of bodies
C G - real, the gravitational constant
C
C are defined in the file nbmode.h
C
C File: nbmode.fcm Header: nbmode.h
C------- setup() ---------------------------------------
subroutine setup()
INCLUDE 'nbmode.h'
C Initialize coordinates
c = SPREAD( [1:3], 2, N)
i = SPREAD([1:N], 1, 3)
end subroutine setup
C------- init() ----------------------------------------
C Load initial positions, velocities, also masses,
C gravitational const. All are read from logical
C unit 10. Replicates masses over x,y,z.
subroutine init()
INCLUDE 'nbmode.h'
integer k
real mx
X0 = 0.0 ! For unused particles
V0 = 0.0
do, k = 1, NB
read (10, *) X0(1, k), X0(2, k), X0(3, k)
end do
do, k = 1, NB
read (10, *) V0(1, k), V0(2, k), V0(3, k)
end do
M = 0.0 ! For unused particles
do, k = 1, NB
read (10, *) mx
M(:, k) = mx
end do
read (10, *) G
end subroutine init
C------- mkast(nx, ar, am) -----------------------------
C Insert nx asteroids at radius ar with total mass am.
C This assumes that an orbit of radius 1 about origin
C corresponds to speed 2*pi (e.g., au-earth-yr units).
subroutine mkast(nx, ar, am)
integer nx ! Number of new asteroids
real ar, am ! Orbit radius, total mass
INCLUDE 'nbmode.h'
real, parameter:: TWOPI = 6.2831853
real, array(1:3, 1:N) :: sx, cx, arg
real spd
integer nn
print 21, nx, ar, am
21 FORMAT (" Making ", I4, " asteroids, radius ",
* F10.5, ", mass ", F8.5)
nn = NB + nx
spd = TWOPI/sqrt(ar)
where ( (NB < i) .and. (i <= nn) )
arg = (TWOPI*(i - NB))/nx
sx = sin(arg)
cx = cos(arg)
end where
where ( (NB < i) .and. (i <= nn) .and. (c == 1) )
X0 = cx
V0 = -sx
end where
where ( (NB < i) .and. (i <= nn) .and. (c == 2) )
X0 = sx
V0 = cx
end where
where ( (NB < i) .and. (i <= nn) .and. (c == 3) )
X0 = 0.0
V0 = 0.0
end where
where ( (NB < i) .and. (i <= nn) )
X0 = X0*ar
V0 = V0*spd
M = am/nx
end where
NB = nn
end subroutine mkast
C------- print_state(PX) -------------------------------
C Print current state of system according to PX.
C See comments for main program for PX settings.
subroutine print_state(PX)
integer PX ! Print control
INCLUDE 'nbmode.h'
integer k, np
np = NB
if ( (PX >= 0) .and. (PX < NB) ) np = PX
print 21
21 FORMAT (/ " x y z ",
* " vx vy vz m")
do, k = 1, np
print 22, k, X0(1,k), X0(2,k), X0(3,k),
* V0(1,k), V0(2,k), V0(3,k), M(1,k)
end do
22 FORMAT (1X, I3, 6F10.5, E10.2)
end subroutine print_state
C------- print_x(ax, na, nf) -------------------------------
C Print ax(:, na:nf), nf-na <= 9 (for testing)
subroutine print_x(ax, na, nf)
INCLUDE 'nbmode.h'
real, array(1:3, 1:N) :: ax
CMPF ONDPU ax
CMF$ LAYOUT ax(,)
integer na, nf
integer u, lm
lm = min(nf, na+9)
do, u = 1, 3
print 21, ax(u, na:lm)
21 FORMAT ( 10(1X : E9.2 :) )
end do
end subroutine print_x
C------- print_c(Xf, Xc, Af, Ac, na, nb) -----------------------
C Print conveyor status in positions na:nb (for testing)
subroutine print_c(Xf, Xc, Af, Ac, na, nf)
INCLUDE 'nbmode.h'
real, array(1:3, 1:N) :: Xf, Xc, Af, Ac
integer na, nf
integer kk
print 22, (kk, kk = na, nf)
22 FORMAT ( 10(1X : I9 :) )
print *, " - Xf ----------"
call print_x(Xf, na, nf)
print *, " - Xc ----------"
call print_x(Xc, na, nf)
print *, " - Af ----------"
call print_x(Af, na, nf)
print *, " - Ac ----------"
call print_x(Ac, na, nf)
end subroutine print_c
C------- Grav(Xf) --------------------------------------
C Calculate accelerations (Newtonian gravitation)
C Uses systolic ring calculation, taking advantage
C of anti-symmetric force. Data being circulated
C on `conveyor' are marked `c', fixed items are `f'.
function Grav(Xf)
INCLUDE 'nbmode.h'
real, array(1:3, 1:N) :: Xf, Grav
CMPF ONDPU Xf
CMF$ LAYOUT Xf(,)
real, array(1:3, 1:N) :: Af, Ac, D, R, Xc, Mc
integer k, lm
Af = 0.0 ! Fixed accelerations
Ac = 0.0 ! Circulating accelerations
Xc = Xf ! Circulating positions
Mc = GM ! Circulating scaled masses
lm = min(NB, N/2) - 1 ! For loop control
do, k = 1, lm ! 2-way accel. steps
Xc = CSHIFT(Xc, SHIFT=-1, DIM=2) ! Step conveyor
Mc = CSHIFT(Mc, SHIFT=-1, DIM=2)
if ( k > 1 ) then
Ac = CSHIFT(Ac, SHIFT=-1, DIM=2)
end if
C Xc(:,i) = Xf(:,i-k), Similarly Mc and MG
C Ac(:,i) contributes to particle i-k
D = Xc - Xf
R = sqrt(SPREAD(SUM(D*D, 1), 1, 3))
WHERE ( GM > 0.0 .and. Mc > 0.0 )
D = D/R**3
Af = Af + Mc*D
Ac = Ac - GM*D
end WHERE
end do ! for 2-way accel.
C Ac(i) corresponds to Af(i-lm)
if ( NB > N/2 ) then ! Final 1-way accel.
Xc = CSHIFT(Xc, SHIFT=-1, DIM=2) ! Step conveyor
Mc = CSHIFT(Mc, SHIFT=-1, DIM=2)
Ac = CSHIFT(Ac, SHIFT=-1, DIM=2)
lm = N/2 ! lm gives shift of Ac
D = Xc - Xf
R = sqrt(SPREAD(SUM(D*D, 1), 1, 3))
WHERE ( GM > 0.0 .and. Mc > 0.0 )
D = D/R**3
Af = Af + Mc*D
end WHERE
end if ! for final accel.
D call print_c(Xf, Xc, Af, Ac, 29, 35)
Ac = CSHIFT(Ac, SHIFT=lm, DIM=2) ! Back to correct position
D call print_c(Xf, Xc, Af, Ac, 1, 7)
C Combine accels for final result.
Grav = Af + Ac
end function Grav
C------- modeul(h, ns) --------------------------------
C Perform ns steps of modified Euler with time-step h.
C X0, V0 are taken to be current state, and are updated
C with final state. X1, V1 are used in the process.
subroutine modeul(h, ns)
real h
integer ns
INCLUDE 'nbmode.h'
INTERFACE
function Grav(Xf)
INCLUDE 'nbmode.h'
real, array(1:3, 1:N) :: Xf, Grav
CMPF ONDPU Xf
CMF$ LAYOUT Xf(,)
END INTERFACE
real h2
integer k
h2 = 0.5*h ! h/2, for convenience
do k = 1, ns ! ns steps
X1 = X0 + h2*V0 ! Half-way positions
V1 = V0 + h2*Grav(X0) ! Half-way velocities
X0 = X0 + h*V1 ! Next positions
V0 = V0 + h*Grav(X1) ! Next velocities
end do
end subroutine modeul
C------- nbmode -------------------------------
C Main program: reads problem and control information
C from data file named by user, performs and times
C calculation. Print control:
C 0 - just performance data
C < 0 - all
C > 0 - number given
program nbmode
INCLUDE 'nbmode.h'
integer ns, np, nv, k, PX
real dt, ar, am
character*50 fnm
print '(" Data file> "$)' ! Request data file
read *, fnm
print '(" Print ctl> "$)' ! Request print control
read *, PX
C Begin initialization
open(10, file=fnm, status = "OLD")
read (10, *) NB
call setup() ! Setup geometric structure
call init() ! Read data for individual bodies
C Read "asteroid" descriptions from file
read (10, *) ns
do while ( ns > 0 )
read (10, *) ar, am
call mkast(ns, ar, am)
read (10, *) ns
end do
read (10, *) ns, np, dt ! Integration parameters
close(10)
print 21, NB, N
21 FORMAT (/, X, I4, " bodies, geometry allows ", I4)
print 22, G, dt
22 FORMAT (" G = ", E13.6, " delta-t = ", E12.5)
print *, " simulated time = 0"
call print_state(PX)
C Scale masses
GM = G*m
C np runs of ns steps, timestep dt
do, k = 1, np
call CM_timer_clear(1)
call CM_timer_start(1)
call modeul(dt, ns)
call CM_timer_stop(1)
call CM_timer_print(1)
print 23, k*ns*dt
23 FORMAT (/ " simulated time = ", F9.3)
call print_state(PX)
end do
end program nbmode
C*******************************************************
C Data files:
C
C N No. of explicit bodies
C
C x y z Position of body 1
C : :
C x y z Position of body N
C
C vx vy va Velocity of body 1
C : :
C vx vy va Velocity of body N
C
C m Mass of body 1
C : :
C m Mass of body N
C
C G Gravitational constant
C
C nx Number,
C ar am radius, mass of asteroid group
C : : 0 or more groups
C nx Number,
C ar am radius, mass of asteroid group
C -1 Stopper
C
C ns Number of time steps per phase
C np Number of phases
C dt delta-t (per time step)
C
C Further data in the file will be ignored.
C*******************************************************