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

problem='2D-Diagonal flow field.';
%% visualization

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

%% 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={
    'no sorption'
    'Linear'
    'Freundlich'
    'Langmuir'
    'First order kinetic'
    'Dual domain mass transfer'
    'Dual dom. mass transf. +sorption'};

%% Get advection method used from MT3D worksheet
[MT3Dparnams,MT3Dparvals]=getExcelData([basename,'.xls'],'MT3D','vertical');
MIXELM=MT3Dparvals(strmatchi('MIXELM',MT3Dparnams));

%%

figure; hold on;

[c,h]=contourf(xm,ym,C(end).values,[0.1,0.2 0.4 0.8 1.6 3.2 6.4]);
clabel(c,h);

title(sprintf('%s t=%g d, advection solution method = %s',...
    problem,C(end).time,ADVmethod{MIXELM+2}));
xlabel('xGr [m]'); ylabel('yGr [m]');
grid('on');

%% Observation points

figure; hold on;

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

Obs=readOBS(d(1).name);
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);
title(sprintf('%s, observation points, advection solution method = %s',...
    problem,ADVmethod{MIXELM+2}));
xlabel('time [s]'); ylabel('Rel. conc mg/l'); grid on
