%Math 22B, SS1 2007, damped mass-spring system solver
%Christopher L. Wheaton, 4250

%This MATLAB script plots the solution to a damped mass-spring system.
%This script also draws two phase portraits and makes a trajectory plot.
%This system is assumed to initially be at rest and at maximum amplitude.
%This system is modeled by the diff. eq. mx'' + bx' + kx = 0.
%There is no external forcing function for this system (g(t) = 0).
%In the overdamped case, the discriminant (b^2-4km) is greater than 0.
%In the critically damped case, the discriminant is equal to 0.
%In the underdamped case, the discriminant is less than 0.
%In this script, parameters such as, time duration, number of time
%intervals, spring constant, mass of object, and initial displacement,
%can be modified by the user to fit a specific application.

%In general, as the number of time intervals increases, the smoother the
%plots. Although, since this script is based on Euler's method, the error
%compounds as subsequent points are computed. This is because the method 
%employed requires the previous point to calculate the next.

clear %clear all variables
clc %clear command window

%input parameters
t = 20 %time duration (in seconds)
N = 1000 %number of time intervals
b = 2 %damping coefficient (in Newtons/(meters/second))
%assuming k and m = 1
%choose b > 2 for overdamped
%choose b = 2 for critical damping
%choose 0 < b < 2 for underdamped
%choose b = 0 for no damping
%choose b < 0 for negative damping
k = 1 %spring constant (in Newtons/meter)
m = 1 %mass of object (in kilograms)
L = 2 %initial displacement (in meters)

%definitions
time = linspace(0,t,(N+1)); %array of times
dt = t/N %time interval size (in seconds)
x = zeros(1,N+1); %fill x array with zeros
x(2) = L; %initial displacement
x(1) = L; %since slope is zero (specified condition), use same value 
x_dot = zeros(1,N+1); %fill x_dot array with zeros
x_double_dot = zeros(1,N+1); %fill x_double_dot array with zeros
x_double_dot(1) = -k*L/m; %known accel. because v = 0 at max amplitude

%perform computations to fill x (position) array
for n = 2:N %run loop from n = 2 to n = N
    x_dot(n) = (x(n)-x(n-1))/dt; %time derivative with respect to x
    x_double_dot(n) = (-b*x_dot(n)-k*x(n))/m; %second time derivative
    x(n+1) = x(n)+x_dot(n)*dt+(dt^2/2)*x_double_dot(n); %insert next value
end

%plot position of mass-spring system versus time
figure(1) %bring figure window 1 to the foreground
clf %clear figure window
hold on %perform subsequent plots in the same figure window
plot(time,x) %plot the array of values
plot(time,0) %plot constant function zero
title('Position of Mass-Spring System vs. Time')
xlabel('time (seconds)')
ylabel('position (meters)')
grid on

%plot phase portraits
figure(2) %bring figure window 2 to the foreground
clf %clear figure window
plot(x,x_dot) %plot x_dot vs. x
title('Velocity vs. Position')
xlabel('position (meters)')
ylabel('velocity (m/s)')
grid on

figure(3) %bring figure window 3 to the foreground
clf %clear figure window
plot(x_dot,x_double_dot) %plot x_double_dot vs. x_dot
title('Acceleration vs. Velocity')
xlabel('velocity (m/s)')
ylabel('acceleration (m/s^2)')
grid on

%plot trajectory
figure(4) %bring figure window 4 to the foreground
clf %clear figure window
plot3(x,x_dot,x_double_dot)
title('Trajectory Plot')
xlabel('position (meters)')
ylabel('velocity (m/s)')
zlabel('acceleration (m/s^2)')
grid on

%end of script