%{
Code adapted from Michael Hart's "Projectile Motion":
https://www.youtube.com/watch?v=TrUa2gP1Tc8

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 sequentially 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_collision(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_collision(x0,y0,z0)

close all;

nParticles=10;      % 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 --- %
for i=1:nParticles
    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
        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 ----- %
        x = v0*t*cosd(theta) + x0;
        y = v0*t*sind(theta) + y0;
        z = v0*t*sind(phi) + -0.5*a.*t.^2;
        
        % ----- Plot Over Time Step ----- %
        for k = 1:length(t);
            if x(k) >=0 && y(k)>= 0 && z(k) >= 0 ;
                plot3(x(k),y(k),z(k),'.k');
                title('Sequentially Launched Projectiles','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, ['volc_',sprintf('%03d',f),'.png'])
                pause(0.1);
            end;
        end;
    end;
end
end
