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

% after loading the basename and basename.mat we have retrieved all
% parameters used to assemble the model. This allows cleaning the workspace
% before the load command.

%% retrieve grid nessesities

[xGr,yGr,xm,ym,DELX,DELY,Nx,Ny]=modelsize(xGr,yGr);

[LAYparnams,LAYparvals]=getExcelData([basename,'.xls'],'LAY');
beta=LAYparvals(:,strmatchi('SP2_1',LAYparnams));

%% Reading unformattted files
%H=readDat([basename,'.hds']);

%% Relative concentration (given in PNTSRC of SSM module)
% To get relative concentrations, it may be more practicle to use an input
% concentration of 1. But then make sure this does not interfere with the
% non-linear concentrion dependent behavior
C0=0.05;  % max concentration to compute relative concentration

%% Concentrations from files, each file is for one species
% (we have only one in this example though)
% But in general w may have an almost arbitrary number (but no more than
% one hundred) of different species,
% whose results MT3DMS saved in files MT3Dxxx.UCN where xxx is the species'
% number 001 002 003 etc.
%
% To read in at most the files of hundred species:
%
C=cell(100,1);
for iComp=1:100
    try  % if this fails we have read in the last species conc file
        C{iComp}=readMT3D(sprintf('MT3D00%d.UCN',iComp));
    catch ME  % which is ok, so remove the non-used cells in C
        C=C(1:iComp-1);
        break;
    end
end

% Now C is a cell array of lenght equal to the nr of species.

NT=length(C{1});  % number of transprot steps in these files
NCOMP=length(C);  % number of species

%% B=readBud([basename,'.bgt'],'FLOWRIGHTFACE');
% Don't need this here
%
%% Visualization

isothm={
    'no sorption'
    'Linear'
    'Freundlich'
    'Langmuir'
    'First order kinetic'
    'Dual domain mass transfer'
    'Dual dom. mass transf. +sorption'};

clr='brgkmc';   % define line colors
mrk='^sp*Vx+';  mrk=repmat(mrk,1,ceil(NT/length(mrk)));  % enough markers for curves

for iComp=1:NCOMP
    figure; hold on
    ttile=sprintf('Species %d, time from %g to %g s',...
        iComp,C{iComp}(1).time,C{iComp}(end).time);
        time=[C{iComp}.time];
        data=reshape([C{iComp}.values],[Nx,NT,NLAY]);
        data=data/C0; % Relative concentration
        for iLay=1:NLAY
            subplot(NLAY,1,iLay); hold on
            for it=1:NT
                ic=rem(it,length(clr)); if ic==0, ic=4; end 
                plot(xm,data(:,it,iLay),clr(ic));
            end
            title(sprintf('%s %d species process=%s SP1=%g, SP2=%g, time=%.0f to %.0f s',...
                basename,NCOMP,isothm{rct.ISOTHM+1},rct.SP1{iLay},rct.SP2{iLay},...
                C{iComp}(1).time,C{iComp}(end).time));
            xlabel('x [cm]'); ylabel('rel conc');
            grid('on');            
        end
        
        figure; hold on
        ix=find(xm>=8,1,'first');
        leg=cell(NLAY,1);
        for iLay=1:NLAY
            plot(time,permute(data(ix,:,iLay),[2,3,1]),clr(iLay));
            leg{iLay}=sprintf('beta=%g',beta(iLay));
        end
        legend(leg);
        title(sprintf('%s breakthrough at x=%.0f of species Nr %d, process=%s',...
            basename,xm(ix),iComp,isothm{rct.ISOTHM+1}));
        xlabel('time [s]'); ylabel(sprintf('Concentration at xm=%.0f cm',xm(ix)));
        grid('on');
        
end


%% Observation points
% this is the alternative to the approach above, we have defined
% an obsrvation point in the worksheet, so MT3DMS saved the concentration
% over time for this point in the file MT3D*.OBS

figure; hold on;

d=dir('MT3D*.OBS');

Obs=readOBS(d(1).name);

plot(Obs.t,Obs.values/C0);  % relative to Co, right-most column is exit point
title(sprintf('Obs point at x=%g cm, iX=%d, process=%s',...
    xm(Obs.Loc(end,3)),Obs.Loc(end,3),isothm{rct.ISOTHM+1}));
xlabel('time [s]'); ylabel('Rel. conc mg/l'); grid on

