function f = spring_mass_simulation(m,k,b)

% spring_mass_simulation(m,k,b)
%
% This function solves the equations of motion of a spring-mass system
% and the animation showing the solution is created.
% The input variables are: mass, m, [kg], stiffness, c, [N/m], and
% friction coefficient, b, [N*s/m].
%
% The simulation is nominally runs for T=2 seconds. Feel free to change
% this parameter if necessary.
%
% This code is for educational purposes only. Do not use or distribute
% outside of BME 130
%
%Prof. Zoran Nenadic
%UC Irvine
%08/02/2005

close all

disp(['natural frequency: ' num2str(sqrt(k/m))])
disp(['damping factor: ' num2str(sqrt(b/(2*sqrt(k*m))))])

g = 9.81;   % define gravity acceleration [m/s^2]
T = 10;      % define the time of simulation [s]
Npt = 30;  % define the number of frames per second of animation

t = linspace(0,T,Npt*T);  % define the time vector

% COORDINATES OF THE SPRING
x_spr = [0.5 0.5; 0.5 1.0; 1 0; 0 1; 1 0; 0 1; 1 0.5; 0.5 0.5];
y_spr = [0 -1; -1 -1.5; -1.5 -2.5; -2.5  -3.5; ...
    -3.5 -4.5; -4.5 -5.5; -5.5 -6.0; -6.0 -7.0]*0.05;

ylim = [-2.75*m*g/k+y_spr(end,2) 0];  % define limit for y-axis

% SCALES OF THE MOTION OF DIFFERENT SPRING PARTS WHICH MAKE THE SPRING
% APPEAR ELASTIC
c = [0 0; 0 1/6; 1/6 2/6; 2/6 3/6; 3/6 4/6; 4/6 5/6; 5/6 6/6; ...
    1 1];

% PLOT THE INITIAL POSITION OF THE SPRING AND MASS
hold on
for i = 1:length(x_spr)
    sh(i) = plot(x_spr(i,:),y_spr(i,:),'b-','LineWidth',5);
end
set(gcf,'Position',[500 100 800 600])

% PLOT THE INITIAL POSITION OF THE MASS
mh = plot(x_spr(end,2),y_spr(end,2),'ko','LineWidth',5);
msize = round(0.1*2*m*g/(k*diff(ylim))*420); % adaptive marker size
set(mh,'MarkerSize',msize,'MarkerFaceColor','r');

% PLOT THE UPPER AND LOWER BOUNDS OF THE MASS MOTION
plot([-0.5 1.5],y_spr(end,2)*[1 1],'k--','LineWidth',2)
plot([-0.5 1.5],y_spr(end,2)*[1 1]-2*m*g/k,'k--','LineWidth',2)


xh = xlabel(['Time: ']);
yh = ylabel('Position [m]');
th = title(['m=' num2str(m) '  [kg]        k=' num2str(k) ...
    '  [N/m]        b=' num2str(b) '  [kg/s]']);

if b^2 < 4*k*m
    msg ='periodic';
else
    msg = 'aperiodic';
end 

set([xh yh th],'FontName','Times','FontSize',16, ...
    'FontWeight','b','FontAngle','i');

set(gca,'XLim',[-1.5 2.5],'YLim',ylim, ...
    'XTick',[],'FontName','Times','FontWeight','b', ...
    'FontSize',16);    % set the axis properties
set(gcf,'Color',[0.99 0.87 0.0])

uh = uicontrol('Style','text','Position',[390 10 150 25], ...
    'String',upper(num2str(msg)),'FontSize',16,'FontName','Times', ...
    'FontWeight','b','FontAngle','i','ForeGroundColor','w', ...
    'BackGroundColor','k');

disp('hit any key to continue')
pause

disp('solving ...');    % solving should go fast

% MATLAB ODE SOLVER
[TOUT,YOUT] = ode45(@smf,t,[0 0]);  % the function is spring-mass fcn
    function dx = smf(t,x)
        dx = zeros(2,1);    % a column vector
        dx(1) = x(2);
        dx(2) = -k/m*x(1) - b/m*x(2) + g;
    end

x1 = YOUT(:,1);     % first column position (second column velocity)

disp('plotting ...')
str = '.';
for i = 1:length(t)
    
    %update time
    set(xh,'String',['Time: ' num2str(t(i),'%2.2f') ' [s]']);
    for j = 1:length(x_spr)

        %update spring segment (vertical) positions
        set(sh(j),'YData',y_spr(j,:)-x1(i)*c(j,:));
    end

    %update the vertical position of the mass
    set(mh,'YData',y_spr(end,2)-x1(i)*c(j,end))
    drawnow
    
    %comment this out for super fast animation
    pause(0.01)
end

f = 'EXIT_SUCCESS';
end