%GET_SAC
%
%	get_sac
%
%	MATLAB M-file to read data directly from a SAC file.
%	Once the file is selected, it brings the following
%	data in the workspace:
%		npts:	number of data points
%       data:   time series data from sac file (double)
%		sampdelta:	sampling interval
%		date:	year month and day of the begining of file
%		hour:		time of
%		minu:		begining
%		sec:		of file
%		stat:   station name
%		<station name>:	trace data
%
%	Interactive form.  See also fget_sac for programable form.
%
%	                                     Guy Tytgat, 17 feb 95
%   Modified Bob Smalley, Nov 2009 to add default file select window name
%      if SAC_file_prompt does not exsit, and to determine endian automatically.
%      SAC header is 158 4 byte words long.

%[file,path] = uigetfile('*','Select a SAC file',300,300);
%[file,path] = uigetfile('*','Select a SAC file');
if(exist('SAC_file_prompt','var'))
    [file,path] = uigetfile('*',SAC_file_prompt);
else
    [file,path] = uigetfile('*','Enter name of SAC file');
end

filename = sprintf('%s%s',path,file);

%get file size
D = dir(filename);
sizeoffile=D.bytes/4;

%careful here -- if data comes from iris dmc - use unix
%if made on pc-linux (eg output sac on pc-linux) use intel!!
%[fid,message] = fopen(filename,'r','b');   %unix
%[fid,message] = fopen(filename,'r','l'); %intel
%start by assuming endian of data is same as computer
[comptype,maxsizd,endian]=computer;
%but in end, need to match to data, not machine where running.
[fid,message] = fopen(filename,'r',endian);

if fid ~= -1
  [fheader,count] = fread(fid,70,'float');
  [nheader,count] = fread(fid,35,'long');
  npts = nheader(10);
%check size from header against size of file - if same we are OK, if not
%change endian of data file and re-read.
  if npts+158~=sizeoffile
      %switch endian
      if endian=='l'||endian=='L',endian='b';else,endian='l';end;
      [fid,message] = fopen(filename,'r',endian);
      [fheader,count] = fread(fid,70,'float');
      [nheader,count] = fread(fid,35,'long');
      npts = nheader(10);
  end

  [lheader,count] = fread(fid,5,'long');
  [kheader,count] = fread(fid,[8,24],'char');

  sampdelta = fheader(1);
  date = nheader(1:2);
  hour = nheader(3);
  minu = nheader(4);
  seco = nheader(5) + (nheader(6)/1000);

  data = zeros(npts,1);
  [data,count] = fread(fid,npts,'float');

% added by bodin 19June98 to enable the reading of
% the file name in the calling routine, it will be
% called "fname"......
  fname=filename;
%
  close_stat = fclose(fid);
 % clear fheader nheader lheader kheader station
 % clear file path filename fid message count close_stat
else
  disp(message)
  clear file path filename fid message
end
