%% Analyzing output of this model tracer dilution test
% TO 091205

% Interpret the results of the test. Clearly, the method of injection of
% the tracer either as uniform distribution of increase concentration at
% the entry face of the model or as mass_loading does change the tracer
% distribution in the model at different times, but does not alter the
% method of interpretation. Therefore, this intepreting script does not
% need any reference to the question on how the tracer was injected (but we
% should make sure that in either case the same overall mass has been
% injected (which we did, see mf_adapt).
%
% Clearly, the mass_laoding the tracer distribution is moredepdendent on the
% persion in all three directions. Without dispersion, there is no way for
% the mass_loaded tracer to reach all cells of the model.
%

%% load model name and basename contained in name.mat
load name
load(basename);  % this yields the stored matrices that make up the grid and the model

[NROW,NCOL,NLAY]=size(IBOUND);

%% compute coordinates of axes and cell centers

[xGr,yGr,xm,ym,DELX,DELY,NCOL,NROW]=modelsize(xGr,yGr); % grid housekeeping

ZM=0.5*(Z(:,:,1:end-1)+Z(:,:,2:end));   % elevation of cell centers

[XM,YM,ZM]=meshgrid(xm,ym,(ZM(1,1,:))); % cell centers for 3D visualization  
[DX,DY,DZ]=meshgrid(DELX,DELY,diff(Z(1,1,:)));

%% load the unformatted files with the heads, the concentration and the budget terms

% if ~exist('H','var')
%     H=readDat([basename,interm,'.hds']); H=maskHC(H,IBOUND);   % read heads and mask inactive cells with NaN
%     for i=1:length(H)
%         H(i).values(IBOUND==0)=NaN;
%     end
% end

if ~exist('C','var')
    C=readMT3D('MT3D001.UCN');     C=maskHC(C,(ICBUND & IBOUND));   % read concentrations and mask inactive with NaN
end

if ~exist('B','var')
  B=readBud([basename,'.bgt']);  B=maskHC(B,IBOUND); % read budgetfile and mask
end

%% make movie
film =0;

cM=max(max(max(~isnan([C.values])))); dc=cM/10;
%isovalue=dc:dc:cM-dc;
if film
        mov=avifile(basename,'compression','none',...
        'fps',15,'colormap',colormap,...
        'videoname','Mt3d');
else
    figure;
end

h=title(sprintf('t=%g s',0));
xlabel('xGr [m]'); ylabel('yGr [m]'); zlabel('z [m]');

for i=1:length(C)
   %[conc
   %t]=plotCNC(xm,zm,C,period,tstp,'row',ixSec,isovalue,'w',@contourf,'none');    % plot the concentrations
   
   slice(XM,YM,ZM,C(i).values,[],0,0); % slice is nice but colors cannot be managed
   
   title(sprintf('t=%g s',C(i).time));
   drawnow;

   if film
       F=getframe(gcf);
       mov = addframe(mov,F);
   end
end

if film
    mov=close(mov);
end


%% read Observaton point concentration output
% We have defined some observation points in the Worksheet BTNOBS. These
% point require no package and will be used automatically when specified.
% MT3DMS generates files MT3D???.OBS where ??? is the species number
% The data of all observation points are in the same file, one file per
% species

d=dir('MT3D*.OBS');  % get names of all files with observations

clrs='brgkmc';  % colors to be used in plot

for iObs=1:length(d)
    figure; hold on;
    Obs=ReadObs(d(iObs).name);     % reads file into Obs struct
    Obs.values(Obs.values<0)=NaN;  % mask inactive cells
    
    clr=repmat(clrs,1,ceil(Obs.NOBS/length(clrs)));  % make sure there are enough colors
    
    leg=cell(Obs.NOBS,1);   % initialize legend texts
    for iLoc=1:Obs.NOBS     % for all obsrvation points (locations)
        leg{iLoc}=sprintf('pnt [%2d,%2d,%2d]',Obs.Loc(iLoc,:));
        plot(Obs.t,Obs.values(:,iLoc),clr(iLoc));
    end
    legend(leg);
    title('Observation point concentrations');
    xlabel('time [s]'); ylabel('concentration'); grid on
end

%% Now compute back the flow by the observations

% First get the applied mass
load Q0.mat  % contains Q0 and M as saved in mf_adapt

% The computed flow from the salt dilution test equals the intgral of the
% tracer concentation over the time from the start of the test ot infinity,
% divied by the total injected mass
FlowsRetrieved=M./trapz(Obs.t,Obs.values);

% Give the result
fprintf('The dischage (flow [m3/s]) computed form the obs points is\n');
fprintf('  %g',FlowsRetrieved);
fprintf('\n');
   
% Now you may experiment by placing the observation wells at different
% lcoations.

%% Of course, one could also compute the retrieved fraction from every cell
% of the model. This is usefull to investigate at what possible observation
% points the salt dilution test would fail due to the flow distribution in
% the system, i.e. in the box and the end pipe.
% To compute this, we extract from the model output the tracer
% concentration over time for all active cells

% for plotting results get
Iactive=find(IBOUND~=0);  % get all active cells of model

NT  =length(C);                % the number of time steps saved in concentration file
Data=NaN(length(Iactive),NT);  % allocate to store conc data for all active cells
for it=1:NT
    Data(:,it)=C(it).values(Iactive);  % fill this matrix
end
ct=trapz([C.time],Data,2);    % integrate over time using trapezium rule
Qt=M./ct;                     % this is the retrieved discharge that we
                             % we would get by sensoring any model cell
% because in the test we know the real Q, i.e. sum(Q0) from mf_adapt and
% saved in Q0.mat together with total mass M, we can now conveniently
% compute the fraction of the original discharge that we would measure by
% doing a tracer dilution test in any cell of the model:

fprintf(['Minimum relative flow from test = %.3f\n',...
         'Maximum relative flow form test = %.3f\n'],...
         [min(Qt),max(Qt)]/sum(Q0));
     
 %% Visualization in 3D

 % first put the rsults back into their 3D position
Qrel =NaN(size(IBOUND)); Qrel(Iactive)=Qt/sum(Q0);

figure; hist(Qrel); xlabel('Qrel'); ylabel('frequency');

%% Show this in 3D using transparent isosurfaces, which shows the relative
% discharge we would obtain if we put our sensor in any cell
isovalue=0:0.1:10; isovalue(isovalue<min(Qrel(:)) | isovalue>max(Qrel(:)))=[];
    
figure; xlabel('x [m]'); ylabel('y [m]'); zlabel('z [m]');
title('3D view on relative discharge measured at any point from dilution test');

% first show the walls of the pipe and box by contouring its boundary,
% which is easily found by IBOUND==0 and coutouring value 0.5
patch(isosurface(XM,YM,ZM,IBOUND==0,0.5),'edgecolor','b','facecolor','b');
view(3);     % convenient standard 3D view

% Than plot the isosurfaces for the relative disccharge
for i=1:length(isovalue);
    patch(isosurface(XM,YM,ZM,Qrel,isovalue(i)),'facecolor','r','edgecolor','none');
    hold on;
end
alpha(0.3); % make isosurfaces trasnparent
        
%% Total mass in the model during the simulation
I=IBOUND~=0;
fprintf('\nTotal tracer [g] in the model over time:');
for i=1:length(C)
 if rem(i,10)==1; fprintf('\nTime step %3d Mass=',i); end
     fprintf(' %12.3f',sum(C(i).values(I).*abs(DX(I).*DY(I).*DZ(I))));
end
fprintf('\n');

%% Max conc in the model after say 20 s
fprintf('\nMax tracer conc [g/l] in the model over time:');
for i=1:length(C)
 if rem(i,10)==1; fprintf('\nTime step %3d Conc=',i); end
     fprintf(' %12.3f',max(C(i).values(I)));
end
fprintf('\n');
