%Unit Step Response for a Cantilever Beam

clc
clear all
close all

%Beam Properties
p = 2700; %kg/m^3
E = 69e9; %Pa
L = 0.1; %m
B = 0.01; %m 
H = 0.01; %m 
V = L*B*H; %m^3
m = p*V; %kg
I = (1/2)*B*H^3; %m^4
k = (3*E*I)/(L^3); %N/m
C = 10; %kg/s (UNKNOWN)
F = 10000; %N
t = 0:0.000025:0.05; %s
wn = sqrt(k/m); %rad/s
z = C/(2*sqrt(k*m)); %unitless
wd = wn*sqrt(1-z^2); %rad/s
T = (2*pi)/wd; %s
Fs = 1/0.000025; %1/s  sampling frequency
Ts = 1/Fs; %s sample time


x = (F/k)*((1-exp(-z*wn*t).*(cos(wd*t)+(wn/wd)*sin(wd*t)))-(1-exp(-z*wn*(t-T)).*(cos(wd*(t-T))+z*(wn/wd)*sin(wd*(t-T)))).*heaviside(t-T)); %unit step response

figure, plot(t,x,'r'),xlabel('t [s]'),ylabel('x(t) [m]'),title('Unit Step Response'); 

Nsamps = length(x);
t = (1/Fs)*(1:Nsamps);          %Prepare time data for plot

%Do Fourier Transform
y_fft = abs(fft(x));            %Retain Magnitude
y_fft = y_fft(1:Nsamps);      
f = Fs*(0:Nsamps-1)/Nsamps;   %Prepare freq data for plot

%Plot in Frequency Domain
figure
plot(f, y_fft)
xlim([0 1000])
xlabel('Frequency (Hz)')
ylabel('Amplitude')
title('Frequency Response')