function hw2
%% ENGN0040 Dynamics and Vibrations 2018
%  Homework #2
%  Problem 2

close all

m = 4.5e-03;    % Pill mass (kg)

%% 2.1
% Read the data
data_Fed = csvread('Fed_Rat2.csv',1,0);
data_Fasted = csvread('Fasted_Rat2.csv',1,0);


%% 2.2
% Plot the trajectories
%
plot3(data_Fed(:,2),data_Fed(:,3),data_Fed(:,4))
 grid('on');
 ylabel('y (mm)');
 xlabel('x (mm)');
 zlabel('z (mm)');
 title('Trajectory of a pill inside a fed rat');
figure
plot3(data_Fasted(:,2),data_Fasted(:,3),data_Fasted(:,4))
 grid('on');
 ylabel('y (mm)');
 xlabel('x (mm)');
 zlabel('z (mm)');
 title('Trajectory of a pill inside a fasted rat');

%% 2.3
% Calculate the velocities, accelerations, and deduce the magnitude of the
% force as a function of time
%
%  Here's the force calculation using the diff function and MATLAB dot notation
% The accelerations are in mm/s^2 and mass in grams so we get forces in
% microNewtons
time = data_Fed(:,1);
v_Fed = diff(data_Fed(:,2:4))/0.1;
a_Fed = diff(v_Fed(:,1:3))/0.1;
Fmag_Fed = m*sqrt(a_Fed(:,1).^2 + a_Fed(:,2).^2 + a_Fed(:,3).^2);
figure
plot(time(1:length(Fmag_Fed(:))),Fmag_Fed(:),'Color',[1,0,0],'LineWidth',1);
 grid('on');
 ylabel('Force (m N)');
 xlabel('Time (sec)');
 title('Forces on a pill inside a fed rat');

time = data_Fasted(:,1);
v_Fasted = diff(data_Fasted(:,2:4))/0.1;
a_Fasted = diff(v_Fasted(:,1:3))/0.1;
Fmag_Fasted = m*sqrt(a_Fasted(:,1).^2 + a_Fasted(:,2).^2 + a_Fasted(:,3).^2);

figure
plot(time(1:length(Fmag_Fasted(:))),Fmag_Fasted(:),'Color',[1,0,0],'LineWidth',1);
 grid('on');
 ylabel('Force (m N)');
 xlabel('Time (sec)');
 title('Forces on a pill inside a fasted rat');
 
 
% Here's the calculation using loops (working component by component is
% fine too)
[rows,cols] = size(data_Fed);
v_Fed = zeros(rows-1,3);
for i=1:rows-1
   v_Fed(i,1:3) = (data_Fed(i+1,2:4) - data_Fed(i,2:4))/(data_Fed(i+1,1)-data_Fed(i,1)); 
end
Fmag_Fed = zeros(rows-2,1);
for i=1:rows-2
   a_Fed(1:3) = (v_Fed(i+1,1:3) - v_Fed(i,1:3))/(data_Fed(i+1,1)-data_Fed(i,1));
   Fmag_Fed(i) = m*sqrt(a_Fed(1)^2 + a_Fed(2)^2 + a_Fed(3)^2);
end

[rows,cols] = size(data_Fasted);
v_Fasted = zeros(rows-1,3);
for i=1:rows-1
   v_Fasted(i,1:3) = (data_Fasted(i+1,2:4) - data_Fasted(i,2:4))/(data_Fasted(i+1,1)-data_Fasted(i,1)); 
end
Fmag_Fasted = zeros(rows-2,1);
for i=1:rows-2
   a_Fasted(1:3) = (v_Fasted(i+1,1:3) - v_Fasted(i,1:3))/(data_Fasted(i+1,1)-data_Fasted(i,1));
   Fmag_Fasted(i) = m*sqrt(a_Fasted(1)^2 + a_Fasted(2)^2 + a_Fasted(3)^2);
end

%% Find the average forces on the pills by integrating the forces....

avF_fed = trapz(time(1:length(Fmag_Fed)),Fmag_Fed)/(time(length(Fmag_Fed))-time(1));
avF_fasted = trapz(time(1:length(Fmag_Fasted)),Fmag_Fasted)/(time(length(Fmag_Fasted))-time(1));

str = compose('Average force on the pill in a fed rat %5g',avF_fed);
disp(str)
str = compose('Average force on the pill in a fasted rat %5g',avF_fasted);
disp(str)


end