 %function [t,a,p]=readsac(sacfile,npts,'many') 
 % 
 %Required input:  sacfile 
 %Optional input:  npts,'many' 
 %Required output: t,a 
 %Optional output: p 
 % 
 %t = independent variable ; a = dependent variable. 
 %sacfile= file containing sac data (or list of files). 
 %   
 % 
 %The sac headers are stored in p as: 
 %p(1:3)     delta,depmin,depmax, 
 %p(4:6)     B,E,Event Origin Time 
 %p(7:16)    T0 - T9 
 %p(17:24)   STLA,LO,EL,DP,EVLA,LO,EL,DP 
 %p(25:28)   Dist,Az,BAz,GCArc 
 %p(29:34)   Year,Day,Hour,Min,Sec,MillSec 
 %p(35)      Npts 
 %p(36:40)   IFType,IDep,IZType,IevTyp,ISynth 
 %p(41)      user0 - number of days in a cross correlation
 % 
 %With npts presribed, only the first npts of the data 
 %file are read in.  
 % 
 %With the Option 'many' the filenames are read from 
 %the file specified. Npts must be specified because 
 %the files may have different npts. 
 % 
 %Ex: [t,a]=readsac('sac.data'); -read data from 1 binary sac file 
 %Ex: [t,a]=readsac('mysac.dat',2400,'many') - read several  
 %datafiles. (The list of files is in mysac.dat.) 
 % 
 %Written by Greenfield&Battenhouse 
   

 function [t,a,p]=readsac(sacfile,varargin) 

 if length(varargin)==0  
  [t,a,p]=readsacfile(sacfile); 
 end 

 if length(varargin)==1 
  [t,a,p]=readsacfile(sacfile); 
  npts=varargin{1};  
  t=t(1:npts); 
  a=a(1:npts); 
  p(35)=npts; 
  p(5)=(npts-1)*p(1); 
 end 

 if length(varargin)==2 
  npts=varargin{1}; 
  fid=fopen(sacfile,'r'); 
  j=1; while feof(fid)==0, file(j,:)=fgetl(fid); j=j+1; end;   
 [numfile,dull]=size(file); 
  fclose(fid);  
  for j=1:numfile 
   [t_,a_,p_]=readsacfile(file(j,:)); 
   p(j,:)=p_; 
   p(j,35)=npts; 
   p(j,5)=(npts-1)*p(1); 
   t(j,:)=t_(1:npts); 
   a(j,:)=a_(1:npts); 
  end 
 end 
 return 
   

 function [t,a,p]=readsacfile(sacfile) 
 sacfid=fopen(sacfile,'r'); 
 fparam=zeros(1,70);  iparam=zeros(1,35);  p=zeros(1,41); 
 % Read Parameters 
  fparam=fread(sacfid,70,'float32'); 
  iparam=fread(sacfid,35,'int32');
  p=assignparam(fparam,iparam); 
 for j=1:40, if p(j)==-12345 p(j)=nan; end, end; 
 % Read data 
  npts=p(35); 
  a=zeros(1,npts); 
  t=zeros(1,npts); 
  t=p(4):p(1):p(4)+(npts-1)*p(1); 
  fseek(sacfid,158*4,-1);  %(used to be 154?) 
  a(1:npts)=fread(sacfid,npts,'float32'); 
 if isNAN(p(4)) | isNAN(p(1)), 
  t(1:npts)=fread(sacfid,npts,'float32'); 
 end 
 fclose(sacfid); 
 return 

 function p=assignparam(fparam,iparam) 
 p=zeros(1,40); 
 p(1:3)=fparam(1:3);            %delta,depmin,depmax, 
 p(4:5)=fparam(6:7);            %B,E 
 p(6)=fparam(8);                %Event Origin Time 
 p(7:16)=fparam(11:20);         %T0 - T9 
 p(17:24)=fparam(32:39);        %STLA,LO,EL,DP,EVLA,LO,EL,DP 
 p(25:28)=fparam(51:54);        %Dist,Az,BAz,GCArc 
 p(29:34)=iparam(1:6);          %Year,Day,Hour,Min,Sec,MillSec 
 p(35)=iparam(10);              %Npts 
 p(36:38)=iparam(16:18);        %IFType,IDep,IZType 
 p(39:40)=[iparam(23) iparam(25)];   %IevTyp,ISynth 
 p(41) = fparam(41);            %user0 
 return 

 function l=isNAN(a) 
 l=1- (a>0 | a<0 | a==0); 
 return 
