readTSnims.m

function [y,SysTF,dt,t0,dataZero,ch,sta,error] = readTSnims(nimsFiles,decFac,SCALE)

% Usage: [y,SysTF,dt,t0,dataZero,ch,sta,error] = readTSnims(nimsFiles,decFac,SCALE);
%
%   Loads and aligns time series data for a series of stations;
%     modified (and simplified) from rdTS.m
%   In contrast to rdTS, this does not allow yet for timing errors
%     (with gps, these should not occur); also system TFs are not really
%     set co
%
%   Inputs: cell array giving pathnames for NIMS data files for

%
%   Outputs:
%	y = data array (all channels, for all sites merged into a
%		single array
%	TF = structure containing system parameter table for all 
%		channels at all sites (presently trivial ... instrument
%                responses are being ignored for NIMS so far)
%       dt = sampling interval
%	t0 = time of first data point in array y 
%		(relative to dataZero, in seconds)
%   dataZero = the reference date (string),
%            this is usually first day of the same year
%            e.g., datestr([2011 1 1 0 0 0]) = '01-Jan-2011'
%	ch = character string giving channel ids
%	sta(nchT) = site names (one for each channel), derived 
%            from data file names
%	error flag < 0 for error reading files

if nargin == 1
   decFac = 1;
   SCALE = 0;
end
if nargin  == 2
   SCALE = 0;
end

error = 0;
%  these are now hard-coded for NIMS
nch = 5;
freq = [8,1e-5];

nfiles=length(nimsFiles);
y=[];
ch=[];
samp=[];
STA={};
sta=[];
TF = [];

for k=1:nfiles
   %% Read NIMS asc or bin time series files
   if(nimsFiles{k}(end-2:end) == 'asc')
      [data,dt,dataStart,dataZero,ch_id,STA{k}] = ...
		readNIMSasc(nimsFiles{k},decFac,SCALE);
   elseif(nimsFiles{k}(end-2:end) == 'bin')
      %   binary NIMS file
      [data,metaData,dataZero] = ...
		readNIMSbin(nimsFiles{k},decFac,SCALE);
      dt = metaData.dt;
      dt = dt*decFac;
      startDate = metaData.startDate;
      startTime = metaData.startTime;
      t0 = datenum([startDate ' ' startTime],'yyyy-mm-dd HH:MM:SS');
      t0 = t0-datenum(dataZero);
      t0 = round(t0*24*3600/dt);
      ch_id = ['Hx';'Hy';'Hz';'Ex';'Ey'];
      STA{k} = metaData.runID;
   elseif(nimsFiles{k}(end-2:end) == 'dbn')
      %   decimated NIMS file
      [data,metaData,dataZero] = ...
		readNIMSbin(nimsFiles{k},decFac,SCALE);
      dt = metaData.dt;
      dt = dt*decFac;
      startDate = metaData.startDate;
      startTime = metaData.startTime;
      t0 = datenum([startDate ' ' startTime],'yyyy-mm-dd HH:MM:SS');
      t0 = t0-datenum(dataZero);
      t0 = round(t0*24*3600/dt);
      ch_id = ['Hx';'Hy';'Hz';'Ex';'Ey'];
      STA{k} = metaData.runID;
   elseif(nimsFiles{k}(end-2:end) == 'mat')
      %   matlab mat file produced by matlab NIMS reader
      [data,dataStart,dataZero,dt,ch_id,STA{k}] = ...
		readNIMSmat(nimsFiles{k},decFac,SCALE);
      %  compute relative starting scan number
      dt = dt*decFac;
      t0 = round((datenum(dataStart)-datenum(dataZero))*24*3600/dt);
   end
   %  here t0 is in samples, not seconds
   [SysTF] = setUnitSysTF(nch,freq); 
   ch=[ch; ch_id];
   samp=[samp;dt];
   %% Check to see if start times and length of files are equal
   [dum,npts] = size(data);
   if k==1
      nchT = nch;
      nptsMax = npts;
      t0All = t0;
      y=data;      
      numch = nchT;
   else
      nOffset = t0All - t0;
      if(nOffset > 0)
        y = [zeros(nchT,nOffset)/0 y];
        nptsMax = nptsMax + nOffset;
        t0All = t0;
      elseif (nOffset < 0)
        data = [zeros(nch,-nOffset)/0 data];
        npts = npts - nOffset;
      end
      nExtra = npts - nptsMax;
      if(nExtra > 0)
        y = [y zeros(nchT,nExtra)/0];
        nptsMax = nptsMax + nExtra;
      elseif (nExtra < 0)
        [NN1,idum] = size(data);
        data = [data zeros(NN1,-nExtra)/0 ];
      end
      y=[y;data];      
      nchT = nchT+nch;
      numch=[numch; nch];

      if samp(k)~=samp(1)
         error=-2;
         return
      end
   end
   TF = [TF ; SysTF.TF];
end
SysTF.TF = TF;

K = 0;
for k = 1:nfiles
   for ic = 1:numch(k)
      K = K+1;
      sta(K,:) = sprintf('%10s',char(STA{k}));
   end 
end
%   convert t0 to seconds
t0 = t0All*dt;