clear all

global M G N coll_distance a

N=100;
G=6.67408*10^-11; %m3 kg-1 s-2
M=zeros(N,1);

c_pluto=10^22; %sets amount to vary masses by
c_jupiter=10^27;
c_asteroid=.5*10^15;

M_earth=5.972*10^24; %kg
M_pluto=1.30900*10^22; %kilograms
M_jupiter=1.89813*10^27; %kg 
M_asteroid=10^15;  %kg


%initial conditions for 5 earths, 2 jupiters, 92 asteroids (N=100)
M(1)=1.989*10^30; %kg- mass of our sun
M(2:6)=M_earth -c_earth + 2*c_earth*rand(size(M(2:6),1),1);%mass of earths plus some variation
M(7:8)=M_jupiter - c_jupiter + 2*c_jupiter*rand(size(M(7:8),1),1); %mass of jupiters plus some variation
M(9:N)=M_asteroid - c_asteroid +2*c_asteroid*rand(size(M(9:N),1),1); %mass of asteroids plus variations

%initial conditions for mix of jupiters plutos earths (use N=50)
%M(2:N/2-10)=M_earth -c_earth + 2*c_earth*rand(size(M(2:N/2-10),1),1);%mass of earth plus some variation
%M(N/2+1-10:N-2)=M_pluto -c_pluto+ 2*c_pluto*rand(size(M(N/2+1-10:N-2),1),1); %masss of pluto plust some variation
%M(N-2:N)=M_jupiter - c_jupiter + 2*c_jupiter*rand(size(M(N-2:N),1),1); %mass of jupiter plus some variation

%initial conditions for half earths half plutos (use n=50)
%M(2:N/2)=M_earth/5 -c_earth + 2*c_earth*rand(size(M(2:N/2),1),1);%mass of earth plus some variation
%M(N/2+1:N)=M_pluto/5 -c_pluto+ 2*c_pluto*rand(size(M(1+N/2:N),1),1); %masss of pluto plust some variation


AU=149597870700; %meters
coll_distance=0.01*AU; %distance for inellastic collisions  (I don't call the function in the code that utilizes the variable
                       %because it slows things down, and doesn't affect
                       %outcomes all that much

a=AU*5; %5 AU in meters, size of box that particles will be distributed in
b=10000; % 10km/s  orbital velocity in range of initial velocities that particles will begin with

theta=360*rand(N,1);  %used to set initial positions and velocities
r=a*rand(N,1)+2*AU;
omega=(5.0*10^-8)*rand(N,1); %rad/s


%set initial positions and velocities
X=r.*cosd(theta);
Y=r.*sind(theta);
Z=-r/2+r.*rand(N,1);
X(1)=0; %put sun at origin
Y(1)=0;
Z(1)=0;
Vx=zeros(N,1);
Vy=zeros(N,1);
for i=1:N
    if theta(i)<=90 | 270 < theta(i) <= 180
        Vx(i)=r(i)*omega(i)*sind(theta(i));
        Vy(i)=-r(i)*omega(i)*cosd(theta(i));
    else
        Vx(i)=-r(i)*omega(i)*sind(theta(i));
        Vy(i)=r(i)*omega(i)*cosd(theta(i));
    end
end
Vz=zeros(N,1);


R=[X Y Z];
V=[Vx Vy Vz];
V(1,:)=[0 0 0]; %set sun velocity to zero
YY=[R V]  

store_YY(1:N,1:6,1)=YY; %stores all the entries in YY for each iteration for plotting later
t=0;
dt=24*60*60*4; %timestep is 4 days

clearvars c_earth c_pluto theta omega r  b c AU R V

for i=1:20000
    i
    
    %YY=inelastic_collisions3D(YY); %commented out because it slows down
    %the code and doesn't affect the outcome of the simulation very much
    
    
    F1= Nbody_ode3D(t      , YY             ); %runge-kutta method
    F2= Nbody_ode3D(t+dt/2 , YY + F1*(dt/2) );
    F3= Nbody_ode3D(t+dt/2 , YY + F2*(dt/2) );
    F4= Nbody_ode3D(t+dt   , YY + F3*dt     );
    YY = YY + (F1 + F2.*2 + F3.*2 + F4).*(dt/6);
    
    store_YY(1:N,1:6,i+1)=YY;
    clearvars F1 F2 F3 F4
    t=t+dt;

end

save('earthsjupitersasteroids1.mat') %save variables to file to plot later



