function varargout = readcsi(file,channels,Samples,N2Skip)
%function [header,data,datacell] = readcsi(file,channels,Samples,N2Skip)
%
%=========================================================================
% READCSI.M - reads binary or ASCII data logged by Campbell Scientific
% 'table-based' dataloggers (e.g., CR1000, CR3000, CR5000, CR6).
%
%    Copyright (C) Scott Miller 2002-2017.
%    Last Revised: Jan 28, 2017
%
% USAGE: [header,data] = readcsi(file,channels,Samples,N2Skip)
%
% INPUTS
% ------
% file - full file path
% channels - optional input string array that specifies the
%            channels to be read (to extract data); for example
%            channels={'channel1';'channel2'}
% Samples - optional input to read a specific subset of the
%           data.  Input as [firstsample lastsample].
%           If Samples is empty, reads all of the data.
%           If lastsample is NaN, reads from firstsample to end
% N2Skip - optional input to skip N2Skip samples.  Does not
%          do averaging, simply skips samples
%
% OUTPUTS
% -------
% header - character arry of variable names from datalogger program
%
% data [optional] - data array
%
% datacell [optional] - if the data include strings, datacell contains 
%           is the same size as data with strings in cell format 
%           and numeric fields empty
%
% Note: a useful use of this script it to call the function with
% no output variables specified - this will read the header and
% scroll it to the Matlab terminal window, and is useful to query
% what is in a file; i.e.,  >>read3000(file).  Or, specifying only
% one output variable will return just the header as a variable
%
% AUTHOR:  Scott Miller
%          Atmospheric Sciences Research Center
%          SUNY Albany (smiller@albany.edu)
%          http://asrc.albany.edu/people/faculty/miller
%
% DISCLAIMER:
% This software is provided "as is" without warranty of any kind.
%
% Modified 01/27/2017
% - found that some tables are outputing commas within an output string
% that are not separating fields. Replaced my comma search with a regexp
% search that excludes cases when the comma is between double quotes
%
% - all ASCII tables appear forced to use a string format for the timestamp
% like 2017-01-08 01:00:00. Assuming that is true, modifying this code to
% unpack the timestamp during the read process, rather than do it in later
% steps. So, if the timestamp was the only ascii field and the rest are all
% numeric, do not need to output the cell array
%
% Modified 12/31/2016
% - added capability to read Campbell FP2 format
% - for ASCII file format (TOA), all data are no longer put into cells;
% instead, if there are mixed strings and numeric output, there are 2 data
% arrays of identical size returned by the script: 'data' contains numeric
% data with string fields empty, and 'datacell' contains strings in cell
% format with numeric fields empty
%
% Modified 7/7/2013
% - changed the search string in the header that
% separates fields from a comma to "," since commas are used in the
% fieldnames in some of the fields of the status table
%

if nargin < 4
    N2Skip=0;
end

if nargin < 3
    Samples = [];
end

if nargin < 2
    channels=[];
end

CR=13;
LF=10;
COMMA=44;

% added deblank Aug 9, 2006
fid=fopen(deblank(file),'r','ieee-le');
d=fread(fid,inf,'uchar');
%d=fread(fid,10000,'uchar');

% find Line Feeds and Carriage Returns

icr = find(d==CR);
ilf = find(d==LF);

%d(icr)

%size(icr)
%size(ilf)

% disp(['line 2: ' int2str(d([icr(2) icr(2)+1 icr(2)+2])')]);
% disp(['line 3: ' int2str(d([icr(3) icr(3)+1 icr(3)+2])')]);
% disp(['line 4: ' int2str(d([icr(4) icr(4)+1 icr(4)+2])')]);
% disp(['line 5: ' int2str(d([icr(5) icr(5)+1 icr(5)+2])')]);
% disp(['line 6: ' int2str(d([icr(6) icr(6)+1 icr(6)+2])')]);
% disp(['line 7: ' int2str(d([icr(7) icr(7)+1 icr(7)+2])')]);

%char(d(4))

switch d(4)
    case abs('A')
        
        EOH=ilf(4);
        
        % read last line of header to get the file structure
        HLine2 = d(ilf(1)+1:ilf(2)-1)';
        HLine3 = d(ilf(2)+1:ilf(3)-1)';
        HLine4 = d(ilf(3)+1:ilf(4))';
        
        % disp(char(HLine2));
        % disp(char(HLine3));
        % disp(char(HLine4));
        
        % Jul 7 2013, rework this search to account for commas used in the field
        % name. This is done in the datalogger status files
        
        begfields2 = [1 strfind(HLine2,'","')+2];
        endfields2 = [strfind(HLine2,'","') length(HLine2)-1];
        
        %begfields2 = [1 find(HLine2==COMMA)+1];
        %endfields2 = [find(HLine2==COMMA)-1 length(HLine2)-1];
        
        begfields3 = [1 find(HLine3==COMMA)+1];
        endfields3 = [find(HLine3==COMMA)-1 length(HLine3)-1];
        
        begfields4 = [1 find(HLine4==COMMA)+1];
        endfields4 = [find(HLine4==COMMA)-1 length(HLine4)-1];
        
        Nfields = length(begfields2);
        
        for i=1:Nfields
            % don't read the quotes at beginning and end of each field
            FieldName{i} = char(HLine2(begfields2(i)+1:endfields2(i)-1));
            FieldUnits{i}= char(HLine3(begfields3(i)+1: endfields3(i)-1));
            disp([int2str(i) ' of ' int2str(Nfields) ', name: ' FieldName{i}]);
            disp([int2str(i) ' of ' int2str(Nfields) ', units: ' FieldUnits{i}]);
        end
        
    case abs('B')
        % nov 10, 2006: redoing this because it chokes when the first
        % character of the binary data after the last line of the header
        % turns out to be a double quote in ascii - this happens some
        % times!
        %
        % find the end of the binary header - find the last occurrence of
        % the ascii combination [34 13 10] which is
        % ['double quote' 'carriage return' 'line feed' ]
        % this line is the last line of the header
        
        %  strfind(char(d'),char([34 13 10]))
        nrowsa=length(strfind(char(d'),char([34 13 10])));
        
        % double check the end of the binary header by looking for the ascii combination
        % [34 13 10 34] which is
        % ['double quote' 'carriage return' 'line feed' 'double quote' ]
        % this line is the second last line of the header
        
        %  strfind(char(d'),char([34 13 10 34]))
        
        %disp('11-24-06: using the search for the second-last header row, instead of last row');
        nrowsb=length(strfind(char(d'),char([34 13 10 34])))+1;
        nrows=min([nrowsa nrowsb]);
        
        EOH=ilf(nrows);
        
        % $$$   % find the end of the binary header - find the last occurrence of
        % $$$   % the ascii combination [13 10 34] which is
        % $$$   % ['carriage return' 'line feed' 'double quote']
        % $$$   % The line following this line is the last line of the header
        % $$$
        % $$$   nrows=length(strfind(char(d'),char([13 10 34])))+1
        % $$$
        % $$$   EOH=ilf(nrows)
        
        % read last line of header to get the file structure
        HLineName = d(ilf(nrows-4)+1:ilf(nrows-3)-1)';
        HLineUnits = d(ilf(nrows-3)+1:ilf(nrows-2)-1)';
        HLinePrecision = d(ilf(nrows-1)+1:ilf(nrows))';
        
        % $$$   disp(char(HLineName));
        % $$$   disp(char(HLineUnits));
        % $$$   disp(char(HLinePrecision));
        
        begfieldsName = [1 find(HLineName==COMMA)+1];
        endfieldsName = [find(HLineName==COMMA)-1 length(HLineName)-1];
        
        begfieldsUnits = [1 find(HLineUnits==COMMA)+1];
        endfieldsUnits = [find(HLineUnits==COMMA)-1 length(HLineUnits)-1];
        
        begfieldsPrecision = [1 find(HLinePrecision==COMMA)+1];
        endfieldsPrecision = [find(HLinePrecision==COMMA)-1 length(HLinePrecision)-2];
        
        Nfields = length(begfieldsName);
        
        for i=1:Nfields
            % don't read the quotes at beginning and end of each field
            FieldName{i}=char(HLineName(begfieldsName(i)+1:endfieldsName(i)-1));
            FieldUnits{i}=char(HLineUnits(begfieldsUnits(i)+1: endfieldsUnits(i)-1));
            FieldPrecision{i}=char(HLinePrecision(begfieldsPrecision(i)+1: endfieldsPrecision(i)-1));
        end
end

%FieldPrecision

%EOHA = ilf(4);
%EOHB = ilf(5);

if nargout<2
    header=[];
    data=[];
    for i=1:Nfields
        switch d(4)
            case abs('A')
                disp([char(FieldName{i}) ', (' char(FieldUnits{i}) ')']);
                header=strvcat(header,char(FieldName{i}));
            case abs('B')
                disp([char(FieldName{i}) ', (' char(FieldUnits{i}) '), ('  char(FieldPrecision{i}) ')']);
                header=strvcat(header,char(FieldName{i}));
        end
    end
    return
end

% determine if file is ASCII or binary
switch d(4)
    case abs('A')
        header=[];
        units=[];
        for ifield=1:Nfields
            header=strvcat(header,char(FieldName{ifield}));
            units=strvcat(units,char(FieldUnits{ifield}));
        end
        
        disp(['reading ASCII file ' file ' ...']);
        fseek(fid,0,'bof');
        ftell(fid);
        
        d=fread(fid,inf,'uchar');
        d=d(EOH+1:end); % remove header
        
        % find Line Feeds and Carriage Returns
        icr = find(d==CR);
        ilf = [0; find(d==LF)];
        
        nlines=length(ilf);
        
        data=nan(nlines,Nfields);
        %datacell=cell(nlines,Nfields);
        
        % after the header, read each line of data
        for iline=1:nlines-1
            
            TextLine = d(ilf(iline)+1:ilf(iline+1)-1)';
            
            %disp([int2str(iline) ': ' char(TextLine)]);
            
            pattern = ',(?=(?:[^\"]*\"[^\"]*\")*[^\"]*$)';
            begfields=[1 regexpi(char(TextLine),pattern)+1];
            endfields=[regexpi(char(TextLine),pattern)-1 length(TextLine)-1];
                
            %begfields = [1 find(TextLine==COMMA)+1];
            %endfields = [find(TextLine==COMMA)-1 length(TextLine)-1];
                        
            %Nfields
            NfieldsInLine = length(begfields);
            %pause
            
            if NfieldsInLine==Nfields
                for ifield=1:Nfields % for each field in line
                                    
                    if strcmp(upper(char(FieldName{ifield})),'TIMESTAMP')
                        
                        %disp([int2str(ifield) ' is timestamp']);
                        
                        %generate timestamp from first element of each record
                        %1234567890123456789
                        %2017-01-08 01:00:00
                        
                        tstmp=char(TextLine(begfields(ifield)+1: endfields(ifield)-1));
                        
                        if ~isempty(tstmp)
                            yyyy=str2num(tstmp(:,1:4));
                            mm=str2num(tstmp(:,6:7));
                            dd=str2num(tstmp(:,9:10));
                            
                            hh=str2num(tstmp(:,12:13));
                            minutes=str2num(tstmp(:,15:16));
                            seconds=str2num(tstmp(:,18:19));
                            
                            data(iline,ifield)=datenum(yyyy,mm,dd,hh,minutes,seconds)-datenum(yyyy,1,0);
                        end
                        
                    elseif isempty(findstr(char(TextLine(begfields(ifield): endfields(ifield))),'"')) % check for quotes - quotes read as string
                        
                        % disp('no quotes')
                        ThisValue=char(TextLine(begfields(ifield): endfields(ifield)));
                        %disp(char(TextLine(begfields(ifield): endfields(ifield))));
                        
                        % if it is a purely numberic value, put into data table
                        if all(ismember(ThisValue, '0123456789+-.eEdD'))
                            data(iline,ifield) = str2num(ThisValue); %char(TextLine(begfields(ifield): endfields(ifield)));
                        else
                            % if it is mixed with non numeric characters, put into a cell
                            % arra
                            datacell{iline,ifield} = char(TextLine(begfields(ifield): endfields(ifield)));
                        end
                        
                    else
                        % don't read the quotes at beginning and end of each field
                        datacell{iline,ifield} = char(TextLine(begfields(ifield)+1: endfields(ifield)-1));
                    end
                    %         disp(data(iline,ifield));
                end
            else
                %       disp(['Line ' int2str(iline) ' has ' int2str(NfieldsInLine) ' fields: needs to have ' int2str(Nfields)]);
            end
            %      pause
        end
        data=data';
        if exist('datacell');
            datacell=datacell';
        end
        
    case abs('B')
        % Calculate the number of bytes in a record and get the
        % corresponding matlab precision
        %  FieldPrecision
        for i=1:size(FieldPrecision,2)
            if strcmp(char(FieldPrecision(i)),'ULONG')
                NBytes(i) = 4;
                MatlabPrec{i}='uint32';
            elseif strcmp(char(FieldPrecision(i)),'LONG')
                NBytes(i) = 4;
                MatlabPrec{i}='uint32';
            elseif strcmp(char(FieldPrecision(i)),'FP2')
                NBytes(i) = 2;
                MatlabPrec{i}='uint16';
            elseif  strcmp(char(FieldPrecision(i)),'IEEE4L')
                NBytes(i) = 4;
                MatlabPrec{i}='float32';
            elseif  strcmp(char(FieldPrecision(i)),'IEEE4')
                NBytes(i) = 4;
                MatlabPrec{i}='float32';
            elseif  strcmp(char(FieldPrecision(i)),'IEEE8') % guessing at this - it is in loggernet manual but not cr6 manual
                NBytes(i) = 8;
                MatlabPrec{i}='float64';
            elseif  strcmp(char(FieldPrecision(i)),'SecNano')
                NBytes(i) = 8;
                % Jan 12, 2013: this might need to be changed.  Apparently SecNano is
                % 2 longs, the first is seconds and the second is Nanoseconds
                MatlabPrec{i}='uint64';
            elseif  strfind(char(FieldPrecision(i)),'ASCII')
                tmp=char(FieldPrecision(i));
                NBytes(i) = str2num(tmp(7:strfind(tmp,')')-1));
                % Jan 12, 2013: this might need to be changed.  Apparently SecNano is
                % 2 longs, the first is seconds and the second is Nanoseconds
                MatlabPrec{i}=tmp;
            end
        end
        
        % MatlabPrec
        
        %%% Start reading the channels
        header = [];
        units = [];
        
        % first position pointer at the end of the header
        fseek(fid,EOH,'bof');ftell(fid);
        
        BytesPerRecord=sum(NBytes)*ones(size(NBytes)) ;
        BytesCumulative = [0 cumsum(NBytes(1:length(NBytes)-1))];
        
        % check if a subset of the data is to be read, or the entire file
        % should be read.  Or, if Samples is nonzero, but the second
        % element of Samples is NaN, then read to the end of the file, by
        % specifiying N2Read as inf
        if isempty(Samples)
            N2Read=inf;
            Samples=1;
        elseif isnan(Samples(2))
            N2Read = inf;
        else
            N2Read = (Samples(2)-Samples(1)+1)/(N2Skip+1);
        end
        
        if isempty(channels)
            channels = char(FieldName);
        end
        
        fprintf(1,'Reading %4.0f channels from [ %s ]:',size(channels,1),file);
        
        if size(channels,1) >= 10
            fprintf(1,'\n');
        end
        
        NColumns=30+length(file);
        
        for j=1:size(channels,1)

            ich=strmatch(deblank(channels(j,:)),FieldName,'exact');
            
            Beg = EOH + BytesCumulative(ich) + BytesPerRecord(1)*(Samples(1)-1);
            
            fseek(fid,Beg,'bof');
            
            %disp([int2str(j) ':  ' char(MatlabPrec(ich))]);
            
            % unpack FP2 format if necessary
            if strcmp(char(FieldPrecision(j)),'FP2')
                % 12/31/16 - note: explicitly specifying 'b' for big endian (CR6 datalogger).
                % If the logger is CR3000/CR1000, this might need to be modified to
                % 'l'. In that case, would need to return the datalogger type with
                % the header
                data(:,j)= rd2byte(fread(fid,N2Read,char(MatlabPrec(ich)),(N2Skip+1)*BytesPerRecord-NBytes(ich),'b')) ;
            % deal with ascii
            elseif strfind(char(FieldPrecision(j)),'ASCII')
                data(:,j)= NaN ;
            
            else 
                %size(data)
                %j
                %N2Read
                %char(MatlabPrec(ich))
                %(N2Skip+1)*BytesPerRecord-NBytes(ich)
                data(:,j)= fread(fid,N2Read,char(MatlabPrec(ich)),(N2Skip+1)*BytesPerRecord-NBytes(ich)) ;
            end
            
            %data(:,j)
            
            header=strvcat(header,char(FieldName{ich}));
            units=strvcat(units,char(FieldUnits{ich}));
            
            fprintf(1,' .');
            
            NColumns=NColumns-2;
            if NColumns <= 0
                fprintf(1,'\n');
                NColumns=30+length(file);
            end
        end
        
        fprintf(1,'\n');
        data=data';
        %datacell=[];
end

varargout{1}=header;
varargout{2}=data;
varargout{3}=units;

if exist('datacell')
    varargout{4}=datacell;
end

fclose(fid);

end

function d2byte = rd2byte(words)

%6/1/2000

FirstBit  = 32768; %bin2dec('1000000000000000')
Bits2And3 = 24576; %bin2dec('0110000000000000')
Bits2Not3 = 16384; %bin2dec('0100000000000000')
Bits3Not2 = 8192;  %bin2dec('0010000000000000')
Bits4to16 = 8191;  % bin2dec('0001111111111111')

d2byte=NaN*ones(size(words));

Mult=ones(size(words));

Mult(find(bitand(words(:),FirstBit)==FirstBit)) = -1*...
    Mult(find(bitand(words(:),FirstBit)==FirstBit));

Mult(find(bitand(words,Bits2And3) == Bits3Not2 )) = 0.1 * ...
    Mult(find(bitand(words,Bits2And3) == Bits3Not2 )) ;

Mult(find(bitand(words,Bits2And3) == Bits2Not3 )) = 0.01 * ...
    Mult(find(bitand(words,Bits2And3) == Bits2Not3 )) ;

Mult(find(bitand(words,Bits2And3) == Bits2And3 )) = 0.001 * ...
    Mult(find(bitand(words,Bits2And3) == Bits2And3 ));

d2byte = Mult.*bitand(words,Bits4to16);

% disp([num2str(toc,4) ' sec vectorized']);

end


