%% Analyzing model output
load('name'); load(basename);
load underneath  % extra parameters from mf_adapt

%% 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');
%C=readMT3D('MT3D001.UCN',[],1,1,1:100);
B=readBud([basename,'.bgt'],'FLOWRIGHTFACE');  % get only flow rightface

NT=length(C);

%% Contour each result or make movie

figure;  hold on; grid on

clr=repmat('rbgkmc',1,100);

set(gca,'ylim',[0 C0max],'xlim',[0 xGr(end)]);

leg={};
for it=1:NT
    plot(xm,C(it).values(1,:,1),clr(it));
    leg{end+1}=sprintf('%.0f y',C(it).time/365);
    %leg{end+1}=sprintf('%g d',C(it).time);
  
    xmrk=xm(1:5:end);

    sigma=sqrt(2*D/R*C(it).time);
    switch iCase
        case {1,2},         T=C0max*erfc(1/sqrt(2));
            plot(sigma,C0max*erfc(1/sqrt(2)),['*' clr(it)]);
        case 3,             %T=C0max*exp(0.5)*peff*R*(W/2)/(sigma*sqrt(2*pi));
            T=C0max/2*(erfc((xmrk-W/2)/(sigma*sqrt(2)))-erfc((xmrk+W/2)/(sigma*sqrt(2))));
            plot(xmrk,T,['.' clr(it)]);
        case 4,
            T=C0max*(W)/(sigma*sqrt(2*pi))*exp(-xmrk.^2/(2*sigma^2));            
            plot(xmrk,T,['.' clr(it)]);
       case 5,
            T=C0max*(W)/(sigma*sqrt(2*pi))*exp(-xmrk.^2/(2*sigma^2));            
            plot(xmrk,T,['.' clr(it)]);
    end
    leg{end+1}=sprintf('%.0f y',C(it).time/365);
    %leg{end+1}=sprintf('%g d',C(it).time);
end
legend(leg);

xlabel('x [m]');
switch iCase
    case 1, ylabel('conc [ug/L]');
            title('Diffusion');
    case {2,3,4}, ylabel('temp [C]');
            title('Heat flow');
    case {5}, ylabel('temp [C]');
            title('Heat flow due to heat injection pulse');
end
 