C N-body calculation using the modified Euler method.
C Uses O(N^2) processors to calculate accelerations.
C
C Data arrays (position, velocity, acceleration, mass)
C are (N, N, 3), organized:
C
C ________________________________
C /| z /_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/|
C (c) y /_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/ |
C / x /_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/_/ /|
C 1 |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C | |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C | |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C P |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C a |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C r |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C t |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C i |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C c |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C l |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C e |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C (i) |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//|
C | |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|///
C v |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|//
C N |_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|/
C
C 1 --- Particle (j) --> N
C
C The parameter N and common block /NBODYPE/
C with arrays
C
C X0, X1 - Current, next positions
C V0, V1 - velocities
C MI,GJ - mass(i), G*mass(j)
C i,j,c - array coordinates
C diag - diagonal entries
C act - active in acceleration
C
C and common block /NBODYFE/ with
C
C NB - integer <= N, actual no. of bodies
C G - real, the gravitational constant
C
C are defined in the file nb2mode.h
C
C System wide conventions regarding data:
C
C Entries (i, :, :) - data for particle i
C except GJ (:, j, :) = G*mass(j)
C
C This is intentionally redundant, and the programs
C perform many redundant calculations.
C
C File: nb2mode.fcm Header: nb2mode.h
C------- setup() ---------------------------------------
C Initialize coordinates, logical flags.
subroutine setup()
INCLUDE 'nb2mode.h'
i = SPREAD(SPREAD([1:N], 2, N), 3, 3)
j = SPREAD(SPREAD([1:N], 1, N), 3, 3)
c = SPREAD(SPREAD([1:3], 1, N), 2, N)
diag = (i == j)
actv = .FALSE.
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.
C Data for particle i are left in (i, :, :) .
subroutine init()
INCLUDE 'nb2mode.h'
integer k
real mx
X0 = 0.0 ! For unused particles
V0 = 0.0
do, k = 1, NB
read (10, *) X0(k, 1, 1), X0(k, 1, 2), X0(k, 1, 3)
end do
do, k = 1, NB
read (10, *) V0(k, 1, 1), V0(k, 1, 2), V0(k, 1, 3)
end do
M = 0.0 ! For unused particles
do, k = 1, NB
read (10, *) mx
MI(k, 1, :) = mx
end do
read (10, *) G
X0 = SPREAD(X0(:, 1, :), 2, N)
V0 = SPREAD(V0(:, 1, :), 2, N)
MI = SPREAD(MI(:, 1, :), 2, N)
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).
C Data for particle i are left in (i, :, :) .
subroutine mkast(nx, ar, am)
integer nx ! Number of new asteroids
real ar, am ! Orbit radius, total mass
INCLUDE 'nb2mode.h'
real, parameter:: TWOPI = 6.2831853
real, array(1:N, 1:N, 1:3) :: 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
MI = 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 'nb2mode.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(k,1,1), X0(k,1,2), X0(k,1,3),
* V0(k,1,1), V0(k,1,2), V0(k,1,3), MI(k,1,1)
end do
22 FORMAT (1X, I3, 6F10.5, E10.2)
end subroutine print_state
C------- print_x(ax, nx, cx) -------------------------------
C Print ax(1:nx, 1:nx, cx), nx <= 10 (for testing)
subroutine print_x(ax, nx, cx)
INCLUDE 'nb2mode.h'
real, array(1:N, 1:N, 1:3) :: ax
CMPF ONDPU ax
CMF$ LAYOUT ax(,,)
integer nx, cx
integer u, lm
lm = min(nx, 10)
do, u = 1, lm
print 21, ax(u, 1:lm, cx)
21 FORMAT ( 10(1X : E9.2 :) )
end do
end subroutine print_x
C------- Grav(Xf) --------------------------------------
C Calculate accelerations (Newtonian gravitation)
C Uses completely parallel calculation, ignoring
C anti-symmetry of force. Acceleration for body i is
C left in (i, :, :) of result.
C Acceleration on body i due to body j is calculated
C in entries (i, j, :) .
function Grav(Xf)
INCLUDE 'nb2mode.h'
real, array(1:N, 1:N, 1:3) :: Xf, Grav
CMPF ONDPU Xf
CMF$ LAYOUT Xf(,,)
real, array(1:N, 1:N, 1:3) :: Ac, D, R, Xj
integer k, lm
C Spread Xf across j from diagonal:
where ( diag )
Xj = Xf
elsewhere
Xj = -1E30 ! Enormously negative
end where
Xj = SPREAD(MAXVAL(Xj, 1), 1, N)
C Now Xj(i,j,:) = X0(j,j,:)
Ac = 0.0
where ( actv )
D = Xj - Xf ! Displacement
elsewhere
D = 1.0 ! Protects from div by 0
endwhere
R = sqrt(SPREAD(SUM(D*D, 3), 3, 3))
where ( actv ) Ac = GJ*D/R**3
Grav = SPREAD(SUM(Ac, 2), 2, N)
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.
C Note that X's, V's are calculated redundantly.
subroutine modeul(h, ns)
real h
integer ns
INCLUDE 'nb2mode.h'
INTERFACE
function Grav(Xf)
INCLUDE 'nb2mode.h'
real, array(1:N, 1:N, 1:3) :: 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------- nb2mode -------------------------------
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 nb2mode
INCLUDE 'nb2mode.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 and distribute, set actv
where ( diag )
GJ = G*MI
elsewhere
GJ = -1E30 ! Enormously negative
end where
GJ = SPREAD(MAXVAL(GJ, 1), 1, N)
actv = (i <= NB) .and. (j <= NB) .and. .not. diag
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 nb2mode
C*******************************************************
C Data files: See nbmode.f
C*******************************************************