% damped_spring_response.m

%

% This MATLAB script file demonstrates the

% form of the response function for a

% damped spring.

%

% K. Beers. MIT ChE. 12/14/2002

% For discrete Fourier transform, use the

% MATLAB convention

epsilon = -1;

% Now, set the spring constants.

zeta = 1;  % drag constant

K_s = 1;  % spring constant

mass= 1;  % mass

% Generate an appropriate time step and sampling time vector

omega_s = sqrt(K_s/mass);  % natural frequency of spring

P_s = pi/omega_s;  % half-period of spring

if(zeta > sqrt(eps))

    lambda = mass/zeta;

else

    lambda = mass/sqrt(eps);

end

P_sample = max(P_s, min(lambda,10*P_s));

dt = P_s / 100;  % set sufficient sampling period

t0 = 0;  % initial time

t_end = 6*P_s;

t_val = [t0 : dt : t_end];

N = length(t_val);  % # of sampling points

if(~rem(N,2))  % if N is odd

    N = N + 1;

    t_val = [t_val, t_end + dt];

end

P_data = N*dt;

% set frequency values for this data set

d_omega = pi/P_data;

omega_max = (N-1)*d_omega;

omega = [0 : d_omega : omega_max];

% calculate Fourier transform of response function

r_FT_real = zeros(size(omega));

r_FT_imag = zeros(size(omega));

for m=1:length(omega)

    omega_m = omega(m);

    harmonic = omega_s^2 - omega_m^2;

    viscous = (zeta/mass)*omega_m;

    denom = harmonic^2 + viscous^2;

    r_FT_real(m) = harmonic/mass/denom;

    r_FT_imag(m) = epsilon*viscous/mass/denom;

end

r_FT = complex(r_FT_real,r_FT_imag);  % construct complex F.T.

% Now, convert to discrete Fourier transform format.

r_FT = (sqrt(2*pi)/dt)*r_FT;

% calculate response function from inverse Fourier transform

r = real(ifft(r_FT));

% plot response function in Fourier and time space

figure;

% power spectrum of response function.  Only plot up to

% a certain magnitude of the characteristic spring

% frequency.

subplot(2,1,1);

r_PS = real(conj(r_FT).*r_FT);

plot(omega,r_PS);

omega_plot_max = 2*omega_s;

axis([0 omega_plot_max min(r_PS) max(r_PS)]);

xlabel('\omega');

ylabel('|R(\omega)|^2');

phrase_val = ['\zeta = ', num2str(zeta), ...

        ', K_s = ', num2str(K_s), ...

        ', m = ', num2str(mass)];

title(['Response of damped spring, ', phrase_val]);

% time plot of response function

subplot(2,1,2);

plot(t_val,r);

xlabel('t');

ylabel('r(t)');