% Showing the results (positon of the interface) of example 1
% TO 091224

load name            % get basename of this model
load(basename);      % load data for this model (generaged in mf_adapt

stitle=sprintf('%s: Interf. at T=',basename);  % initial title text

daysperyear=365.25;

%% Grid  xGr,yGr and zGr, ZETA and IBOUND are loaded

[NROW,NCOL,NLAY]=size(IBOUND);

H  =readDat([basename,'.HDS']);  % use only H.totim
ZTA=readBud([basename,'.ZTA']);  % get interface positions

NT   =length(ZTA);      % total number of time steps in ZTA file
NSURF=length(ZTA(1).label);

%% Figure of rotating interface

figure; hold on

xlabel('x [m]'); ylabel('y [m]'); zlabel('z [m]'); grid on;

set(gca,'xlim',xGr([1 end]),'ylim',sort(yGr([1 end])),'zlim',sort(zGr([1 end])));

%% Initial interface positions using ZETA
% gridsurf is used to match xGr yGr with values between them at xm and ym

for iSurf=1:length(ZETA)
    for iLay=1:NLAY  % we only have one layer in this case
        gridsurf(xGr,yGr,ZETA(iSurf).values(:,:,iLay)); hold on
    end
end

%% Computed positions using ZTA
for it=[1,NT]          % for all records in bud file (and H-file)
    for iSurf=1:length(ZTA(1).label)
        idx=strmatch(sprintf('ZETAPLANE%d',iSurf'),ZTA(it).label);
        for iLay=1:NLAY  % we only have one layer in this case
            gridsurf(xGr,yGr,ZTA(it).term{idx}(:,:,iLay));
            stitle=[stitle, sprintf(' %d',round(H(it).totim/daysperyear))];
        end
        stitle=[stitle,' y'];
    end
end
title(stitle);

view(10,5);
