%MSAC_NEW - Return a new SAC structure
%
%  [sac_trace_structure] = msac_new(timeseries,delta) ;
%
%     Returns a new SAC (time-series) structure with data x, and sampling delta.
%     This sets a default b and e of 0 and (nx-1)*delta respectively. All other
%     headers are set to SAC Null values. X must be a vector, delta a scalar.
%
%  [sac_trace_structure] = msac_new(timeseries,delta,HEADER,VALUE,...) ;
%     
%     Further arguments are interpreted as header, value pairs: e.g.,
%     [tr] = msac_new(x,delta,'evdp',600) sets the event depth header
%     to be 600 (kilometres). There (obviously) must be an even number of
%     these. 
%
%  NB. No 'smart' processing of header values is undertaken (e.g., setting b
%  or e will not set the other; setting event and station location will not
%  set gcarc etc etc). 
%
%  NB. Default byte-order is big-endian, even on a little endian machine. To
%      change, set the endian header to 'l'. 
%

%------------------------------------------------------------------------------%
%                                                                              %
%   (c) James Wookey, University of Bristol, 2006                              %
%                                                                              %
%   Version 1.0                                                                %
%                                                                              %
%   This program is free software; you can redistribute it and/or modify       %
%   it under the terms of the GNU General Public License as published by       %
%   the Free Software Foundation; either version 2 of the License, or          %
%   (at your option) any later version.                                        %
%                                                                              %
%   This program is distributed in the hope that it will be useful,            %
%   but WITHOUT ANY WARRANTY; without even the implied warranty of             %
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              %
%   GNU General Public License for more details.                               %
%                                                                              %
%   You should have received a copy of the GNU General Public License          %
%   along with this program; if not, write to:                                 %
%                                                                              %
%   Free Software Foundation, Inc., 59 Temple Place, Suite 330,                %
%   Boston, MA  02111-1307  USA                                                %
%                                                                              %
%   See http://www.gnu.org/licenses/gpl.txt for full GPL terms.                %
%                                                                              %
%------------------------------------------------------------------------------%
% Version history

% v1.0 - Original version
% 2006-06-15

function [tr] = msac_new(timeseries,delta,varargin) ;

% check input parameters
if nargin < 2 
   error('Not enough arguments for MSAC_NEW')
end
if ~isvector(timeseries)
   error('Argument timeseries must be a vector')
end
if ~isscalar(delta)   
   error('Argument delta must be a scalar')
end

%  setup SAC null values

   SAC_rnull = -12345.0      ;
   SAC_inull = -12345        ;
   SAC_lnull = -12345        ;
   SAC_cnull = '-12345  '    ;

   sachdrs = ...
          {'delta','depmin','depmax','scale','odelta','b','e','o','a',...
           'internal0','t0','t1','t2','t3','t4','t5','t6','t7','t8','t9',...
           'f','resp0','resp1','resp2','resp3','resp4','resp5','resp6',...
           'resp7','resp8','resp9','stla','stlo','stel','stdp','evla',...
           'evlo','evel','evdp','mag','user0','user1','user2','user3',...
           'user4','user5','user6','user7','user8','user9','dist','az',...
           'baz','gcarc','internal1','internal2','depmen','cmpaz',...
           'cmpinc','xminimum','xmaximum','yminimum','ymaximum',...
           'unused1','unused2','unused3','unused4','unused5','unused6',...
           'unused7','nzyear','nzjday','nzhour','nzmin','nzsec',...
           'nzmsec','nvhdr','norid','nevid','npts','internal3','nwfid',...
           'nxsize','nysize','unused8','iftype','idep','iztype',...
           'unused9','iinst','istreg','ievreg','ievtyp','iqual',...
           'isynth','imagtyp','imagsrc','unused10','unused11',...
           'unused12','unused13','unused14','unused15','unused16',...
           'unused17','leven','lpspol','lovrok','lcalda','unused18',...
           'kstnm','kevnm','khole','ko','ka','kt0','kt1','kt2','kt3',...
           'kt4','kt5','kt6','kt7','kt8','kt9','kf','kuser0','kuser1',...
           'kuser2','kcmpnm','knetwk','kdatrd','kinst'} ;
   
   for i=[1:70]
      eval(sprintf('tr.%s = SAC_rnull ;',sachdrs{i})) ;
   end
      
   tr.delta     = delta ;
   tr.b         = 0.0 ;

   for i=[71:105]
       eval(sprintf('tr.%s = SAC_inull ;',sachdrs{i})) ;
   end

   tr.nvhdr     = 6 ;% default
   tr.npts      = length(timeseries) ;% number of samples
   tr.iftype    = 1 ;% default (time series file)
   tr.idep      = 5 ;% default
   tr.iztype    = 9 ;% default
   tr.ievtyp    = 5 ;% default

   tr.leven     = 1 ;% default
   tr.lpspol    = 0 ;
   tr.lovrok    = 1 ;
   tr.lcalda    = 1 ;
   tr.unused18  = 0 ;

   for i=[111:133]
          eval(sprintf('tr.%s = SAC_cnull ;',sachdrs{i})) ;
   end

%  setup the .endian field (used for writing)
   tr.endian = 'b' ; % default

   tr.trace = timeseries ;
   
% now process the optional arguments (if there are any)
   if nargin == 2 % no, there aren't
      tr.time = [tr.b:delta:(tr.npts-1)*delta] ;
      return ;
   end
   
%  check that we have an even number of arguments
   if mod(length(varargin),2)~=0
      error('Header-Argument pairs must be specified') ;
   end
   
%  now process the list
   for i=1:2:length(varargin)-1

%     check that the header is valid
      tmp = cell(1,133) ;
      [tmp{:}] = deal(varargin{i}) ;
      index1 = strcmp(sachdrs,tmp) ;
      index2 = find(index1==1) ;
      if isempty(index2)
         error('Invalid header specified: %s',varargin{i}) ;
      end

%     check that the value specified is valid, and set it
%     REAL HEADERS
      if index2 >= 1 & index2 <= 70
         if ~isnumeric(varargin{i+1})
            error('Header %s must take a numeric value',varargin{i})
         end
         eval(sprintf('tr.%s = varargin{i+1} ;',varargin{i})) ;
         
%     INTEGER HEADERS
      elseif index2 >= 71 & index2 <= 105
         if ~isnumeric(varargin{i+1})
            error('Header %s must be an integer',varargin{i})
         end
         if mod(varargin{i+1},1)~=0
            error('Header %s must be an integer',varargin{i})
         end
         eval(sprintf('tr.%s = varargin{i+1} ;',varargin{i})) ;

%     LOGICAL HEADERS
      elseif index2 >= 106 & index2 <= 110
         if ~islogical(varargin{i+1})
            error('Header %s must take a logical (0 or 1) value',varargin{i})
         end
         eval(sprintf('tr.%s = varargin{i+1} ;',varargin{i})) ;
        
%     CHARACTER HEADERS
      elseif index2 > 110
         if ~ischar(varargin{i+1})
            error('Header %s must take a string value',varargin{i})
         end
         if index2 == 112
            ncharmax = 16 ;
         else
            ncharmax = 8 ;
         end
         
         if length(varargin{i+1}) > ncharmax
            eval(sprintf('tr.%s = varargin{i+1}(1:ncharmax) ;',varargin{i})) ;
            warning('Header value for %s was truncated',varargin{i}) ;
         else
            eval(sprintf('tr.%s = varargin{i+1} ;',varargin{i})) ;
         end
      end
   end
   
%  setup the .time field used by some MATLAB scripts   
   tr.time = [tr.b:delta:(tr.npts-1)*delta] ;


% end of MSAC_NEW

