%
% Demo of Newtonian N-body problem 
% (c) 2006, Howard E. Motteler
%
% read initial data as nbody_init.dat
% writes final data as nbody_last.dat
% see "mkbod.m" for the file format
%
% script variables
%   M		 n x 1, object masses
%   Px, Py, Pz	 n x 1, object coordinates
%   Vx, Vy, Vz	 n x 1, object velocitys
%   D		 n x n, distance between objects
%   Fx, Fy, Fz	 n x n, forces between objects
%   G		 gravitational constant (scalar)
%   dT		 time increment (scalar)

% call mkbod.m to generate a new set of bodies
mkbod

wtype = 'ieee-be';   % for IBM Power PC
% wtype = 'ieee-le';     % for Intel/AMD

% defaults
% k  = 1;
% dT = 0.01;
G  = 1;
k  = input('number of iterations > ');
dT = input('enter dT > ');

% read the initial dataset
fp = fopen('nbod_init.dat', 'r', wtype);
m = fread(fp, 1, 'int32');  
n = fread(fp, 1, 'int32');
tmp = fread(fp, 7 * n, 'float32');
fclose(fp);
fprintf(1, 'read data for %d bodies\n', n);

% assing names to fields
tmp = reshape(tmp, 7, n)';
M  = tmp(:,1);
Px = tmp(:,2);
Py = tmp(:,3);
Pz = tmp(:,4);
Vx = tmp(:,5);
Vy = tmp(:,6);
Vz = tmp(:,7);

% set initial projection axes
v = 2 * [min(Px),max(Px),min(Py),max(Py)];

% loop on dT steps
tic
for i = 1:k

  % distance between all pairs of objects
  Fx =  ones(n,1) * Px' - Px * ones(1,n);
  Fy =  ones(n,1) * Py' - Py * ones(1,n);
  Fz =  ones(n,1) * Pz' - Pz * ones(1,n);
  
  Dsq = Fx .* Fx + Fy .* Fy + Fz .* Fz + eye(n);
  D = sqrt(Dsq);
  
  % mutual forces between all pairs of objects
  F =  G * (M * M') ./ Dsq;
  F =  F - diag(diag(F));   %  set 'self attraction' to 0
  Fx = (Fx ./ D) .* F ;
  Fy = (Fy ./ D) .* F ;
  Fz = (Fz ./ D) .* F ;
  
  % net force on each body
  Fnet_x = Fx * ones(n,1);
  Fnet_y = Fy * ones(n,1);
  Fnet_z = Fz * ones(n,1);
  
  % change in velocity:
  dVx = Fnet_x * dT ./ M;
  dVy = Fnet_y * dT ./ M;
  dVz = Fnet_z * dT ./ M;
  
  Vx = Vx + dVx;
  Vy = Vy + dVy;
  Vz = Vz + dVz;
  
  % change in position
  Px = Px + Vx * dT;
  Py = Py + Vy * dT;
  Pz = Pz + Vz * dT;
  
  plot(Px,Py,'o')
  axis(v)
  pause (.1);
  
end 

% concatenate data
tmp = [M,Px,Py,Pz,Vx,Vy,Vz];
tmp = tmp';

% write the dataset
fp = fopen('nbod_last.dat', 'w', wtype);
fwrite(fp, 7, 'int32');
fwrite(fp, n, 'int32');
fwrite(fp, tmp, 'float32');
fclose(fp);
fprintf(1, 'wrote data for %d bodies\n', n);
