getwdsdata(wdsfile, code, id)

function [x, t] = getwdsdata(wdsfile, code, id)
% getwdsdata
%
% Get water distribution system data from an EPANET input file
% 
% Syntax
%
%   [x] = getwdsdata(wdsfile, code)
%   [x] = getwdsdata(wdsfile, code, id)
%   [x,t] = getwdsdata(...)
%
% Description
%
% [x] = getwdsdata(wdsfile, code) returns static water distribution system
% data for all nodes/links from an EPANET 2 input file. The type of data is
% specified by the character array 'code'.
%
% [x] = getwdsdata(wdsfile, code, id) returns static water distribution
% system data for the nodes/links specified in cell array of strings 'id'
% from an EPANET 2 input file.
%
% [x,t] = getwdsdata(...) returns dynamic water distribution
% system data with the time steps from an EPANET 2 input file.
%
% Example
%
% [cl, t] = getwdsdata('Net1.inp', 'en_quality', {'2', '12'})
%
% Author: Philip Jonkergouw
% Email:  [e-mail]
% Date:   July 2007


% Initialise a few variables ...
code = upper(code);
errorcode = 0;
x = []; t = [];
if (nargin < 3) id = []; end

% Find out what is required ...

% Static node value?
s = 'EN_ELEVATIONEN_BASEDEMANDEN_PATTERNEN_EMITTEREN_INITQUALEN_SOURCEQUALEN_SOURCEPATEN_SOURCETYPEEN_TANKLEVEL';
if ~isempty(findstr(s, code))
    isstatic = 1;
    location = 'node';
end

% Dynamic node value?
s = 'EN_DEMANDEN_HEADEN_PRESSUREEN_QUALITYEN_SOURCEMASS';
if ~isempty(findstr(s, code))
    isstatic = 0;
    location = 'node';
end

% Static link value?
s = 'EN_DIAMETEREN_LENGTHEN_ROUGHNESSEN_MINORLOSSEN_INITSTATUSEN_INITSETTINGEN_KBULKEN_KWALL';
if ~isempty(findstr(s, code))
    isstatic = 1;
    location = 'link';
end

% Dynamic link value?
s = 'EN_FLOWEN_VELOCITYEN_HEADLOSSEN_STATUSEN_SETTINGEN_ENERGY';
if ~isempty(findstr(s, code))
    isstatic = 0;
    location = 'link';
end

% Water quality data?
ishydraulic = 1;
if strcmp(code, 'EN_QUALITY') ishydraulic = 0; end

% Load the EPANET 2 dynamic link library ...
if ~libisloaded('epanet2') loadlibrary('epanet2', 'epanet2.h'); end

% Open the water distribution system ...
s = which(wdsfile);
if ~isempty(s) wdsfile = s; end

[errorcode] = calllib('epanet2', 'ENopen', wdsfile, '', '');
if (errorcode)
    fprintf('Could not open network ''%s''.\nReturned empty array.\n', wdsfile);
    return;
end

% Retrieve the indices ...
nlocations = length(id);
if ~nlocations
    % Retrieve data for all nodes/links
    countcode = getenconstant(['EN_', upper(location), 'COUNT']);
    [errorcode, nlocations] = calllib('epanet2', 'ENgetcount', countcode, nlocations);
    if (errorcode)
        fprintf(['Problem retrieving number of ', location, 's.\n']);
        [errorcode] = calllib('epanet2', 'ENclose');
        unloadlibrary('epanet2');
        return;
    end
    nlocations = double(nlocations);
    index = 1:nlocations;
else
    getindexfunc = ['ENget', location, 'index'];
    for n = 1:nlocations
        index(1,n) = 0;
        [errorcode, id{n}, index(1,n)] = calllib('epanet2', getindexfunc, id{n}, index(1,n));
        if (errorcode)
            fprintf(['Problem retrieving index for ', location, ' ''', id{n}, '''.\n']);
            [errorcode] = calllib('epanet2', 'ENclose');
            unloadlibrary('epanet2');
            return;
        end
    end
end

getvaluefunc = ['ENget', location, 'value'];
% Obtain correct epanet code ...
epanetcode = getenconstant(code);

if isstatic
    % Static output required. No need to run simulation ...
    for n = 1:nlocations
        x(n) = 0;
        [errorcode, x(n)] = calllib('epanet2', getvaluefunc, index(n), epanetcode, x(n));
        if (errorcode)
            s = '                   ';
            getidfunc = ['ENget', location, 'id'];
            [errorcode, s] = calllib('epanet2', getidfunc, index(n), s);
            if (errorcode)
                fprintf(['Problem retrieving value for ', location, ' ''', x(n), ''' (index).\n']);
            else
                fprintf(['Problem retrieving value for ', location, ' ''', s, '''(id).\n']);
            end
            [errorcode] = calllib('epanet2', 'ENclose');
            unloadlibrary('epanet2');
            return;
        end
    end
else
    % Non-static output required. Run simulation ...
    nextfunc = 'ENnextH';
    initflag = 11;
    enflag = 'H';
    if ~ishydraulic
        enflag = 'Q';
        nextfunc = 'ENstepQ';
        % Solve hydraulics prior to running WQ simulation ...
        [errorcode] = calllib('epanet2', 'ENsolveH');
        initflag = 0;
    end
    
    % Assign function calls (WQ or Hydraulic) ...
    initfunc  = ['ENinit', enflag];
    openfunc  = ['ENopen', enflag];
    runfunc   = ['ENrun', enflag];
    closefunc = ['ENclose', enflag];
    
    % Open and initialise the (WQ/Hydraulic) solver ...
    [errorcode] = calllib('epanet2', openfunc);
    [errorcode] = calllib('epanet2', initfunc, initflag);
    
    % Initialise some variables ...
    tval = 0; tstep = 1; value = 0.0; pt = 0;
    
    % Loop through simulation ...
    while tstep && ~errorcode
        [errorcode, tval] = calllib('epanet2', runfunc, tval);
        pt = pt + 1;
        for n = 1:nlocations
            [errorcode, value] = calllib('epanet2', getvaluefunc, index(n), epanetcode, value);
            x(pt, n) = double(value);
        end
        t(pt) = double(tval)/3600;
        % Continue to the next time step ...
        [errorcode, tstep] = calllib('epanet2', nextfunc, tstep);
    end
    % Close the solver
    [errorcode] = calllib('epanet2', closefunc);
end
% Close EPANET ...
[errorcode] = calllib('epanet2', 'ENclose');
if (errorcode) fprintf('EPANET error occurred. Code %g\n', num2str(errorcode)); end
% Convert to double precision ...
x = double(x);

if libisloaded('epanet2') unloadlibrary('epanet2'); end