function bud=readBud(fname,userLabels,userPeriods,userTstps,userLayers,userRows,userCols)
% bud=readBud([-]fname [,userLabels [,userPeriods [,tsteps [,userLayers [,userRows [,userCols]]]]]])
% ---- read binary MODFLOW output into an Data(nP,nT) array of structs
% you may subselect any part of the file by qualification using the optional
% arguments. These arguments may be entirely omitted or empty as shown
%
% put a - in front of fname to reduce screen output
%
% use 
%    readBud(fname,-1);
% to get infor about the contents of the file
%
% here as examples of usage:
%
%    readbud(fname,'',[1,3,6],'',3)
%
% this extracts all userLabels, only stress userPeriods 1,3 and 6 and only layer 3.
%
%    readbud(fname,{'STORAGE','FLOWRIGHTFACE'},[1 3 5],'',5,3:10,12:21),
%
% this extracts only the 'STORAGE' and 'FLOWRIGHTFACE' flow terms and these only
% for stress userPeriods 1,3 and 6, all time steps and layer 5 and userRows 3
% through 10 and columns 12 through 31.
% If only 1 flow term is desired you may omit the braces { }.
% The layers, userRows and columns can be extracted in any order by adequat
% indexing as is standard in Matlab.
%
% Readbud, therefore, lets you extract any portion of the budget file, no
% matter how big it is, and fast.
%
% The resulting variable, here named bud, is an struct array. Its length
% equals the number of records in the file (or leff if you limited your
% selection by qualifying userPeriods and or timesteps).
% This struct has the following fields
%   .label        cell array of the label names in the file, which refer to
%                 the cell by cell flow terms it contains
%   .term(NLabel) cell array with same length and order as label in which
%                 each cell has a full 3D array cell by cell flow terms
%                 corresponding to the label in the same position
%   .pertim       time within this period
%   .totim        total time from start of simulation
%   .period       current stress period number
%   .tstp         current time step number in this stress period
%   .NROW         total number of userRows    in budget file
%   .NCOL         total number of columns in budget file
%   .NLAY         total number of layers  in budget file
%   .userCols         list of the extracted columns
%   .userRows         list of the extracted userRows
%   .userLayers         list of the extracted layers
%   
% for instance to get to a row and column in a layer of a given flow term
% which corresponds to the second label:
% B(iRecord).term{iLabel}(yourrows,yourcols,yourlays)
% which is the standard Matlab way of indexing.
%
% TO 090104 091214
% TO 100602 make budlabels unique because labels may vary between stress userPeriods
% TO 100818 logic redone to clean up in accordance with readDat and readMT3D
%
% Copyright 2009 Theo Olsthoorn, TU-Delft and Waternet, without any warranty
% under free software foundation GNU license version 3 or later

if fname(1)=='-', fname=fname(2:end); verbose=0; else verbose=1; end

fp=fopen(fname); if fp<1, error('READBUD: Can''t find or open file %s!',fname); end

if exist('userLabels','var')
    if ischar(userLabels), userLabels={userLabels}; end  % must be cell of strings
    userLabels=unique(upper(userLabels));
end

%% Advance a single cell-by-cell record & compute the number or records in file

fprintf('Trying to read %s as BINARY file...',fname); 

bytes=contrRec(fp);

nbyt=ftell(fp); fseek(fp,0,'eof'); Nbyt=ftell(fp); nRec=Nbyt/nbyt; % number of records in file
if rem(Nbyt,nbyt),
   fprintf('failed\n');
   error('Budget file %s unknown binary format! Can''t help it, use a better compiler\n',fname);
end
fprintf('it works!\n');

%% Get the information about the contents of the file from record headings

periodsInFile = zeros(nRec,1);  % here goes the stress period number of each record
tstpsInFile   = zeros(nRec,1);  % and here the time step number
iLabInFile    = zeros(nRec,1);  % here the label numer (we'll see how this works)
labelsInFile  = cell(nRec,1);   % here we store the actual label of the records
layersInFile  = zeros(nRec,1);  % each record in the budget files holds the NLAY
rowsInFile    = zeros(nRec,1);  % same for NROW
colsInFile    = zeros(nRec,1);  % same for NOL  we store them all, no exceptions

%% Store record heading data, but do not yet the values,
% because we may perhaps only need a small selection

fprintf('Scanning headers\n');
for i=1:nRec
    fseek(fp,(i-1)*nbyt,'bof');  % move to start of last layer in file
    [periodsInFile(i),...
        tstpsInFile(i),...
        labelsInFile{i},...
        layersInFile(i),...
        rowsInFile(i),...
        colsInFile(i)]=contrRec(fp,bytes);
    if rem(i, 100)==0, fprintf('.'); end
    if rem(i,2500)==0, fprintf('%d records read\n',i); end
end
fprintf('finished, %d records read\n',i);
    
%% use indices into unique list of labels instead of labels themselves
% iLabInFile is index into the list of unique labels in the budget file
budlabels=unique(labelsInFile);
for i=1:nRec, iLabInFile(i)=strmatchi(labelsInFile{i},budlabels); end

%% Overview of what's in the file
NPER = length(unique(periodsInFile));
NSTP = length(unique(tstpsInFile));
NLBL = length(budlabels);     % already unique
NCOL = max(colsInFile);
NROW = max(rowsInFile);
NLAY = max(layersInFile);

pif=unique(periodsInFile);
tif=unique(tstpsInFile);

fprintf('File contains the following:\n');
fprintf('Stress Periods           : '); fprintf('% 3d',pif); fprintf('\n');
fprintf('Time steps               : '); fprintf('% 3d',tif); fprintf('\n');
fprintf('Number of records in file: %10d\n',nRec);
fprintf('Number of stress periods : %10d\n',length(pif));
fprintf('Number of time steps     : %10d\n',length(tif));
fprintf('Number of layers         : %10d\n',NLAY);
fprintf('Number of Rows           : %10d\n',NROW); 
fprintf('Number of columns        : %10d\n',NCOL);
fprintf('Number of unique labels  : %10d\n',NLBL);
for i=1:NLBL, fprintf('%s\n',budlabels{i}); end

if exist('userLabels','var') && isnumeric(userLabels) && userLabels(1)<=0, return; end

%% select user specified selecton or all

% translate userLabels into indices into list of unique labels
if exist('userLabels','var') && isempty(strmatch('',userLabels,'exact')) % verify if user specified non-empty labels
    % transl
    iUserLabels=zeros(size(userLabels));
    for iL=1:length(userLabels)
        iUserLabels(iL)=strmatchi(userLabels{iL},budlabels); % automatically generates an error if label illegal
    end
else % user want all
    userLabels=budlabels;              % all userLabels in the file in order
    iUserLabels=1:length(budlabels);   % their indices for storage
end


%% We create the same sort of selectors for the stress userPeriods
if  exist('userPeriods','var') && ~isempty(userPeriods)
    if ~isempty(userPeriods(userPeriods<1 | userPeriods>NPER))
        error('Requested userPeriods beyond range(1..%d) in budget file %s!',NPER,fname);
    end
    userPeriods=unique(userPeriods);
else
    userPeriods=unique(periodsInFile);
end

%% and also for the time steps
if exist('userTstps'  ,'var') && ~isempty(userTstps)
    if ~isempty(userTstps(userTstps<1 | userTstps>NSTP))
        error('Requested userPeriods beyond range(1..%d) in buget file %s!',NSTP,fname);
    end
    userTstps=unique(userTstps);
else
    userTstps=unique(tstpsInFile);
end


%% Selectors have values for all records of the inptut file. There
%  contents are userPeriod index, the userTstep index, the userLayers index
%  the combined userPeriod and userTstep index plus
%  the output record index and the output layer index
ULABC=1; UPCOL=2; UTCOL=3; UPTCOL=4;

Select=zeros(nRec,UPTCOL);  
for i=1:length(userLabels),  Select(iLabInFile    == iUserLabels(i),ULABC)=i; end
for i=1:length(userPeriods), Select(periodsInFile == userPeriods(i),UPCOL)=i; end
for i=1:length(userTstps),   Select(tstpsInFile   == userTstps(i)  ,UTCOL)=i; end

%% records in inFile to deal with (LPT=Label Period TimeStep Layer combinations)

%  = inFile records to deal with
LPTL=find(Select(:,ULABC)>0 & Select(:,UPCOL)>0 & Select(:,UTCOL)>0);

% unique period tstp combinations
userPT = unique([periodsInFile(LPTL) tstpsInFile(LPTL)],'rows');

% output records in Select(:,UPTCOL)
for i=1:size(userPT,1),
    Select(periodsInFile ==userPT(i,1) & tstpsInFile==userPT(i,2),UPTCOL)=i;
end

%% Rows and columns
if exist('userRows'   ,'var') && ~isempty(userRows)
    if ~isempty(userRows(userRows<1 | userRows>NROW))
        error('Requested userRows are beyond the range(1..%d) in budget file %s!',NROW,fname);
    end
else
    userRows=1:NROW;
end

if exist('userCols'   ,'var') && ~isempty(userCols)
    if ~isempty(userCols(userCols>NCOL | userCols<1))
        error('Requested columns are beyond the range(1..%d) in budget file %s!',NCOL,fname);
    end
else
    userCols=1:NCOL;
end

if exist('userLayers'   ,'var') && ~isempty(userLayers)
    if ~isempty(userLayers(userLayers<1 | userLayers>NLAY))
        error('Requested layers are beyond range(1..%d) in budget file %s!',NLAY,fname);
    end
else
    userLayers=1:NLAY; % all layers
end

%% Get the unique period and tstep combinations for output records

nRecOut = length(unique(Select(LPTL,UPTCOL)));

for i=1:length(userLabels)
    % generate a dummy matrix for preallocation. Notice that the matrix may
    % be substantially smaller than the full 3D grid array because of user
    % selections of userPeriods, time steps, layers, userRows and columns !
    B.term{i,1}=NaN(length(userLayers),length(userRows),length(userCols));
    B.label{i,1}='';
end

bud=repmat(B,[nRecOut,1]); % size of output sttruct array is allocated here

%% Get the actual data values, here we go ...
for i=1:length(LPTL)
    iRecIn=LPTL(i);
    fseek(fp,nbyt*(iRecIn-1),'bof');
    iROut =Select(iRecIn,UPTCOL);   % output record Nr
    iLabel=Select(iRecIn,ULABC);

    [bud(iROut).period,...
        bud(iROut).tstp,...
        bud(iROut).label{iLabel},...
        bud(iROut).NLAY,...
        bud(iROut).NROW,...
        bud(iROut).NCOL,...
        values]=contrRec(fp,bytes);

    bud(iROut).term{iLabel}=values(userRows,userCols,userLayers); % values is a 3D array in the Budget file
    bud(iROut).lays=userLayers;
    bud(iROut).rows=userRows;
    bud(iROut).cols=userCols;

    if verbose
        fprintf('Label(per=%3d,tstp=%3d) =%s\n',bud(iROut).period,bud(iROut).tstp,bud(iROut).label{iLabel});
    else
        fprintf('.');
        if rem(iRecIn,50)==0;
            fprintf('%d\n',iRecIn);
        end
    end
end
if ~verbose && rem(iRecIn,50)~=0, fprintf('\n'); end

% That's all ... TO 090105

fclose(fp);
end

function [kper,kstp,label,nlay,nrow,ncol,values]=contrRec(fp,bytes)
% [kstp,kper,text,nlay,ncol,nrow,data]=contrRec(fid,bytes)
% --- reads a complete layer from a UNFORMATTED MODFLOW budget file
% --- if nargin==1, only the number of bytes at the beginning and enf
% of every record is returned. This is compiler dependent.
% once we know the number of bytes, a full record is returend.
% TO 070703 091214

if nargin==1, % just get the offset of this unformatted file
    fread(fp, 2,'int');  % kper, kstp
    label =char(fread(fp,16,'char')');
    
    % find offset from last null byte in label
    bytes=find(label==0,1,'last'); if isempty(bytes), bytes=0; end
    fread(fp,bytes,'int8'); % offset
    
    ncol=fread(fp, 1,'int');    % ncol nrow nlay
    nrow=fread(fp, 1,'int');    % ncol nrow nlay
    nlay=fread(fp, 1,'int');    % ncol nrow nlay
    
    fread(fp,bytes,'int8'); % offset
    fread(fp,bytes,'int8'); % offset
    
    n=ftell(fp); fread(fp,1,'float'); floatlen=ftell(fp)-n;
    fseek(fp,floatlen*(ncol*nrow*nlay-1),0);
    
    fread(fp,bytes,'int8'); % offset
    
    kper=bytes;
    return;
end

% offset is supposed to be give in bytes if narout>1
fread(fp,bytes,'int8');

kstp  =fread(fp, 1,'int');
kper  =fread(fp, 1,'int');
label =char(fread(fp,16,'char')');
label(label==' ')='';

ncol  =fread(fp, 1,'int');
nrow  =fread(fp, 1,'int');
nlay  =fread(fp, 1,'int');

fread(fp,bytes,'int8');
fread(fp,bytes,'int8');

if nargout<7, % don't fetch but skip
    n=ftell(fp); fread(fp,1,'float'); floatlen=ftell(fp)-n;
    fseek(fp,floatlen*(ncol*nrow*nlay-1),0);
else % fetch
    values=permute(reshape(fread(fp,ncol*nrow*nlay,'float'),[ncol,nrow,nlay]),[2,1,3]);  % is now nrow,ncol,nlay  (looking from above, Matlab style)    
end

fread(fp,bytes,'int8');

end
