% This is the Nesl code for calculating the pairwise interactions between
a set of n bodies, for a single timestep. It is a brute-force method that
calculates the interaction between every pair and takes in O(n^2) time.
%
% Calculate acceleration at (x0,y0,z0) due to particle (x1,y1,z1,m1) %
function force((x0,y0,z0),(x1,y1,z1,m1)) =
let
dx = x1-x0;
dy = y1-y0;
dz = z1-z0;
rsq = (dx^2 + dy^2 + dz^2);
mrcubed = m1/(rsq*sqrt(rsq))
in dx*mrcubed,dy*mrcubed,dz*mrcubed $
% functions to manipulate 3D vectors %
function unzip3(v) =
let
(xs,yzs) = unzip(v);
(ys,zs) = unzip(yzs);
in
(xs,ys,zs) $
function sumv(vec) =
let
(vecx,vecy,vecz) = unzip3(vec);
in
(sum(vecx),sum(vecy),sum(vecz)) $
% acceleration at a point due to a vector of particles %
function acc_p(((x0,y0,z0),m),v) =
let
v1 = {(x1,y1,z1,m1) : ((x1,y1,z1),m1) in v | (x1 /= x0) or (y1 /= y0
) or (z1 /= z0)};
forces = {force((x0,y0,z0),v1): v1};
in sumv(forces) $
% Acceleration for every particle due to all the rest %
% It takes as an argument the vector of particles. Each particle is
represented as ((x,y,z),mass).
It returns the vector acclerations on each particle as [(ax, ay, az)].
%
function get_dir_forces(v) = {acc_p(v1,v) : v1 in v} $