% PIX_2.m
% particle drop on a spring 

clear; clc; close all

g = 9.81; % gravitational acceleration (m/s^2)
m = 10; % mass (kg)
k1 = 294e3; % spring constant (N/m)
h = 1; % drop height (m)
v0 = sqrt(2*g*h); % initial velocity (m/s)

% equation of motion
% m x'' = m g-k1*x

syms t x1 x2

dx1 = x2;
dx2 = g-k1*x1/m;

dx = [dx1; dx2];

% g=matlabFunction(dx,'file','my_dx','vars',{t, [x1; x2]});
% matlabFunction creates a MATLAB file or 
% anonymous function from a sym object

g = matlabFunction(dx,'vars',{t, [x1; x2]});

% options for ODE solver
optionx1 = odeset('RelTol',1e-6,'MaxStep',1e-6,...
     'Events',@eventx1);
% odeset creates ODE options structure
% RelTol - relative error tolerance 
% MaxStep - upper bound on step size 
% Events - detects events  
 
ic0 = [0 v0]; % initial conditions

t0 = 0; tf = 1; ti = [t0 tf]; % time

[tte, xse, te, ye] = ode45(g, ti, ic0, optionx1);

tf = real(te(2));
vf = real(ye(4));

fprintf(' \n')
fprintf('tf = %6.3e (s)\n',tf)
fprintf('vf = %6.3f (m/s)\n',vf)

time = real(tte);
dis = real(xse(:,1));
vel = real(xse(:,2));

Pen = k1*dis;
xmax = max(dis);
Pmax = max(Pen);
fprintf('x max = %6.3e (m)\n',xmax)
fprintf('P max = %6.3e (N)\n',Pmax)

% plots
figure(1)
grid on
subplot(3,1,1)
plot(time,dis,'LineWidth',1.5)
ylabel('x (m)')
grid on
subplot(3,1,2)
xlabel('t (s)')
plot(time,vel,'LineWidth',1.5)
ylabel('dx/dt (m/s)')
grid on
subplot(3,1,3)
plot(time,Pen,'LineWidth',1.5)
ylabel('Pe (N)'), xlabel('t (s)')
grid on

print -depsc f2.eps 



function [value,isterminal,direction] = eventx1(~,x)

% y  = x(1);
% y' = x(2);

value = x(1);

isterminal = 1;
direction = 0;
end

% if isterminal vector is set to 1, integration will halt 
% when a zero-crossing is detected.

% The elements of the direction vector are -1, 1, or 0, 
% specifying that the corresponding event must be decreasing, 
% increasing, or that any crossing is to be detected






% tf = 1.847e-02 (s)
% vf = -4.429 (m/s)
% x max = 2.617e-02 (m)
% P max = 7.694e+03 (N)


% total time = 1.847e-02 (s)
% max displacement  = 2.617e-02 (m)
% max velocity =  4.430 (m/s)
% max elastic force = 7693.653 (N)
