%% mf_analyze, print and plot results fo simulation

%% Analyzing output of xsecton model
% Plot /shows the results of the smulation for a given stress period and
% time step to be set by the user
% TO 090223
clear variables
close all

%% 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
load underneath

[xGr,yGr,Z,xm,ym,zm,Dx,Dy,Dy,Nx,Ny,Nz]=modelsize3(xGr,yGr,zGr);

film=1;  % if film==1, make avi file else not

%% load the unformatted files with the heads, the concentration and the budget terms
interm='';
H=readDat([basename,interm,'.hds']);
H=maskHC(H,IBOUND & H(1).values<1000);   % read heads and mask inactive cells with NaN

C=readMT3D(['MT3D001',interm,'.UCN']);
C=maskHC(C,IBOUND);   % read concentrations and mask inactive with NaN
C=maskHC(C,ICBUND);   % read concentrations and mask inactive with NaN

B=readBud([basename,interm,'.bgt'],'FLOWRIGHTFACE');  B=maskHC(B,IBOUND); % read budgetfile and mask
for i=1:length(B)
  for j=1:length(B(i).term)
      B(i).term{j}(IBOUND==0 | H(1).values>1000)=NaN;
  end
end

%% Photo to plot on
foto='Series1_05_p18.jpg';
load(['.' filesep 'Photos' filesep 'photostruct.mat']); % yields struct d
i=strmatchi(foto,{d.name});
fig=ImageModel(['Photos' filesep foto],d(i).uMdl,d(i).vMdl,d(i).DX,d(i).DY); 

%set(fig,'position',get(0,'ScreenSize'));

%% Set up contours and first plot to initialize animation

ifrf=strmatch('FLOWRIGHTFACE',B(1).label);
minPsi=Inf; maxPsi=-Inf;
for i=1:length(B)
    FRF=XS(B(1).term{ifrf});  % Flow right face
    Psi=flipud([zeros(size(FRF(1,1:end-1))); cumsum(flipud(FRF(:,1:end-1)))]);
    minPsi=min([minPsi;Psi(:)]); maxPsi=max([maxPsi;Psi(:)]);
end

concs= 0.05:0.1:0.95;      % Conc is always between 0 and 1
chead=40:0.05:100;      % Heads is always between base and top of model
cpsi =minPsi:(maxPsi-minPsi)/50:maxPsi; % Psi-range is computed above

% Stream function mask, to not plot stream lines above sand body
PSIMASK=NaN(Nz+1,Nx-1);
PSIMASK(inpolygon(ones(size(zGr(:)))*xGr(2:end-1),zGr(:)*ones(size(xGr(2:end-1))),xSand,zSand))=1;
PSIMASK(inpolygon(ones(size(zGr(:)))*xGr(2:end-1),zGr(:)*ones(size(xGr(2:end-1))),xCanL,zCanL))=NaN;

if exist('M','var'), clear('M'); end
if film
    if ismac
        aviobj=avifile(basename,'compression','none',...
            'fps',10,...
            'videoname','Pennink, series 5, March 1905');
    else
        aviobj=avifile(basename,'compression','Cinepak',...
            'fps',10,...
            'videoname','Pennink, series 5, March 1905');
        end
end
hold on;

ifrf=strmatch('FLOWRIGHTFACE',B(1).label);
FRF=XS(B(1).term{ifrf});  % Flow right face
Psi=flipud([zeros(size(FRF(1,1:end-1))); cumsum(flipud(FRF(:,1:end-1)))]).*PSIMASK;

[c,hc]=contourf(xm,squeeze(zm),XS(C(1).values),concs,'g');
[c,hh]=contour (xm,squeeze(zm),XS(H(1).values),chead,'r');
[c,hp]=contour (xGr(2:end-1),squeeze(zGr),Psi,cpsi,'w');

title(sprintf('time = %g min',C(1).time));
xlabel('x [m]');                 % set x-label
ylabel('z [m]');                 % set y-label
set(gca','color',[1 1 0.8])      % light blue background color

drawnow

F=getframe(fig); 

if film, aviobj=addframe(aviobj,F); end

%% Show all the figures and put them in the avi file

for i=2:length(C)
    period=C(i).period;
    tstp  =C(i).tstp;                   % choose a timestep period present in the file
    
    ih=find([H.period]==period & [H.tstp]==tstp);
    ib=find([B.period]==period * [B.tstp]==tstp);
    set(hh,'zdata',XS(H(ih).values));
    set(hc,'zdata',XS(C(i).values));
    
    FRF=XS(B(ib).term{1}); FRF(isnan(FRF))=0;
    Psi=flipud([zeros(size(FRF(1,1:end-1))); cumsum(flipud(FRF(:,1:end-1)))]).*PSIMASK;
    set(hp,'zdata',Psi);
    
    title(sprintf('time = %g min',C(i).time));

   drawnow

    F = getframe(fig);
    if film, aviobj = addframe(aviobj,F); end
end

if film, aviobj=close(aviobj); end
%save movie.mat F