function [data,fs,time,Aux]=Acousonderead(start,duration,datapath);
% Function [data,fs,time,Aux]=Acousonderead(start_time,duration,datapath);
%
% This function will read requested data, given a start time & duration,
% from a directory of continuously sampled, single channel Acousonde 
% MT data files, spanning raw data files as necessary to build the 
% requested data vectors.
%
% Multi-channel LP[HF] combination, or duty cycle recordings are not yet
% supported.
%
% NOTE: Assumes all acoustic data *.MT files are co-located at the 'datapath'
% folder, and Auxiliary channel data files are in a sub-directory of the
% acoustic data, named 'Auxiliary'.
%
% USAGE: [data,fs,time,Aux]=Acousonderead(start_time,duration,datapath);
% INPUTS:
%               start_time: Start time of requested data.  Can be input in
%                           MatLab datenumber format, or as an ASCII formatted 
%                           string 'DD-MMM-YYYY hh:mm:ss' which will be
%                           converted to datenumber automatically
%               duration:   length of requested data (in seconds)
%               datapath:   (opt) string containing directory path, if the
%                           Acousonde data is not in the current directory
% OUTPUTS:
%               data:       Requested data vector, uPa units (1e-6 Pascals)
%               fs:         Sample rate of acoustic data (Hz)
%               time:       Time vector for acoustic data, in MatLab datenumber format
%               Aux:        Auxiliary channel structure data (optional)
% 
% Example:
%   [data,fs,time]=Acousonderead(datenum(2012,1,24,17,0,0),60,'C:\data\');
%       Will read 1 minute of data, starting at 1/24/12 17:00:00 GMT from
%       Acousonde files located in the 'C:\data\' directory
%
% ver 2.1
% Chris Miller
% Naval Postgraduate School
%
% sub-function MTRead.m courtesy of Bill Burgess, Acoustimetrics

% 8/21/09 Version 1.0 - Chris Miller, Naval Postgraduate School
% 2/27/12 Version 1.1 - Added remote path as optional input parameter (cm)
% 6/27/12 Version 1.2 - Added Auxiliary channel data as optional output
% 6/25/13 Version 2.0 - Added full Auxiliary channel extraction for start/duration
% 7/03/13 Version 2.0.1 - CM, corrected concatenation beween files for correct timing
% 9/04/13 Version 2.0.1a - CM, fixed bug where requested data stradled multiple MT files
% 6/19/18 Version 2.1 - Corrected mixed sample rates for Aux channels and missing AUX files (if not desired during backup)

DEBUG=0;
if nargin<3 || isempty(datapath)
    datapath=pwd;
    datapath=[datapath,filesep];
end
if ~isnumeric(start)
    start=datenum(start);
end
%% Auxiliary channel extraction
if nargout>3
    D2=dir([datapath,filesep,'Aux*']);
    Auxfiles=dir([datapath,filesep,D2.name,filesep,'*.MT']);
    % First pull temperature data files from Auxiliary folder
    names=struct2cell(Auxfiles);
    names2=cell2mat(names(1,:).');      % Auxiliary channel files present (char array)
    %% pull Temperature data
    ind=find(names2(:,3)=='T');
    if ~isempty(ind)
    T=[];
    for i=1:length(ind),
        [Temp,header(i),info(i)]=MTRead([datapath,filesep,D2.name,filesep,names2(ind(i),:)]);
        T=[T; Temp];    % Concatenate aux data from multiple files
    end;
    taux=(0:length(T)-1)/info(1).srate/3600/24+info(1).datenumber;  % time vector at channel sample rate
    index=find(taux>=start & taux<=start+duration/24/3600);
    Aux.Temp.t=taux(index);
    Aux.Temp.T=T(index);
    Aux.units.Temp=header(1).calunits;
    end;
    
    %% Pull pressure data
    ind=find(names2(:,3)=='P');
    if ~isempty(ind)
    P=[];
    for i=1:length(ind)
        [Pres,header(i),info(i)]=MTRead([datapath,filesep,D2.name,filesep,names2(ind(i),:)]);
        P=[P; Pres];    % Concatenate aux data from multiple files
    end
    paux=(0:length(P)-1)/info(1).srate/3600/24+info(1).datenumber;  % time vector at channel sample rate
    index=find(paux>=start & paux<=start+duration/24/3600);
    Aux.Pressure.t=paux(index);
    Aux.Pressure.P=P(index);
    Aux.units.Pressure=header(1).calunits;
    end
    
    %% Pull compass data (x,y,z) = (Pos North, Pos East, Pos Down)
    % X (North)
    ind=find(names2(:,3)=='X');
    if ~isempty(ind)
    X=[];
    for i=1:length(ind)
        [Compx,header(i),info(i)]=MTRead([datapath,filesep,D2.name,filesep,names2(ind(i),:)]);
        X=[X; Compx];    % Concatenate aux data from multiple files
    end
    Xaux=(0:length(X)-1)/info(1).srate/3600/24+info(1).datenumber;  % time vector at channel sample rate
    index=find(Xaux>=start & Xaux<=start+duration/24/3600);
    Aux.compass.t=Xaux(index);
    Aux.compass.N=X(index);
    Aux.units.compass.N=header(1).calunits;
    % Y (East)
    ind=find(names2(:,3)=='Y');
    Y=[];
    for i=1:length(ind)
        [Compy,header(i),info(i)]=MTRead([datapath,filesep,D2.name,filesep,names2(ind(i),:)]);
        Y=[Y; Compy];    % Concatenate aux data from multiple files
    end
    Yaux=(0:length(Y)-1)/info(1).srate/3600/24+info(1).datenumber;  % time vector at channel sample rate
    index=find(Yaux>=start & Yaux<=start+duration/24/3600);
    Aux.compass.W=Y(index);
    Aux.units.compass.W=header(1).calunits;
    % Z (Down)
    ind=find(names2(:,3)=='Z');
    Z=[];
    for i=1:length(ind)
        [Compz,header(i),info(i)]=MTRead([datapath,filesep,D2.name,filesep,names2(ind(i),:)]);
        Z=[Z; Compz];    % Concatenate aux data from multiple files
    end
    Zaux=(0:length(Z)-1)/info(1).srate/3600/24+info(1).datenumber;  % time vector at channel sample rate
    index=find(Zaux>=start & Zaux<=start+duration/24/3600);
    Aux.compass.Z=Z(index);
    Aux.units.compass.Z=header(1).calunits;
    end;

    %% Pull accelerometer data (i,j,k) = (Pos North, Pos East, Pos Down)
    % I (along axis)
    ind=find(names2(:,3)=='I');
    if ~isempty(ind)
    I=[];
    for i=1:length(ind)
        [Accelx,header(i),info(i)]=MTRead([datapath,filesep,D2.name,filesep,names2(ind(i),:)]);
        I=[I; Accelx];    % Concatenate aux data from multiple files
    end
    Iaux=(0:length(I)-1)/info(1).srate/3600/24+info(1).datenumber;  % time vector at channel sample rate
    index=find(Iaux>=start & Iaux<=start+duration/24/3600);
    Aux.accelerometer.t=Iaux(index);
    Aux.accelerometer.X=I(index);
    Aux.units.accelerometer.X=header(1).calunits;
    % J (across axis)
    ind=find(names2(:,3)=='J');
    J=[];
    for i=1:length(ind)
        [Accely,header(i),info(i)]=MTRead([datapath,filesep,D2.name,filesep,names2(ind(i),:)]);
        J=[J; Accely];    % Concatenate aux data from multiple files
    end
    Jaux=(0:length(J)-1)/info(1).srate/3600/24+info(1).datenumber;  % time vector at channel sample rate
    index=find(Jaux>=start & Jaux<=start+duration/24/3600);
    Aux.accelerometer.Y=J(index);
    Aux.units.accelerometer.Y=header(1).calunits;
    % Z (Down)
    ind=find(names2(:,3)=='K');
    K=[];
    for i=1:length(ind)
        [Accelz,header(i),info(i)]=MTRead([datapath,filesep,D2.name,filesep,names2(ind(i),:)]);
        K=[K; Accelz];    % Concatenate aux data from multiple files
    end
    Kaux=(0:length(K)-1)/info(1).srate/3600/24+info(1).datenumber;  % time vector at channel sample rate
    index=find(Kaux>=start & Kaux<=start+duration/24/3600);
    Aux.accelerometer.Z=K(index);
    Aux.units.accelerometer.Z=header(1).calunits;
    end

    %% Pull ambient Light sensor data - Model B units only
    ind=find(names2(:,3)=='L');
    if ~isempty(ind)
    Light=[];
    for i=1:length(ind)
        [L,header(i),info(i)]=MTRead([datapath,filesep,D2.name,filesep,names2(ind(i),:)]);
        Light=[Light; L];    % Concatenate aux data from multiple files
    end
    Laux=(0:length(Light)-1)/info(1).srate/3600/24+info(1).datenumber;  % time vector at channel sample rate
    index=find(Laux>=start & Laux<=start+duration/24/3600);
    Aux.Light.t=Laux(index);
    Aux.Light.Light=Light(index);
    Aux.units.Light=header(1).calunits;
    end
    %% Pull Battery engineering data, when sampled
    ind=find(names2(:,3)=='B');
    if ~isempty(ind)
    Bat=[];
    for i=1:length(ind)
        [L,header(i),info(i)]=MTRead([datapath,filesep,D2.name,filesep,names2(ind(i),:)]);
        Bat=[Bat; L];    % Concatenate aux data from multiple files
    end
    Laux=(0:length(Bat)-1)/info(1).srate/3600/24+info(1).datenumber;  % time vector at channel sample rate
    index=find(Laux>=start & Laux<=start+duration/24/3600);
    Aux.Bat.t=Laux(index);
    Aux.Bat.B=Bat(index);
    Aux.units.Bat=header(1).calunits;
    end
end
%% Acoustic data extraction
D=dir([datapath,'*.MT']);  % Get a directory list of all MT data files in current directory

% Test that the requested start lies within the available data
[~,~,infoF]=MTRead([datapath,D(1).name],0,0);      % read info on first file
[~,~,infoL]=MTRead([datapath,D(length(D)).name],0,0);
if start < infoF.datenumber || start>infoL.datenumber + infoL.seconds/24/3600,
    error('Acousonderead: Start_time must be between %s and %s',datestr(infoF.datenumber),datestr(infoL.datenumber+infoL.seconds/24/3600));
end

files = whichMT(datapath,D,start,duration);     % Which files contain requested data

fs=infoF.srate;
filetime=infoF.seconds;   % length of time for file (s)

[~,~,info1]=MTRead([datapath,D(files(1)).name],0,0);

if DEBUG
    disp(['Reading from ',info1.filename]);
end

if length(files)==1
    % File contains all requested data, proceed with single MTRead.m call
    [data,~,info]=MTRead([datapath,D(files(1)).name],(start-info1.datenumber)*24*3600,duration);
    if DEBUG
        disp(['Read from ',info.filename,', ',datestr(info.datenumber)]);
    end
    t=info.datenumber;
else
    % Need to keep track of successive file reads and combine output data
    [data,~,info]=MTRead([datapath,D(files(1)).name],(start-info1.datenumber)*24*3600);
    if DEBUG
        disp(['Read ',num2str(info.count),' from ',info.filename,', ',datestr(info.datenumber)]);
    end
    t=info.datenumber;
    duration=duration-info.count/info.srate; % remove read data from remaining 'duration' to read
    for fileno=2:length(files)-1,
        [datanext,~,info]=MTRead([datapath,D(files(fileno)).name]);    % mPa units
        if DEBUG
            disp(['Read from ',info.filename,', ',datestr(info.datenumber)]);
        end
        duration=duration-info.count/info.srate;
        data=[data; datanext];
    end
    [datanext,~,info]=MTRead([datapath,D(files(length(files))).name],0,duration); % mPa units
    if DEBUG
        disp(['Read ',num2str(info.count),' from ',info.filename,', ',datestr(info.datenumber)]);
    end
    data=[data; datanext];  % mPa units from MTRead function
end

data=data*1e3;      % uPa units  (convert MTread mPa units to uPa)
time=(0:length(data)-1)/fs/3600/24+t;
time=time(:);
%data=detrend(data); % De-mean the data  - leave for the end user to do if desired
%% plotting used to confirm data reads during debug evaluation
if DEBUG==2
    figure(1);clf;
    if nargin>2
      subplot(511);
      plot(Aux.Temp.t,Aux.Temp.T);datetick('x');axis tight;grid on;
      ylabel([char(176),'C']);title('Temperature');
      subplot(512);
      plot(Aux.Pressure.t,Aux.Pressure.P);datetick('x');datetick('x');axis tight;grid on;
      ylabel(['Pressure, (',Aux.units.Pressure,')']);title('Pressure');axis ij;
      subplot(513);
      plot(Aux.compass.t,Aux.compass.N,Aux.compass.t,Aux.compass.W,Aux.compass.t,Aux.compass.Z);legend('N','W','Z');
      title(['Compass bearing (',Aux.units.compass.N,')']);datetick('x');axis tight;grid on;
      subplot(514);
      plot(Aux.accelerometer.t,Aux.accelerometer.X,Aux.accelerometer.t,Aux.accelerometer.Y,Aux.accelerometer.t,Aux.accelerometer.Z);legend('X','Y','Z');
      title(['Accelerometer (',Aux.units.accelerometer.X,')']);datetick('x');axis tight;grid on;
      subplot(515);
      plot(time,data);datetick('x');axis tight;grid on;
      title('Acoustic data');ylabel('\muPa');
    else
      subplot(211);
      plot((0:length(data)-1)/info.srate/3600/24+t,data);
      datetick('x');axis tight;grid on;xlabel('time (HH:MM)');
      title(datestr(t(1)));
      subplot(212);
      specgram(data,2048,fs);
      v=axis;v(4)=10e3;axis(v);
    end
end

%%
function files = whichMT(datapath,D,start,duration)
% function files = whichMT(datapath,D,start,duration)
% Exhaustive search function that will search all MT files in a given path
% to find those fiiles that contain requested data.  Exhaustive search was
% implemented in order to return any requested data, including "ping pong"
% sampling between HF & LF channels
ending=start+duration/3600/24;
starts=zeros(1,length(D));
ends=zeros(1,length(D));

for i=1:length(D)
    [~,~,info]=MTRead([datapath,D(i).name],0,0);
    starts(i)=info.datenumber;
    ends(i)=info.datenumber+(info.seconds)/24/3600;
end
files=min(find(start<ends)) : min(find(start+duration/24/3600<ends));

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [p,header,info] = MTRead(filename,start,len,slMode)
% function [data,HEADER,INFO] = MTRead(FILENAME,START,LENGTH,'P')
% Reads a file that is in MT format.
%
%  INPUTS:
%           FILENAME	   the name of file (as a string variable) to read.
%           START          the number of seconds into the file to begin reading.
%                            (Default = beginning of file)
%           LENGTH         the length of data to read (in seconds)
%                            (Default = to end of file)
%           'P' (optional) Use if START and LENGTH are input using number of
%                            sample points rather than seconds.  When not
%                            present, time input (seconds) are assumed.
%
% OUTPUTS:
%           data           vector containing the data read from the file in 
%                            mPa (1e-3Pa) units.
%           HEADER         Structure variable containing the header information 
%                            (string values) from the MT file.  Information is
%                            primarily engineering values
%           INFO           Structure variable containing data relavent
%                            information
%           INFO.filename    Name of file
%           INFO.filesize    Total size of file, in bytes
%           INFO.when        Date/Time in MATLAB 6-element form
%           INFO.datenumber  Date/Time in MatLab datenumber format
%           INFO.yd          Date/Time in yearday 
%           INFO.whenC       Date/Time in C time_t form
%           INFO.srate       Sample Rate in Hz
%           INFO.nsamp       Number of sample points in file
%           INFO.seconds     Number of seconds in file
%           INFO.count       Number of sample points actually read and
%                            returned in data
%
%  Note that the time reported in HEADER is at the beginning
%  of the file, but that the time reported in INFO.when and
%  INFO.whenC are at the beginning of data read; in other
%  words, the time at the start of the file, plus START seconds.

%  MODIFICATION HISTORY:
%
%  8/18/03 WCB Original Bob Blaylock internal code patched
%  to work properly with strings that may have garbage after
%  the null terminator.
%
%  4/11/04 WCB Renamed to MTRead(), all necessary functions
%  incorporated into a single file, removed Blaylock comments
%  warning about "obsolete" header data -- these only apply
%  to the BGB extensions to MT format that the Bioacoustic
%  Probe does not use.  The necessary functions that had to
%  be appended to this file were:  c2mat_tm.m, mat2c_tm.m,
%  contains.m, limits.m, and strtrim.m.
%
%  4/14/04 WCB Added check for malformed wordsize entry in the
%  MT header; assumes 2 bytes per sample if the original value
%  is not properly formed.
%
%  9/17/09 CWM changed info.when assignment to correct an assignment error,
%  added info.datenumber, info.yearday time entry

if nargin < 4
    slMode = 's';
    if nargin < 3
        len = inf;
        if nargin == 1
            start = 0;
        end
    end
end

[f,msg] = fopen(filename,'r','b');
if f < 1
   fclose all;
   error([10 '  Sorry about this...' 10 '    But I can''t open this file.' 10 '      I get this error:' 10 10 '  "' filename '":  ' msg 10 10]);
end
header.magicstring = MakeString(fread(f,8,'char'));
if ~strcmp(header.magicstring,'DATA')
   fclose all;
   error([10 '  This is the wrong file!' 10 '    It''s not MacTrimpi format.' 10 '      I can do nothing!' 10 10 '  Trying to read "' filename '"' 10]);
end
header.totalhdrs       = MakeString(fread(f, 3,'char'));
header.abbrev          = MakeString(fread(f, 8,'char'));
header.stationcode     = MakeString(fread(f, 3,'char'));
header.title           = MakeString(fread(f,82,'char'));
header.month           = MakeString(fread(f, 3,'char'));
header.day             = MakeString(fread(f, 3,'char'));
header.year            = MakeString(fread(f, 5,'char'));
header.hours           = MakeString(fread(f, 3,'char'));
header.minutes         = MakeString(fread(f, 3,'char'));
header.seconds         = MakeString(fread(f, 3,'char'));
header.msec            = MakeString(fread(f, 4,'char'));
header.sampling_period = MakeString(fread(f,15,'char'));
header.samplebits      = MakeString(fread(f, 3,'char'));
header.wordsize        = MakeString(fread(f, 2,'char'));

if str2num(header.wordsize) < (str2num(header.samplebits)/8)
    warning(['  (Loading file "' filename '")' newl newl '  The samplebits field...' newl '    Does not fit the wordsize field.' newl '      This file may be bad. ' newl]);
end

header.typemark        = MakeString(fread(f, 1,'char'));
header.swapping        = MakeString(fread(f, 1,'char'));
header.signing         = MakeString(fread(f, 1,'char'));
header.caltype         = MakeString(fread(f, 1,'char'));
header.calmin          = MakeString(fread(f,15,'char'));
header.calmax          = MakeString(fread(f,15,'char'));
header.calunits        = MakeString(fread(f,40,'char'));
header.recordsize      = MakeString(fread(f, 6,'char'));
header.sourcevers      = MakeString(fread(f, 9,'char'));
fseek(f,0,'eof');
info.filename = filename;
info.filesize = ftell(f);
fclose(f);
info.srate = 1/str2num(header.sampling_period);
info.when = [str2num(header.year) str2num(header.month) str2num(header.day) str2num(header.hours) str2num(header.minutes) str2num(header.seconds)+str2num(header.msec)/1000];
info.datenumber = datenum([header.year,'/',header.month,'/',header.day,' ',header.hours,':',header.minutes,':',header.seconds,'.',header.msec],'yyyy/mm/dd HH:MM:SS.FFF');
info.yd = info.datenumber - datenum([info.when(1) 1 0]);
if upper(slMode) == 'P',      % Start & Length specified in # Points (samples)
     info.whenC = mat2c_tm(info.when) + start/info.srate;
     info.datenumber = info.datenumber + start/info.srate/24/3600;
else
     info.whenC = mat2c_tm(info.when) + start;        % Corrected start time (with offset)
     info.datenumber = info.datenumber + start/24/3600;
end
info.when = c2mat_tm(info.whenC);
%
% Some MT files have corrupted wordsize (no null terminator) -- if so
% assume 2-byte words
%
if (isempty(header.wordsize))
 	header.wordsize = '2';
end
info.nsamp = (info.filesize - 512*str2num(header.totalhdrs))/str2num(header.wordsize);
info.seconds = info.nsamp/info.srate;

if len > 0  %  Only load data if it's been asked for.
    if any(contains(header.swapping,'SLsl'))
        mode = 'ieee-le';
    else
        mode = 'ieee-be';
    end
    [f,msg] = fopen(filename,'rb',mode);
    if f < 1
        fclose all;
        error([10 '  Sorry about this...' 10 '    But I can''t open this file.' 10 '      I get this error:' 10 10 '  "' filename '":  ' msg 10 10]);
    end

    if upper(slMode) == 'P'    % specified start time in sample 'P'oints rather than time
        status = fseek(f,round(512*str2num(header.totalhdrs) + round(start)*str2num(header.wordsize)),'bof');  % Skip by samples/points
    else
        status = fseek(f,round(512*str2num(header.totalhdrs) + round(start*info.srate)*str2num(header.wordsize)),'bof');   % skip by time (seconds)
    end
    if status == 0  %% If status is nonzero, we probably went past the end of the file.
        if any(header.caltype == 'fF')
            if ~any(str2num(header.wordsize) == [4,8])
                fclose(f);
                error([10 '  Invalid word size!' 10 '    Only valid Float sizes...' 10 '      Are four or eight bytes.' 10 10]);
            end
            binType = ['float' num2str(str2num(header.wordsize)*8)];
        else
            binType = ['bit' num2str(str2num(header.wordsize)*8)];
            if any(contains(header.signing,'Uu'))
                binType = ['u' binType];
            end
        end
        if upper(slMode) == 'P'
            [p,info.count] = fread(f,len,binType);
        else
            [p,info.count] = fread(f,round(len*info.srate),binType);
        end
        fclose(f);
        
        calmax = str2num(header.calmax);
        calmin = str2num(header.calmin);
        if (length(calmin) == 1) && (length(calmax) == 1) && ((calmin + eps) < calmax) && ~any(header.caltype == 'fF')
            calmax = str2num(header.calmax);
            calmin = str2num(header.calmin);
            if any(contains(header.signing,'Uu'))
                bitmin = 0;
                bitmax = (2^str2num(header.samplebits)) - 1;
            else
                bitmin = -(2^(str2num(header.samplebits)-1));
                bitmax = (2^(str2num(header.samplebits)-1)) - 1;
            end
            multiplier = (calmax-calmin)/(bitmax-bitmin);
            p = (p - bitmin).*multiplier + calmin;
        end
    else
        p = [];    % Output an empty matrix if requested data is beyond the length of the current file
    end
else
    p = [];      % Also output an empty matrix of zero length LENGTH input is requested (ie, only return header/info values)
    info.count = 0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = MakeString(s)
% The standard is that the string is null-terminated,
% so we know there will be a null somewhere.  Make the
% string out of everything up to the null.
%
% The standard also says that the rest of a field should
% be nulled, but in case it wasn't and there is garbage
% after the terminating null, this will ignore the garbage.
%
% Note if the length of s is exactly 1 we just leave it
% alone.  Only if it's longer than one do we look for
% the null terminator.
 s = char(s(:)');
 if (length(s) > 1)		% Multi-char string?
 	xx = find(s==0);	% Yes, look for null terminator
 	if (~isempty(xx))	% Null terminator exists
 		if (xx(1) > 1)	% Is it not the first character?
 			s = strtrim(s(1:(xx(1)-1)));
 		else
 			s = ''; 	% Multi-char string but first character was a null!
 		end
 	else
 		s = '';			% Multi-char string but no null terminator found!
 	end
 else					% Single-char string
 	if (s(1) == 0)		% But is that single char a null?
 		s = '';			% If so, this is a null string
 	end
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function s = strtrim(s)
%  function s = strtrim(s)
%
%  Trims blanks and nulls from both ends of the string.
  if ~isempty(s)
   xx = limits(find((~isspace(s))&(s~=0)));
    if isempty(xx)
     s = '';
    else
     s = s(xx(1):xx(2));
    end
  end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function l=limits(x)
%  Returns a two-element vector containing the lowest and highest values
%  found in all elements of the input array.
 l = [min(x(:)) max(x(:))];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tm = c2mat_tm(tc,fs)
% tm = c2mat_tm(tc)
%  Where tc is a scalar containing a date/time as stored by C, as the
%  number of seconds past Midnight, January 1, 1970, this function will
%  return in tm a six-element row vector containing the same date/time
%  in MATLAB's format.
%
%
% tm = c2mat_tm(tc,fs)
%
%  The additional argument, fs, will be added to the seconds field of the
%  result.  This is to allow times to be specifed to fractions of a second.
%  If no second input argument is given, then the fractional part of tc, if
%  present, will be used.
a = version;
 if (a(1) == 4) && (a(2) == '.')
  error('Missing c2mat_tm.mex');
 else
   if nargin > 2
    tc = tc + fs;
   end
  tm = datevec((tc+6.216730560000000e+010)/(24*60*60));
 end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tc,fs] = mat2c_tm(tm)
% tc = mat2c_tm(tm)
%
%  Where tm is a 6-element row vector containing a date/time in MATLAB's
%  format, this function will return in tc a scalar containing the same
%  date/time in the standard C format, as the number of seconds past
%  Midnight, January 1, 1970.
%
%
% [tc,fs] = mat2c_tm(tm)
%
%  If the second output argument, fs, if given, then only an integer value
%  will be returned in tc, with the fractional part of the seconds in fs.
a = version;
 if (a(1) == 4) && (a(2) == '.')
  error('Missing mat2c_tm.mex');
 else
  tc = datenum(tm(1),tm(2),tm(3),tm(4),tm(5),tm(6))*(24*60*60) - 6.216730560000000e+010;
   if nargout > 1
    fs = tc - fix(tc);
    tc = fix(tc);
   end
 end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function x = contains(what,where,dummy)
% X = contains(WHAT,WHERE)
%
%  For each element of WHAT, this function returns a 1 in the corresponding
%  element of X if that element is found in WHERE, and a 0 if it is not.
%
  if nargin >= 3
   dummy = dummy*dummy;
  end

  if isempty(where)
   x = zeros(size(what));
   return;
  elseif isempty(what)
   x = [];
   return;
  end

 x = logical(zeros(size(what)));
  for ii = 1:numel(what)
   xx = any(what(ii) == where);
    if isempty(xx)
     xx = 0;
    end

    while numel(xx) > 1
     xx = any(xx(:));
    end
   x(ii) = xx;
  end