Initial commit.

Author: zcshiner

ballistic.m

@@ -0,0 +1,37 @@
+function [T,motionTotal] = ballistic(p0, theta, phi, velocity) % angles are in radians
+% BALLISTIC(ICS, azimuth, elevation, v0)
+% Calculate the path of a ballistic object with air resistance.
+% aizimuth is in degrees along xy plane from the x-axis, elevation is in
+% degrees from the xy plane to the initial velocity vector.
+%
+% Equations via: http://www.cs.cornell.edu/~bindel/class/cs3220-s12/notes/lec27.pdf
+
+g = [0,0,9.81]'; % m/s % xyz
+c = 0.45; %1/m; % scaled drag coefficient
+
+% decompose heading & pitch into xyz magnitudes times velocity scalar
+v0 = velocity.*[cos(theta)*cos(phi), sin(theta)*cos(phi), sin(phi)]';
+ICS = [p0';v0];
+tFinal = v0(3)/g(3)+sqrt((v0(3)/g(3)).^2+2*p0(3)/g(3));
+% for testing: v = v0
+
+    % Ballistic model, considering air drag
+    function rhs = dragModel(t,motion)
+        v = motion(4:6); % last 3 terms are velocity
+        rhs = [v; (-g-c.*abs(v).*v)]; %[vx,vy,vz,ax,ay,az]'
+    end
+
+    % ODE options: stop when we hit zero
+    function [value,isterminal, direction ] = hitground(t,y)
+        value = y(3); % Check for zero crossings of z position
+        isterminal = 1; % Terminate on a zero crossing
+        direction = -1; % We only care about crossing y > 0 to y < 0
+    end
+
+odeopt = odeset('Events', @hitground);
+[T,XV] = ode45(@dragModel,[0,tFinal.^2],ICS,odeopt);
+
+% add in acceleraton data.. because I can.
+motionTotal = horzcat(XV,[-c.*abs(XV(:,4:5)).*XV(:,4:5),-g(3)-c.*abs(XV(:,6)).*XV(:,6)]);
+
+end