
clear
pause off 

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

dir

file1=input(' Give filename for x 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 z 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) ? ');
p=input(' Give dominant period  ? ');
src=ricker(dt,p);

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


nr=length(s1)/nt;
inc=input(' Receiver increment for plotting (>=1) ? ');
off=input(' Receiver offset  ? ');

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

t=(1:nt)*dt;

for i=1:inc:nr,
s1(:,i)=mconv(s1(:,i),src);
s2(:,i)=mconv(s2(:,i),src);
end


% find max
rmax1=max(max(abs(s1)));
rmax2=max(max(abs(s2)));
rmax=max([rmax1 rmax2 ]);

iinc=0;
s1i=zeros([nt,nr/inc]);
s2i=zeros([nt,nr/inc]);


% preparing seismograms for plots
for i=1:inc:nr,
iinc=iinc+1;
s1i(:,iinc)=s1(:,i)+(iinc-1)*rmax*off;
s2i(:,iinc)=s2(:,i)+(iinc-1)*rmax*off;
end

subplot(1,2,1),plot(t,s1i),title(' X-component'),xlabel(' Time (usec) '),
axis([0 max(t) -rmax*off nr/inc*rmax*off]),pause
subplot(1,2,2),plot(t,s2i),title(' Z-component'),xlabel(' Time (usec) ')
axis([0 max(t) -rmax*off nr/inc*rmax*off]),pause

ans=input(' Print plot to PS file (1) ? ');

if ans == 1,
disp(' Writing plot to file seis.eps ')
print -deps seis.eps
end

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