
clear
pause off 

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

dir

file=input(' Give name for run : ');

disp(' Which profile ? ');
disp(' 1) X '); 
disp(' 2) Y ') ;
disp(' 3) Z ') ;
ans=input(' Which ? ');

if ans == 1 
fileout=sprintf('%s_X3c.eps',file);
file1=sprintf('%s_Xx',file);
file2=sprintf('%s_Xy',file);
file3=sprintf('%s_Xz',file);
elseif ans == 2,
fileout=sprintf('%s_Y3c.eps',file);
file1=sprintf('%s_Yx',file);
file2=sprintf('%s_Yy',file);
file3=sprintf('%s_Yz',file);
elseif ans == 3,
fileout=sprintf('%s_Z3c.eps',file);
file1=sprintf('%s_Zx',file);
file2=sprintf('%s_Zy',file);
file3=sprintf('%s_Zz',file);
end


if ~exist(file1), disp(' File (x) not found '),break, end
disp(' Loading x-component ')
eval(['load -ascii ' file1]);
if ~exist(file2), disp(' File (y) not found '),break, end
disp(' Loading y-component ')
eval(['load -ascii ' file2]);
if ~exist(file3), disp(' File (z) not found '),break, end
disp(' Loading z-component ')
eval(['load -ascii ' file3]);

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

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

nr=length(s1)/nt;

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


again=1;

while again==1,


inc=input(' Receiver increment for plotting (>=1) ? ');
off=input(' Receiver offset  ? ');

t=(1:nt)*dt;

disp(' Seismogram plots : ');
disp(' 1) Trace normalised ? ');
disp(' 2) Relative scaling ? ');
ans=input(' Which ? ');

if ans==2,

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

iinc=0;
s1i=zeros([nt,nr/inc]);
s2i=zeros([nt,nr/inc]);
s3i=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;
s3i(:,iinc)=s3(:,i)+(iinc-1)*rmax*off;
end

elseif ans == 1,

% find max
rmax1=max(max(abs(s1)));
rmax2=max(max(abs(s2)));
rmax3=max(max(abs(s3)));
rmax=max([rmax1 rmax2 rmax3 ]);
 
iinc=0;
s1i=zeros([nt,nr/inc]);
s2i=zeros([nt,nr/inc]);
s3i=zeros([nt,nr/inc]);
 
% preparing seismograms for plots
for i=1:inc:nr,
iinc=iinc+1;
s1i(:,iinc)=s1(:,i)/max(s1(:,i))*rmax+(iinc-1)*rmax*off;
s2i(:,iinc)=s2(:,i)/max(s2(:,i))*rmax+(iinc-1)*rmax*off;
s3i(:,iinc)=s3(:,i)/max(s3(:,i))*rmax+(iinc-1)*rmax*off;
end

end





subplot(1,3,1),plot(t,s1i),
title(file1),xlabel(' Time (sec) '),
axis([0 max(t) -rmax*off nr/inc*rmax*off]),pause
subplot(1,3,2),plot(t,s2i),
title(file2),xlabel(' Time (sec) ')
axis([0 max(t) -rmax*off nr/inc*rmax*off]),pause
subplot(1,3,3),plot(t,s3i),
title(file3),xlabel(' Time (sec) ')
axis([0 max(t) -rmax*off nr/inc*rmax*off]),pause




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

if rans == 1,
disp(sprintf(' Writing plot to file %s',fileout));
eval([' print -deps ' fileout ]);
end

again=0;
again=input(' Plot again (1-yes) ');

end

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