% script to plot the typical transfer function for CMG40T

clear all
close all

% there are three poles and three zeros
% in Hz so multiply by 2*pi to get radians/s
pi2=2*pi;

P1=(-23.56e-3 + i*23.56e-3)*pi2;
P2=(-23.56e-3 - i*23.56e-3)*pi2;
P3=-50*pi2;

Z1=0;
Z2=0;
Z3=159*pi2;

% the normalizing constant assures unit amplification over
% the flat part of the response
% note that if np and nz were unequal, the normalizing constant
% would need to be adjusted for the disproporationate number
% factors of 2pi in the conversion to radians.  That is we would
% need to multiple A by (2pi)^(np-nz).
A=-0.314;

%sensor gain is in V/m/s
%G=800;
% use a gain of 1 to verify that Hs without it really is flat
% with unit amplitude.
G=1;

% there is an additional constant in the digitizer that is converts
% volts to computer counts.  Nominally this can be determined by
% dividing the max counts by the max volts though some manufacturers
% provide calibrated measurements of this constant.  For a Guralp
% dm24 its about 7.8E+05 count/V.

% use N points with a sample rate of 100 Hz
N=100000;
Fs=100;

% get a frequency vector
f=linspace(0,Fs-Fs/N,N);

% get a Laplace vector with zero real (ie the fourier equivalent)
sigma=0;
s=sigma+i*pi2*f;

% get the transfer function
num=(s-Z1).*(s-Z2).*(s-Z3);
dem=(s-P1).*(s-P2).*(s-P3);
Hs=num./dem;
% scale by gain and normalizing factor
Hs=Hs*G*A;


% plot amplitude
subplot(2,1,1)
semilogx(f,20*log10(abs(Hs)))
grid on
hold on
% the 3-db point
semilogx(0.0328,-2.6085,'o')
title('transfer function for typical CMG40T')
ylabel('amplitude (db)')

% plot phase response
subplot(2,1,2)
semilogx(f,180*angle(Hs)/pi)
grid on
hold on
% the 3-db point
semilogx(0.0328,92.591,'o')
ylabel('phase (degrees)')
xlabel('frequency (Hz)')

orient('tall')
print -deps lnxtxfer.ps

