%----- example matlab m-file for Chapter 2 -----

load simBS01  % You may choose another simulation case.
              % ex. simBS02,...,simNB

M=4100/980.;  % total weight of the model.  see Fig. 2.2

dir='NS';   % Choose 'NS' or 'EW'
if dir=='NS',  chZ=1; chX=2;,  else,  chZ=7; chX=8;,  end

n1=8192;    % n1 : segment length for FFT
            % You may change this value. 
tmp=fft(eq(1:n1,[chZ,chX]),n1);
spec=tmp(1:n1/2+1,:);
freq=(fs/2)*[0:n1/2]/(n1/2);
omg=2*pi*freq;

fmax=1000;  % fmax : the index for the highest frequency to calculate values  
for L=1:fmax,
  kCn(L)=M*(-omg(L)^2)*(spec(L,2)/(spec(L,1)-spec(L,2))); % see Eq.(2.21),(2.32)
end
kn=real(kCn);

%----- least squares estimation using a parabolic function for Kn
fm1=0.2; fm2=0.4;  % fm1 : the lower frequency for the calculation
                   % fm2 : the upper frequency for the calculation
                   % These values are for model BS.
                   % Change this value, if you choose model NB.
[dum,ifm1]=min(abs(freq-fm1));
[dum,ifm2]=min(abs(freq-fm2));
Zk=kn(ifm1:ifm2)';
Hf=[(freq(ifm1:ifm2).^2)',freq(ifm1:ifm2)',ones(ifm2-ifm1+1,1)];
                                             %  y=ax^2+bx+c      
Pk=[Hf'*Hf]\Hf'*Zk;
kne=Pk(3);    % identified Kn

%----- plotting -- stiffness
figure(1)
subplot(3,1,1)
plot(freq(1:fmax),kn,'r-')   % The lower figure of Fig.2.3
xmax=4.0;
ymax=40;   % Max value for the Y-axis. This value is for model BS.
           % Change this value, if you choose model NB.
axis([0 xmax 0 ymax])
xlabel(['frequency (Hz)'])
ylabel(['stiffness (tonf/cm)'])
text(xmax*0.1,ymax*0.85,['identified Kn =',num2str(kne,3)])
textt=[date,' : ',type,' : ',chname(chX,:),' / ',chname(chZ,:)];
title([textt])
grid

%----- calculating h
hn=0.5*imag(kCn)/kne;   % see Eq.(2.33)

fm1=0.2; fm2=0.4;  % These values are for model BS.
                   % Change this value, if you choose model NB.
[dum,ifm1]=min(abs(freq-fm1));
[dum,ifm2]=min(abs(freq-fm2));
Zh=hn(ifm1:ifm2)';
Hf=[freq(ifm1:ifm2)',ones(ifm2-ifm1+1,1)];   %  y=ax+b
Ph=[Hf'*Hf]\Hf'*Zh;
beta=Ph(2);    % identified beta

cn=imag(kCn);
Zc=cn(ifm1:ifm2)';
Pc=[Hf'*Hf]\Hf'*Zc;
Cne=Pc(1)/(2*pi);     % identified c

%----- plotting -- equivalent-h & viscous damping
subplot(3,1,2)
plot(freq(1:fmax),hn,'r-')   % Fig.2.4
xmax=1.0;
ymax=0.5;  % Max value for the Y-axis. This value is for model BS.
           % Change this value, if you choose model NB.
axis([0 xmax 0 ymax])
xlabel(['frequency (Hz)'])
ylabel(['hysteretic damping ratio'])
text(xmax*0.11,ymax*0.85,['identified beta = ',num2str(beta,3)])
grid

subplot(3,1,3)
plot(freq(1:fmax),cn,'r-')  % Fig.2.5
ymax=25;   % Max value for the Y-axis. This value is for model BS.
           % Change this value, if you choose model NB.
axis([0 xmax 0 ymax])
xlabel(['frequency (Hz)'])
ylabel(['imaginary part (tonf/cm)'])
text(xmax*0.11,ymax*0.85,['identified c =',num2str(Cne,3)])
grid
