%% mf_analyze, print and plot results fo simulation
% TO 090223 100718

clear variables
close all

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

load name
load(basename);  % this yields the stored matrices that make up the grid and the model
load underneath

%% Start


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

isCan=IBOUND<0;


%% load the unformatted files with the heads, the concentration and the budget terms
if ~exist('H','var')
    H=readDat([basename,'','.hds']);
    H=maskHC(H,IBOUND & H(1).values<1000);   % read heads and mask inactive cells with NaN
end
if ~exist('C','var')
    C=readMT3D(['MT3D001','','.UCN']);
    C=maskHC(C,IBOUND);   % read concentrations and mask inactive with NaN
    C=maskHC(C,ICBUND);   % read concentrations and mask inactive with NaN
end
if ~exist('B','var')
  B=readBud([basename,'','.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
end

%% 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(i).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

% Stream function mask, to not plot stream lines above sand body
PSIMASK=(ones(size(zGr,3),1)*interp1(xTop,zTop,xGr(2:end-1))>squeeze(zGr)*ones(size(xGr(2:end-1))));

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

%% Get the photo to plot onto
foto='Series4_01_p64.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'));

% Start plotting
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

FRF=XS(B(1).term{strmatch('FLOWRIGHTFACE',B(1).label)});
Psi=flipud([zeros(size(FRF(1,1:end-1))); cumsum(flipud(FRF(:,1:end-1)))]);

[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),maskHC(Psi,PSIMASK),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 from same stress period as C
    ib=find([B.period]==period);  % B from same stress period as C
    set(hh,'zdata',XS(H(ih).values));  % recontour h
    set(hc,'zdata',XS(C(i).values));   % contour   C
    
    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)))]);
    set(hp,'zdata',maskHC(Psi,PSIMASK));  % contour stream function
    
    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