function HH4
% EN40 Sprintg 2016 Homework 4


close all

%%
data_Reg = csvread('Lax_tests/t2.csv',2,0);   % Read data for Lax ball
Regballmass = 0.149;  % Regulation ball mass, kg

% Plot the graphs
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
% Create plot
plot(data_Reg(:,2),data_Reg(:,3),'LineWidth',2,'Color',[1 0 0]);
xlabel('Increase in diameter (mm)');
ylabel('Force (kN)');
title('Static force-deflection curve for regulation Lacrosse ball');
box(axes1,'on');
set(axes1,'FontSize',14,'XGrid','on','YGrid','on');
annotation(figure1,'textarrow',[0.553061224489796 0.448979591836735],...
    [0.490833030852995 0.584392014519056],'String',{'Loading'},'FontSize',14);
annotation(figure1,'textarrow',[0.342857142857143 0.459183673469388],...
    [0.745916515426497 0.647912885662432],'String',{'Unloading'},'FontSize',14);

% Calculate the integral of force wrt distance during compression
[value,index] = min(data_Reg(:,3));  % Find the max force

% Do the integral - the change in diameter is halved to convert it into the
% relevant deflection
Work_loading = trapz(data_Reg(1:index,2)/2,data_Reg(1:index,3));
impact_velocity_regular = sqrt(2*Work_loading/Regballmass);
fprintf('Regular ball impact velocity %g \n', impact_velocity_regular)

% Calculate the work done during rebound
Work_unloading = -trapz(data_Reg(index:end,2)/2,data_Reg(index:end,3)); 
rebound_velocity_regular = sqrt(2*Work_unloading/Regballmass);
fprintf('Regular ball rebound velocity %g (m/s)\n', rebound_velocity_regular);

%Calculate the restitution coefficient
e_regular = rebound_velocity_regular/impact_velocity_regular;
fprintf('Regular ball restitution coefficient %g (m/s)\n\n', e_regular);

data_Swax = csvread('Lax_tests/ts3.csv',2,0);  % Read data for Lax ball
Swaxballmass = 0.145; % Swax ball mass, kg

figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
% Create plot
plot(data_Swax(:,2),data_Swax(:,3),'LineWidth',2,'Color',[1 0 0]);
xlabel('Increase in diameter (mm)');
ylabel('Force (kN)');
title('Static force-deflection curve for SwaxLax Lacrosse ball');
box(axes1,'on');
set(axes1,'FontSize',14,'XGrid','on','YGrid','on');
annotation(figure1,'textarrow',[0.361512791991101 0.248053392658509],...
    [0.734782608695652 0.895652173913044],'String',{'Unloading'},'FontSize',14);
% Create textarrow
annotation(figure1,'textarrow',[0.345052325713378 0.240970693060317],...
    [0.59735476998343 0.69091375364949],'String',{'Loading'},'FontSize',14);


[value,index] = min(data_Swax(:,3));  % Find the max force

% Do the integral - the change in diameter is halved to convert it into the
% relevant deflection
Work_loading = trapz(data_Swax(1:index,2)/2,data_Swax(1:index,3));
impact_velocity_Swax = sqrt(2*Work_loading/Swaxballmass);
fprintf('Swax ball impact velocity %g (m/s)\n', impact_velocity_Swax)

% Calculate the work done during rebound
Work_unloading = -trapz(data_Swax(index:end,2)/2,data_Swax(index:end,3)); 
rebound_velocity_Swax = sqrt(2*Work_unloading/Swaxballmass);
fprintf('Swax ball rebound velocity %g (m/s)\n', rebound_velocity_Swax)
e_swax = rebound_velocity_Swax/impact_velocity_Swax;
fprintf('Swax ball restitution coefficient %g \n', e_swax);



a = 9.81;
wn = 41000*2*pi;
zeta = 0.05;
wd = wn*sqrt(1-zeta^2);
for i=1:201
   t = (i-1)*100/200/wn;
   x(i) = -a/wn^2 + a/wn^2*exp(-zeta*wn*t)*(cos(wd*t)+(zeta*wn/wd)*sin(wd*t));
   lowerbound(i) = -0.95*a/wn^2;
   upperbound(i) = -1.05*a/wn^2;
   tplot(i) = t;
end

figure
plot(tplot,x);
hold on
plot(tplot,lowerbound);
plot(tplot,upperbound);



end