%% Analyzing model output
load('name'); load(basename);

%% The function XS turns the 3D vertical matrix into a 2D flat matrix for
% easy contouring and convenient showing of cross sections
% fprintf('IBOUND:   '); XS(IBOUND)
% fprintf('HK:       '); XS(HK    )
% fprintf('VK:       '); XS(VK    )
% fprintf('PEFF:     '); XS(PEFF  )
% fprintf('ICBUND:   '); XS(ICBUND)
% fprintf('STCONC:   '); XS(STCONC)
% fprintf('STRTHD:   '); XS(STRTHD)

%% Grid

[xGr,yGr,zGr,xm,ym,zm,DX,DY,DZ,NX,NY,NZ]=modelsize3(xGr,yGr,zGr);

%% Get unformatted data

H=readDat([basename,'.hds']);
C=readMT3D('MT3D001.UCN');
B=readBud([basename,'.bgt'],'FLOWRIGHTFACE');  % get only flow rightface

m=0; M=0;
for i=1:length(C)
    m=min([m;C(i).values(:)]);
    M=max([M;C(i).values(:)]);
end
fprintf('min=%g Max=%g\n',m,M');
caxis([m M]);


%% Stream function

NT=length(B);

Psi=NaN(NZ+1,NX-1,NT);
psi=XS(B(end).term{1}(:,1:end-1,:));  % cut off east most cell
for it=1:length(B)
    Psi(:,:,it)=[flipdim(cumsum(flipdim(psi,1)),1); zeros(size(psi(1,:)))];
end

%% Contour each result or make movied
figure;

load Cbrine

isolines=[0.01 0.02 0.05 0.1 0.2 0.3 0.5]*Cbrine;

if exist('Film','var'); clear('Film'); end
Film(NT)=getframe;  % allocate
for it=1:NT
    hold off
    [c,h]=contour(xm,squeeze(zm),XS(C(it).values),isolines,'linewidth',2);
    hold on
    clabel(c,h);
    contour(xGr(2:end-1),squeeze(zGr),Psi(:,:,it));
    ht =title(sprintf('Hydrocoin problem, t=%.0f',C(it).time));
    xlabel('x [m]');
    ylabel('z [m]');
    Film(it) = getframe;
 end
 movie(Film)
