function makeGalaxy (nP)
% function makeGalaxy ()

% gravitational constant
global Gc
Gc = 6.67;

rC = 500;
mP = 50;
vP = 30;

[M1, X1, V1] = makeCluster(rC, nP/2, mP, vP);
[M2, X2, V2] = makeCluster(rC, nP/2, mP, vP);

ext = ones(nP/2,1);

X1 = X1 + ext*[1000,0];
X2 = X2 - ext*[1000,0];

V1 = V1 + ext*[-10, 100];
V2 = V2 - ext*[-10, 100];

% save
galaxy = [M1 X1 V1; M2 X2 V2];
save 'galaxy.dat' galaxy -ascii -double;


%====================================================================
function [M, X, V] = makeCluster (rC, nP, mP, vP)
% function [M, X, V] = makeCluster (rC, nP, mP, vP)
%
% parameters:
%   rC = radius of cluster
%   nP = number of planets
%   mP = mass of planets
%   vP = velocity of planets

M = zeros(nP, 1);
X = zeros(nP, 2);
V = zeros(nP, 2);

for k = 1:nP
   r = rC * sqrt(rand);
   x = r  * exp(i*(2*pi*rand));
   v = vP * exp(i*(2*pi*rand));
         
   M(k  ) = mP;
   X(k,:) = [real(x) imag(x)];
   V(k,:) = [real(v) imag(v)];
end