
clear
close all

pause off 

disp(' ')
disp(' ------------------------------------------ ')
disp(' Basic plotting for 2-Component seismograms ')
disp(' ------------------------------------------ ')
disp(' ')

dir

file1='pfo_t';
file2='pfo_exx';
% file1=input(' Give filename for z component (in single quotes) ');
% if ~exist(file1), disp(' File (x) not found '),break, end
% disp(' Loading x-component ')
eval(['load -ascii ' file1]);
% file2=input(' Give filename for exz component (in single quotes) ');
% if ~exist(file2), disp(' File (z) not found '),break, end
% disp(' Loading z-component ')
eval(['load -ascii ' file2]);



%nt=input(' Number of time samples (nt) ? ');
%dt=input(' Time increment (dt) in micros ? ');
%nr=input(' Number of receivers (nr) ? ');

nt=20000;
dt=.15;

p=input(' Give cut-off period  ? ');
%src=ricker(dt,p);

s1=eval([file1]);
s2=eval([file2]);


nr=length(s1)/nt;


s1=reshape(s1',nt,nr); 
s2=reshape(s2',nt,nr); 

t=(1:nt)*dt;

% tapering and filtering

ntap=500;
nord=4;
for i=1:nr,
s1(:,i)=taper(s1(:,i),ntap);
s1(:,i)=taper(s1(:,i),-ntap);
s1(:,i)=lowpass(s1(:,i),dt,1/p,nord);
s2(:,i)=taper(s2(:,i),ntap);
s2(:,i)=taper(s2(:,i),-ntap);
s2(:,i)=lowpass(s2(:,i),dt,1/p,nord);
end

dtheta=5;  % receiver distance

for ii=1:nr,

off=5+(ii-1)*dtheta;  % epicentral distance
vz=s1(:,ii);
exz=-int(s2(:,ii))*dt;
%exz=s2(:,ii);
fac=max(vz)/max(exz);

% normalize
vz=vz/max(vz);
exz=exz/max(exz);

plot(t,vz*dtheta+5+(ii-1)*dtheta,'b-',t,exz*dtheta+5+(ii-1)*dtheta,'r-'), hold on 
title(' PREM: black - vx; red - \epsilon_{xx} ')
ylabel(' Epicentral distance (deg)')
xlabel(' Time (s) ')
text(50,6+(ii-1)*dtheta,sprintf('c = %g m/s',fac),'FontSize',8)
axis([0 3000 -Inf  Inf])

end

hold off


% spectral ratios

% Loop over all stations (whole time window)


for i=1:nr,

% spectrum for strains
in=-int(s2(:,i))*dt; n=length(in); 
[fs,f]=myspectrum(in',dt);
sfs=abs(fs(1:length(f)));
% spectrum for translations
in=s1(:,i);n=length(in);
[fs,f]=myspectrum(in,dt);
tfs=abs(fs(1:length(f)));

if i==1, asfs=sfs; else asfs=asfs+sfs; end
if i==1, atfs=tfs; else atfs=atfs+tfs; end

end

% the ratio
sr=atfs./asfs;
Tmax=200;
Tmin=40;

figure
plot(1./f,sr),
xlabel(' Period (s) ')
ylabel(' Velocity (m/s) ')
axis([Tmin Tmax -Inf Inf])


disp(' ... bye ...')
