% 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

daysperyear=365.25;


[xGr,yGr,zGr,xm]=modelsize3(xGr,yGr,zGr); % get xm

[LAYnams,LAYvals,LAYCBD]=getLayers([basename,'.xls'],zGr);

[NROW,NCOL,NLAY]=size(IBOUND);  % get size of model

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);

[MFnams,MFvals]=getExcelData([basename,'.xls'],'MFLOW','Vertical');
ISTRAT=MFvals(index('ISTRAT',MFnams),1);

%% Figure of rotating interface

% setup figure
figure; hold on

leg={};
iLeg=1; % legend index

xlabel('x [m]'); ylabel('z [m]'); grid on;
title(sprintf('SWI example 2, rotation of brackish zone, ISTRAT=%d',ISTRAT));
set(gca,'xlim',xGr([1 end]),'ylim',sort(zGr([1 end])));

% plot all interface positions

clr='brgkcmy'; clr=repmat(clr,1,ceil(NT/length(clr))); % enough colors
stl='-:;.';    stl=repmat(stl,1,ceil(NSURF/length(stl))); % enough styles 

% Color the layers to better see aquifers and aquitards
Layers=NaN(1,NCOL+1,NLAY+sum(LAYCBD)+1);
k=1;
for iLay=1:NLAY
    Layers(:,:,k)=0.9;      % light color if colormap == bone
    k=k+1;
    if LAYCBD(iLay)
        Layers(:,:,k)=0.5;  % dark color if colormap == bone
        k=k+1;
    end
end
caxis([0 1]); colormap('bone');

%hdl=pcolor(xGr,squeeze(zGr),permute(Layers,[3,2,1]));
%set(hdl,'EdgeColor','none');
patch(xGr([1 end end 1 1]),...
      squeeze(zGr([3   3   2 2 3])),'c');

leg{iLeg}='Layers'; iLeg=iLeg+1;

% use the stored interface in ZETA in mf_adapt to plot initial positions
for iSurf=1:length(ZTA(1).label)
    for iLay=1:NLAY
        plot(xm,ZETA(iSurf).values(1,:,iLay),['k' stl(iSurf)]);
        leg{iLeg}=sprintf('IFACE(%d) iLay(%d) T=0',iSurf,iLay);
        iLeg=iLeg+1;
    end
end

% plot computed positions 
for it=1:NT          % for all recoreds 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
            plot(xm,ZTA(it).term{idx}(1,:,iLay),[clr(it) stl(iSurf)]);
            leg{iLeg}=sprintf('IFACE(%d) iLay(%d)T=%.0f y',iSurf,iLay,H(it).totim/daysperyear);
            iLeg=iLeg+1;
        end
    end
end
legend(leg);
% axis equal
% axis tight