function [mass,xyz,vxyz] = gen_solar(n)
% GEN_SOLAR : Generate dataset for nbody simulation
%
% [mass,xyz,vxyz] = gen_solar(n);
%
% Output is a three arrays of initial values for n-body simulation.
%    mass : n-by-1 array of masses in kg
%    xyz  : n-by-3 array of positions in m
%    vxyz : n-by-3 array of velocities in m/sec
%
% If n is omitted, generate the 10-body solar system.
% For 0 < n < 10, generate only the sun and the first n-1 planets.
% 
%
% John R. Gilbert  18 Jan 2010, last updated 12 Jun 2015

% sun mercury venus earth mars jupiter saturn neptune uranus pluto

if nargin == 0
    n = 10;
end;

earth = 4;
earthspeed = 29800;  % m/s

mass = [1.99e30; 3.30e23; 4.87e24; 5.97e24; 6.42e23; ...
    1.90e27; 5.69e26; 1.02e26; 8.68e25; 1.31e22];
dist = [0; 5.8e10; 1.08e11; 1.50e11; 2.28e11; ...
    7.78e11; 1.43e12; 2.87e12; 4.50e12; 5.91e12];

theta = 2*pi*rand(10,1); % start at random position in 2D orbit
xyz = [dist.*cos(theta) dist.*sin(theta) zeros(10,1)];

speed = earthspeed * (dist(earth) ./ dist).^.5;
speed(1) = 0;
direction = [-xyz(:,2) xyz(:,1) zeros(10,1)];
direction = diag(1./sqrt(sum(direction.^2,2))) * direction;
direction(1,:) = [0 0 0];
vxyz = diag(speed) * direction;
vxyz(1,:) = [0 0 0];

mass = mass(1:n, :);
xyz  =  xyz(1:n, :);
vxyz = vxyz(1:n, :);

