% Matlab example: Second order dynamic response of an impulse function

clc
clear all
close all

% Two parameters needed to define the dynamic response of a second order system
% Natural frequency wn (omega n)
% Damping ratio d (xi)
% It is easier for the user to define the natural frequency by defining the period of 
% oscillations in seconds.

% input function broken in Freemat, define T_d for now
T_d = input('Enter the period of observed oscillations (s): ');

% from T_d, find damped natural frequency f_d (Hz) and damped natural frequency w_d (rad/s)
f_d = 1/T_d;
w_d = 2*pi*f_d;

% next, take user input for underdamped and overdamped natural frequencies. 
d_under  = input('Enter a damping ratio between 0 and 1 for an underdamped system: ');
d_over = input('Enter a damping ratio greater than 1 for the overdamped system: ');

% find undamped natural frequency and phase shift phi for the underdamped system.
f_n = f_d * sqrt(1-d_under^2);
w_n = 2*pi*f_d;
phi = asin(sqrt(1-d_under^2));

% for the impulse function, take user input for the value of the impulse spike. 
y_i = input('Enter the initial value of y: ');
y_f = input('Enter the impulse value of y: ');

% define time. 10 periods seems like a good amount of time. 
t = [0:.001:10*T_d];
% natural frequency * time is a good time step for normalized y.
w_nt = w_n.*t;

% Want to display zero damping case, underdamped, critically damped, and overdamped. 
% By definition, the impulse response of a second-order system is simply a differentiation of the step response.
% Shown below are solutions for zero, underdamped, critical damped, and overdamped systems.

% zero damping, d = 0
d_zero = 0;
phi_zero = asin(sqrt(1-d_zero^2));
y_norm_zero_impulse = w_n*exp(-d_zero*w_n.*t).*((d_zero/sqrt(1-d_zero^2))*sin(w_n.*t*sqrt(1-d_zero^2)+phi_zero)-cos(w_n.*t*sqrt(1-d_zero^2)+phi_zero));
y_zero = y_norm_zero_impulse *(y_f-y_i) + y_i;

% user defined underdamped 
y_norm_under_impulse = w_n*exp(-d_under*w_n.*t).*((d_under/sqrt(1-d_under^2))*sin(w_n.*t*sqrt(1-d_under^2)+phi)-cos(w_n.*t*sqrt(1-d_under^2)+phi));
y_under = y_norm_under_impulse * (y_f-y_i) + y_i;

% critical damping (d = 1)
y_norm_crit_impulse = -w_n.*exp(-w_n.*t).*(1-(1+w_n.*t));
y_crit = y_norm_crit_impulse *(y_f-y_i) + y_i;

% user defined overdamped
y_norm_over_impulse = w_n.*exp(-d_over.*w_n.*t).*((d_over.*((cosh(w_n.*t*sqrt(d_over^2 - 1)) + (d_over/sqrt(d_over^2 - 1)).*sinh(w_n.*t*sqrt(d_over^2 - 1)))))-((sqrt(d_over^2 - 1).*sinh(w_n.*t*sqrt(d_over^2-1)))+(d_over.*cosh(w_n.*t*sqrt(d_over^2-1)))));
y_over = y_norm_over_impulse *(y_f-y_i) + y_i;

figure(1)
hold on
% plot(w_nt, y_norm_zero_impulse)
plot(w_nt, y_norm_under_impulse, 'b')
plot(w_nt, y_norm_crit_impulse, 'g')
plot(w_nt, y_norm_over_impulse, 'r')
% title, axis labels, and legend
title('Second order normalized impulse response')
xlabel('nondimensional time, w_n*t')
ylabel('normalized y')
legend('Underdamped','Critically Damped','Overdamped','Location','SouthOutside','Orientation','horizontal')


figure(2)
hold on
%plot(t, y_zero)
plot(t, y_under, 'b')
plot(t, y_crit, 'g')
plot(t, y_over, 'r')
% title, axis labels, and legend
title('Second order impulse response')
xlabel('time (s)')
ylabel('y')
legend('Underdamped','Critically Damped','Overdamped','Location','SouthOutside','Orientation','horizontal')


% What does the response look like if there is another impulse before this one
% settles?

% take user input for the magnitude of a second impulse and the time at which
% it occurs.

t_2 = input('Enter the time in seconds the second step occurs: ');
y_f2 = input('Enter the magnitude of the second step: ');

% define interval. The response to the second impulse needs to begin at t2 = 0
n = (t_2/.001);
t2 = [0:.001:(10*T_d - t_2)];
y_norm_under_impulse2 = w_n*exp(-d_under*w_n.*t2).*((d_under/sqrt(1-d_under^2))*sin(w_n.*t2*sqrt(1-d_under^2)+phi)-cos(w_n.*t2*sqrt(1-d_under^2)+phi));
y_under_impulse2 = [zeros(1,n), y_norm_under_impulse2 * (y_f2-y_i)];
% the actual function is zero until the t_2 when the impulse occurs.
y_norm_crit_impulse2 = -w_n.*exp(-w_n.*t2).*(1-(1+w_n.*t2));
y_crit_impulse2 = [zeros(1,n), y_norm_crit_impulse2 *(y_f2-y_i)];
y_norm_over_impulse2 = w_n.*exp(-d_over.*w_n.*t2).*((d_over.*((cosh(w_n.*t2*sqrt(d_over^2 - 1)) + (d_over/sqrt(d_over^2 - 1)).*sinh(w_n.*t2*sqrt(d_over^2 - 1)))))-((sqrt(d_over^2 - 1).*sinh(w_n.*t2*sqrt(d_over^2-1)))+(d_over.*cosh(w_n.*t2*sqrt(d_over^2-1)))));
y_over_impulse2 = [zeros(1,n), y_norm_over_impulse2 *(y_f2-y_i)];

% to superimpose, simply add the two functions for each damping constant. 
y_under_super = y_under + y_under_impulse2;
y_crit_super = y_crit + y_crit_impulse2;
y_over_super = y_over + y_over_impulse2;

% plot the superimposed responses. 
figure(3)
hold on
plot(t, y_under_super, 'b')
plot(t, y_crit_super, 'g')
plot(t, y_over_super, 'r')
% title and axis labels
title('Second order impulse response')
xlabel('time (s)')
ylabel('y')
% add a plot of the impulses and legend
plot([0,0,0,t_2,t_2,t_2,10*T_d],[y_i,y_f,y_i,y_i,y_f2,y_i,y_i], 'k')
legend('Underdamped','Critically Damped','Overdamped','Impulse Function','Location','SouthOutside','Orientation','horizontal')
