% %{
% Code adapted from Michael Hart's "Projectile Motion":
% https://www.youtube.com/watch?v=TrUa2gP1Tc8
%
% Additional aid for this particular simulation provided by Stephanie
%
% This program serves as a crude simulation of a random dispersion of
% projectiles from a hypothetical collision on a surface (this code was 
% used in the presentation to show simultaneously launched projectiles). 
%
% Some notes:
% - The kinematic equations for 3D are appropriately assigned for
%  the x, y, and z positions.
% - v0, theta, and phi are randomly assigned integers using randi(); the
% programmer can change the max as he/she prefers.
% - The function to enter, 'project_collision2(x0,y0,z0)' requires that x0,
% y0, and z0 are greater than 0; inputting 'project_collision(0,0,0)' makes
% sense conceptually and suffices.
% - The programmer can also change the gravitational constant (i.e. to the
%  Moon's or Mars') to see the resulting difference in projectile motions.
% %}

function [x,y,z] = project_collision2(x0,y0,z0)

%close all;

nParticles=50;     % Launched rocks we choose to use in simulation
                   % from "collision"
f=0;               % In place for saveas function later

% --- Iterate Entire Code for However Many "Launched Rocks" We Choose --- %
if x0>=0 && y0>=0 && z0>=0
    a=9.8;                      % Acceleration due to good ol' gravity
                                % a = 1.6 for the moon; a = 25 for Jupiter
    t=0:0.05:60;                % Time in seconds
    for i=1:nParticles
        v0=randi(50);           % Randomly selects a velocity up to 50 m/s
        theta=randi(90);        % Randomly selects an angle, "theta,"
                                % swept along the x-y plane.
        phi=randi(90);          % Randomly selects a trajectory angle,
                                % "phi," in the z-direction
                                
        % ----- Appropriate Kinematic Formulae ----- %
        % save for each ball
        x(i,:) = v0*t*cosd(theta) + x0;
        y(i,:) = v0*t*sind(theta) + y0;
        z(i,:) = v0*t*sind(phi) + -0.5*a.*t.^2;
    end
    
    % filter for greater than or equal to zero in z direction
    logicaldata = z>=0;
    x = x.*logicaldata;
    y =  y.*logicaldata;
    z =  z.*logicaldata;
    
    % ----- Plot Over Time Step ----- %
    for k = 1:length(t);
        if sum(z(:,k)>0) %stop plotting if all of the balls are at zero 
                         % in z-direction
            figure(2)
            plot3(x(:,k),y(:,k),z(:,k),'.k');
            title('Launched Projectiles from Collision','FontSize',20);
            xlabel('Distance in x Direction (meters)','FontSize',15);
            ylabel('Distance in y Direction (meters)','FontSize',15);
            zlabel('Height (meters)','FontSize',15);
            hold on;
            grid on;
            axis tight;
            f = f + 1;
            saveas(gcf, ['jupcoll_',sprintf('%03d',f),'.png'])
            pause(0.001);
        end;
    end;
end;
end

