%% Analyzing model output
load('name'); load(basename);
load('cIn');  % stores the injection concetration from mf_adapt

problem='2D-radial transport';
%% 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={
    'FDM standard'
    'Linear'
    'Freundlich'
    'Langmuir'
    'First order kinetic'
    'Dual domain mass transfer'
    'Dual dom. mass transf. +sorption'};

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

%% contouring etc
figure; hold on;

clim=[0 1];
[c,h]=contourf(xm,ym,C(end).values);
clabel(c,h);
colormap(flipud(bone));
colorbar;

title([problem,sprintf(' t=%g d, TVD',C(end).time)]);
xlabel('xGr [m]'); ylabel('yGr [m]');
grid('on');

%% Graph along an axis
figure, hold on
plot(xm(xm>=0),C(end).values(ym==0,xm>=0,1)/cIn)

title([problem,sprintf(' t=%g d, TVD',C(end).time)]);
xlabel('radial distance [m]'); ylabel('C/Co');
grid('on');


%% Notice

s='To refine curve near well half DELX and DELY and double NCOL and NROW in mf\_adapt';

xlim=get(gca,'xlim'); ylim=get(gca,'ylim');
text(xlim(1)+0.05*diff(xlim),ylim(1)+0.7*diff(ylim),s);

