
clear
close all


prompt = {'w0 in Hz (1)','Damping (0.5) '};
dlg_title = 'Enter Seismometer characteristics';
num_lines=1;
def ={'1','.5'};
answer  = inputdlg(prompt,dlg_title,num_lines,def);
w0= str2num(char(answer(1)));tsau=str2num(char(answer(2)));

% Simple illustration of the poles and zeros concept of a linear system

% Case 1: the seismometer
%
% Laplace transformed transfer function:
% H(s) = s^2/(s^2+2 tsau omega0 s + omega0^2)
%
% written as ratio of two polynomials and s= Jw
%
%                           (Jw-0) Jw-0)
% H(s)= --------------------------------------------------------------
%       [tsau w0+J(w+ w0sqrt(1-tsau^2)][tsau w0+J(w- w0sqrt(1-tsau^2)]

% seismometer parameters
%w0=1.;           % Eigenfrequency in Hz
%tsau = .1;       % Damping characteristics
w=(.1:.01:10);    % Frequency range

% ok,  let s visualize it

% poles

z1 = 0;
z2 = 0;
p1= -w0*(tsau+j*sqrt(1-tsau^2));
p2= -w0*(tsau-j*sqrt(1-tsau^2));

% Plotting

for i=1:length(w);
H(i)=-((j*w(i)-z1)*(j*w(i)-z2))/((j*w(i)+p1)*(j*w(i)+p2));
end

subplot(221)
text(0,0,'0')
text(real(p1),imag(p1),'P1')
text(real(p2),imag(p2),'P2')
xlabel(' real part ')
ylabel(' imaginary part ')
title(' Poles and zeros ')
axis([ -1.25*abs(real(p1)) 1.25*abs(real(p1)) -1.25*abs(imag(p1)) 1.25*abs(imag(p1))])
grid
axis square

subplot(222),
plot(log10(w),log10(abs(H)))
xlabel(' Log w ')
ylabel(' Log Amplitude ')
title(' Amplitude response ')
axis square

subplot(223),
plot(log10(w),atan(real(H)./imag(H)))
xlabel(' Log w ')
ylabel(' Phase ')
title(' Phase response ')
axis square

subplot(224),
plot(log10(w),log10(1./abs(H)))
xlabel(' Log w ')
ylabel(' Log Amplitude ')
title(' Instrument Correction (Ampl.) ')
axis square


