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

%% visualization

[xGr,yGr,xm,ym,Dx,Dy,Nx,Ny]=modelsize(xGr,yGr);

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

%% Concentrations from files, each file is for one species
%%(we have only one in this example though

C=readMT3D(sprintf('MT3D001.UCN'));

%% B=readBud([basename,'.bgt'],'FLOWRIGHTFACE');

%% Visualization
ADVmethod={
    'ULTIMATE (TVD)'
    'FDM standard'
    'MOC'
    'MMOC'
    'HMOC'};
    
isothm={
    'Linear'
    'Freundlich'
    'Langmuir'
    'First order kinetic'
    'Dual domain mass transfer'
    'Dual dom. mass transf. +sorption'};

%% Get MT3D parameters of this run, as we need MIXELM in figure title
[MT3Dparnams,MT3Dparvals]=getExcelData([basename,'.xls'],'MT3D','vertical');
MIXELM=MT3Dparvals(strmatchi('MIXELM',MT3Dparnams));

%%
figure; hold on;

[c,h]=contour(xm,ym,C(end).values,[0.1,1,10,50]);
clabel(c,h);

title(sprintf(['Well in uniform flow field, Zheng (1999), p139, t=%g d, ',...
    'Solution method = %s'],C(end).time,ADVmethod{MIXELM+2}));
xlabel('xGr [m]'); ylabel('yGr [m]');
grid('on');

%% Observation points
% to show the registered concentration at the concentration points
figure; hold on;

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

for iSpecies=1:length(d)  % for all species
    Obs=readOBS(d(iSpecies).name);  % get obs points conc for species

    leg=cell(Obs.NOBS,1);
    for iObs=1:Obs.NOBS
        leg{iObs}=sprintf('L=%d R=%d C=%d',Obs.Loc(iObs,:));
    end

    plot(Obs.t,Obs.values);  % relative to Co, right-most column is exit point
    legend(leg,2);

    title(sprintf('Conc at obs points, solution method = %s',ADVmethod{MIXELM+2}));
    xlabel('time [d]'); ylabel('Conc mg/l'); grid on
end