ename=strcat(station,'.EGE.NM.1.sac');  %construct the filename
nname=strcat(station,'.EGN.NM.1.sac');
zname=strcat(station,'.EGZ.NM.1.sac');
npts=ttot/delt + 1;
[t,e] = readsac(ename,npts);     % read the data
[t,n] = readsac(nname,npts);
[t,z] = readsac(zname,npts);
e=dmean(e);           % remove the mean from each
n=dmean(n);
z=dmean(z);
ntest=1;
if ntest == 1;
 else;
 % A synthetic test case
 npts=1000;
 ttot=delt*(npts-1);
 azrayl=30.;
 azlove=45;
 amprayl=2;
 amplove=1.;
 [t,e,n,z]=testchirps(npts,delt,amplove,amprayl,azlove,azrayl);
end;
if hilb ==1;     % hilbert transform the vertical component
 zh=hilbert(z);
 z=-imag(zh);
 else;
end;
f1=figure('name','DATA SEISMOGRAMS');  % plot the data
subplot(3,1,1);
plot(t,e);
xlabel('time sec');
ylabel(strcat('EW Comp at ',station));
subplot(3,1,2);
plot(t,n);
xlabel('time sec');
ylabel(strcat('NS Comp at ',station));
subplot(3,1,3);
plot(t,z);
xlabel('time sec');
ylabel(strcat('Z comp at ',station));
%
% filter the data
%
if flp > 0;
   e1=bandpass(e,flp,fhi,npts,delt);
   n1=bandpass(n,flp,fhi,npts,delt);
   z1=bandpass(z,flp,fhi,npts,delt);
   e=e1;
   n=n1;
   z=z1;
   %
 f2=figure('name','FILTERED SEISMOGRAMS');  % plot the filtered data
 subplot(3,1,1);
 plot(t,e);
 xlabel('time sec');
 ylabel(strcat('EW Comp at ',station));
 subplot(3,1,2);
 plot(t,n);
 xlabel('time sec');
 ylabel(strcat('NS Comp at ',station));
 subplot(3,1,3);
 plot(t,z);
 xlabel('time sec');
 ylabel(strcat('Z comp at ',station));
else;
end;

%
%   Moving window loop
%
npts1=fix(ttot/delt) + 1   %  total number of samples to analyze
nwin=fix(twin/delt) + 1    %  number of samples in a time window
npshift=fix(twin/(2*delt))+1  % number of samples to shift over
kfin=fix((npts1-nwin)/(npshift+1))+1 % number of time windows considered
mxde1=0.;
mxde2=0.;
mxde3=0.;
for k=1:kfin;
   nwinst=(k-1)*(npshift-1)+1;  % start of time window
   nwinfn=nwinst+nwin-1;    % end of time window
   a=csigm(e,n,z,nwinst,nwinfn);  % signal matrix
   %a
   c=a'*a;         % correlation matrix
   %c
   [v1,d1]=eig(c);      % eigenvalue/eigenvectors
   %d1
   %v1
   [v,d]=order(v1,d1);     % order the eigenvalues and eigenvectors
   %d
   %v
   ang1(k)=atan2(v(1,1),v(2,1)) * 180/pi; % azimuth for first  eigenvalue
   ang2(k)=atan2(v(1,2),v(2,2)) * 180/pi; % azimuth for second eigenvalue
   ang3(k)=atan2(v(1,3),v(2,3)) * 180/pi; % azimuth for third  eigenvalue
 de1(k)=d(1);
 de2(k)=d(2);
   de3(k)=d(3);
   mxde1=max(mxde1,de1(k));  % find the maximum values
   mxde2=max(mxde2,de2(k));
   mxde3=max(mxde3,de3(k));
 vang1(k)=acos(abs(v(3,1)))* 180/pi; %angle from the vertical
 vang2(k)=acos(abs(v(3,2)))* 180/pi;
 vang3(k)=acos(abs(v(3,3)))* 180/pi;
 t2(k)=delt*(nwinst-1);    % assign time for this window to the window start
end;
%
f3=figure('name','Eigenvalues and Inferred Azimuth');
subplot(3,1,1);
plot(t2,de1,'-or',t2,de2,'-dg',t2,de3,'-+b');
xlabel('time sec');
ylabel('eigenvalues');
subplot(3,1,2);
plot(t2,ang1,'-or',t2,ang2,'-dg',t2,ang3,'-+b');
xlabel('time sec');
ylabel('Azimuth ');
subplot(3,1,3);
plot(t2,vang1,'-or',t2,vang2,'-dg',t2,vang3,'-+b');
xlabel('time sec');
ylabel('incidence angle ');
%
%  Rose plot
f4=figure('name','Azimuth Distribution');
subplot(2,3,1);
title('Azimuth - Largest Eigenvalue');
rose(ang1*pi/180,100);
subplot(2,3,2);
title('Azimuth - Intermediate Eigenvalue');
rose(ang2*pi/180,100);
subplot(2,3,3);
title('Azimuth - Smallest Eigenvalue');
rose(ang3*pi/180,100);
%
%  Compass plots
%f5=figure('name','Compass Plots');
subplot(2,3,4);
title('Azimuth - Largest Eigenvalue');
compass(de1.*cos(ang1*pi/180),de1.*sin(ang1*pi/180));
subplot(2,3,5);
title('Azimuth - Intermediate Eigenvalue');
compass(de2.*cos(ang2*pi/180),de2.*sin(ang2*pi/180));
subplot(2,3,6);
title('Azimuth - Smallest Eigenvalue');
compass(de3.*cos(ang3*pi/180),de3.*sin(ang3*pi/180));
% Do the same thing for the highest amplitudes for each
% eigenvalue
%
nskip=1;
if nskip == 1;
   else;
neig1=0;
neig2=0;
neig3=0;
for k=1:kfin;
   if de1(k) >= 0.5*mxde1;
      neig1=neig1+1;
      angm1(neig1)=ang1(k);
   else;
   end;
   if de2(k) >= 0.5*mxde2;
      neig2=neig2+1;
      angm2(neig2)=ang2(k);
   else;
   end;
 if de3(k) >= 0.5*mxde3;
      neig3=neig3+1;
      angm3(neig3)=ang3(k);
   else;
   end;
end;
%
subplot(2,3,4);
title('Azimuth - Largest Eigenvalue,50% Threshold');
rose(angm1*pi/180,100);
subplot(2,3,5);
title('Azimuth - Intermediate Eigenvalue,50% Threshold');
rose(angm2*pi/180,100);
subplot(2,3,6);
title('Azimuth - Smallest Eigenvalue,50% Threshold');
rose(angm3*pi/180,100);
end;

Only the first part of the  filename is given.  This part reconstructs the rest.
 
 

readsac is a routine written by Roy Greenfield and Tom Battenhouse.
 

Remove the mean.
 
 
 
 

This was used when checking the program results.
 
 
 
 
 
 
 
 
 
 
 
 
 

Plot the raw data.  We are going to plot three seismograms in the window.
 
 
 
 
 
 
 
 
 
 
 

Filter the data between flp and fhi.
 
 
 
 
 
 
 

Plot the filtered seismograms.
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Apply a moving window analysis to the 3 component data.
 
 
 
 
 
 
 
 
 
 
 
 

Load a time portion of the signals into the signal matrix.

Construct the correlation matrix.

Find the eigenvalues and eigenvectors of the signal matrix.
 

Order the resulting eigenvectors/eigenvalues in decending order.
 

Calculate azimuths from the eigenvector components.
 
 
 
 
 
 
 
 

Calculate the apparent incidence angles from the eigenvectors.
 
 
 
 
 

Plot this stuff.
 
 
 
 
 
 
 
 
 
 
 
 

Now make a Rose diagram of azimuths.
The top set is directly from the eigenvectors.
 
 
 
 
 
 
 
 
 

This set of "compass" plots is weighted by the eigenvalues.
 
 
 
 
 
 
 
 
 
 
 
 
 

An experimental portion not invoked by the example shown here.