%% Analyzing model output

load('name'); load(basename);  % Get basename and parameters generated by mf_adapt

%% Retrieve the model grid

xm=0.5*(xGr(1:end-1)+xGr(2:end));
ym=0.5*(yGr(1:end-1)+yGr(2:end));
ZM=0.5*(Z(:,:,1:end-1)+Z(:,:,2:end));
zm=squeeze(ZM(1,1,:));
z =squeeze(Z (1,1,:));

%% Reading unformattted files
H=readDat([basename,'.hds']);  % just have a look (can also look in the LST file

%% Get computed concentrations from files (all species (even if only 1))

% Put the concentrations in a structarray of size NCOMP
C=cell(100,1);  % initialize space for 100 species
for iComp=1:100
    try
        % get struct for species iComp
        C{iComp}=readMT3D(sprintf('MT3D00%d.UCN',iComp));
    catch ME
        NCOMP=iComp-1;
        C=C(1:NCOMP);  % do away with unneeded space
        break;
    end
end


%% B=readBud([basename,'.bgt'],'FLOWRIGHTFACE'); % not used for this problem

%% Visualization of the advancement of the transprot of all species
% subject to dispersion, sorption and decay
% See the header of the figure for the detauls, note that the retardation
% R=1+Rhob*Kd/peff, which, when we set Rhob=1, Kd=1 and peff =0.25, yields R=5

clr='brgkmc'; mrk='^sp*Vx+';

peff=0.25; % taken from mf_adapt !!!

%% Get layer parameters used from workbook to construct a useful figure title

[LAYparnams,LAYparvals,LAYCBD,NLAY]=getLayers([basename '.xls'],Z);
aL     =LAYparvals(:,strmatchi('AL'  ,LAYparnams));
Rho    =LAYparvals(:,strmatchi('RhoB',LAYparnams));
Kd     =LAYparvals(:,strmatchi('SP1_1',LAYparnams));
Lambda =LAYparvals(:,strmatchi('RC1_1',LAYparnams));
R=1+Rho.*Kd./peff;

%% Then plot the results for the 4 layers in separte subplots

for iComp=1:NCOMP
    figure
    sTitle =sprintf('Species %d',iComp);
    tTitle =sprintf('time=%g to %g d',C{iComp}(1).time,C{iComp}(end).time);
    t2Title=sprintf('in %d steps',length(C{iComp}));
    
    for iLay=1:NLAY
        subplot(NLAY,1,iLay); hold on
        lTitle =sprintf('Layer %d',iLay);
        l2Title=sprintf('aL=%gm, R=%g, Lambda=%g/d',aL(iLay),R(iLay),Lambda(iLay));
            for it=1:length(C{iComp})
                plot(xm,C{iComp}(it).values(1,:,iLay),[clr(iComp),'-']);
            end
        set(gca,'xlim',xGr([1 end]),'ylim',[0 1],'xgrid','on','ygrid','on');
        title(['Problem ',basename,', ',sTitle,' ',lTitle,', ', ...
            tTitle,' ',t2Title,', ',l2Title]);
        xlabel('xGr [m]');
        ylabel('Relative conc');
    end
end

% That's all, compare with the manual
% TO 091202