%	***** DEMONSTRATION OF EXPONENTIAL WINDOW METHOD *****
%	 Continuous system: Rod subjected to triangular pulse
%                           Top of rod is free, bottom is fixed

pi2=2*pi;		% 6.28318...

% *** Rod data
C=5;			% speed of waves in rod, in [m/ms]
R=2;			% mass density, in [ton/m^3]   (1 ton=1000 kg)
E=R*C^2;		% Modulus of elasticity
A=0.0001;		% rod cross-section, in [m^2]
L=0.5;			% rod length, in [m]
P=1;			% peak force amplitude
x=0.25;			% output location (distance above support)
td=0.02;		% Duration of forcing function, in [ms]
tp=0.512;		% period of response, in [ms]

% *** FFT data
nn=256;			% number of points in time domain
nf=nn/2+1;		% number of points in frequency domain (0 to Nyquist)
dt=tp/nn;		% time step
df=1/tp;		% frequency step, in kHz
dw=pi2*df;		% frequency step, in krad/s
nyq=0.5/dt;        	% nyquist frequency, in kHz
time=(0:dt:tp-dt);	% time axis vector
freq=(0:df:nyq);	% frequency axis vector (in Hz)
om=(0:dw:pi2*nyq);	%     "       "    "    (in rad/s)

% *** Define forcing function
nr=round(0.5*td/dt);	% number of intervals during which force is rising
td=2*nr*dt;		% redefine duration of force --> even No. of intervals
p=P/nr;
force=p*[(0:nr),(nr-1:-1:0), zeros(1,nn-2*nr-1)];   % define forcing function
plot(time,force);	% plot forcing function
title ('Time history of load')
xlabel('Time (ms)')
grid on
pause

% *** Apply exponential window and FFT to forcing function
a=dw;			% imaginary component of frequency
omc=om-i*a;		% complex frequency array
alfa=omc/C;		% complex wavenumber
e=exp(-time*a);		% decaying exponential window
e1=1./e;		% rising exponential window
force=force.*e;		% apply window to forcing function
F=fft(force*dt);	% FFT of forcing function
plot(freq,abs(F(1:nf)));% plot Fourier transform of forcing function
title('Fourier transform of load')
xlabel('Frequency (kHz)')
grid on
pause

% *** Compute displacement response ***
flex=1/(E*A);		% rod flexibility factor
h=flex*sin(alfa*x)./cos(alfa*L);
h=h./alfa;		% transfer function for displacement
plot(freq,abs(h));	% plot transfer function, absolute value
title('Displacement transfer function')
xlabel('Frequency (kHz)')
grid on
pause

H=[h,conj(h(nf-1:-1:2))];	% creating mirror image for FFT 
u=real(ifft(F.*H*df));	% Inverse FFT
u=e1.*u;		% apply rising exponential window
plot(time,u);		% plot response
title('Point displacement, time history')
xlabel('Time (ms)')
grid on
pause

% *** Compute velocity response ***
h=i*omc.*h;		% velocity transfer function
plot(freq,abs(h));	% plot transfer function, absolute value
title('Velocity transfer function')
xlabel('Frequency (kHz)')
grid on
pause

H=[h,conj(h(nf-1:-1:2))];	% creating mirror image for FFT 
u=real(ifft(F.*H*df));	% Inverse FFT
u=e1.*u;		% apply rising exponential window
plot(time,u);		% plot response
title('Point velocity, time history')
xlabel('Time (ms)')
grid on