% Problem set 1.8 
% A sphere of mass m falls from height h on a vertical spring makes contact
% with the spring as shown in Fig. PS1.8. The spring with the stiffness k
% compresses under the weight of the sphere. Find: 1. Total time of contact
% of the sphere with the spring; 2. the relative displacement of the
% spring; 3. the velocity jump of the sphere (before the contact and after
% the contact with the spring); 4. the maximum elastic force. Verify with
% the results obtained from a computer program. Numerical application: m=10
% kg, h=1 m, k=294*10^3 N/m.

clear all
clc
close all

format long

%input data
m=10; h=1; k=294*10^3; g=9.81;
om=sqrt(k/m);

%initial velocity
v0=sqrt(2*g*h);
%random time set
ti=0.02;

%options for ODE solver to find the zero values crossing
options = odeset('Events',@events);

%ODE solver for function equation
solx=ode45(@equation,[0, ti],[0 v0],options)

%name all the outputs
time=solx.x';
X=solx.y';
t2=solx.xe;
xmin=solx.ye;
t1=t2/2;

%plots
figure(1)
plot(time,X(:,1)); title('The variation of displacement'); % the variation of the displacement
xlabel('Time(s)'); ylabel('Displacement(m)');

figure(2)
plot(time,X(:,2)); title('The variation of velocity'); % the variation of velocity
xlabel('Time(s)'); ylabel('Velocity(m/s)'); 

figure(3)
plot(time,k*X(:,1)); title('The elastic force'); % the variation of the elastic force
xlabel('Time(s)'); ylabel('The elastic force (N)');


%displacements
xt1=deval(solx,t1,1); %the displacement at t1 (maximum displacement);
xt2=xmin(1); %the displacement at t2 (back to initial position, the spring is not compressed);
deltalam=xt1-xt2;


%velocities
vt2=deval(solx,t2,2); %the velocity at t2 (back to initial position, the spring is not compressed);
deltavel=v0-vt2;

%Forces
Pmax=k*deval(solx,t1,1);

%Prints
fprintf('The total time of contact of the sphere with the spring is %f8 \n',t2);
fprintf('\n The relative displacement of the spring is %f8 \n',deltalam);
fprintf('\n The velocity jump is %f8 \n',deltavel);
fprintf('\n The maximum elastic force is %f8 \n',Pmax)