%% mf_analyze, print and plot results fo simulation
% Pennink Series6
%
% A (compressed) avi file can be generated by setting film=1;
% The already made avi files with this modle can be found on the
% download section of http://code.google.com/p/mflab under
% avifiles Pennink/Series6
%
% Compressed avi fles can only be made under Windows (due to Matlab’s lack of
% compression support on other operating systems. By default I set off avifile
% generation on the Mac because the uncompressed avifiles are just huge and a real
% burden to the system, where under Windows with compression they are about 10 MB,
% which is reasonable.
%
%
% TO 090223 100721 100726

clear variables
close all

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

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

%% Start

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

%% load the unformatted files with the heads, the concentration and the budget terms
interm='';

H=readDat([basename,interm,'.hds']);
H=maskHC(H,IBOUND);   % 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']);  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)=0;  % set flow terms inactive cells 0
  end
end

PsiMax=-Inf;
PsiMin=+Inf;

ifrf=strmatch('FLOWRIGHTFACE',B(1).label);
ichd=strmatch('CONSTANTHEAD', B(1).label);
ighb=strmatch('HEADDEPBOUNDS',B(1).label);

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)))]);
    PsiMax=max(PsiMax,max(Psi(:)));
    PsiMin=min(PsiMin,min(Psi(:)));
end

% Stream function mask, to not plot stream lines above sand body
PSIMASK=NaN(size(zGr(:)*xGr(2:end-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;
PSIMASK(inpolygon(ones(size(zGr(:)))*xGr(2:end-1),zGr(:)*ones(size(xGr(2:end-1))),xCanR,zCanR))=NaN;
%PSIMASK(inpolygon(ones(size(zGr(:)))*xGr(2:end-1),zGr(:)*ones(size(xGr(2:end-1))),xCanM,zCanM))=NaN;

%% Set up contours and first plot to initialize animation

concs= 0.05:0.1:0.95;
chead=70:0.25:100;
cpsi =PsiMin:PsiMax/50:PsiMax;

fig=figure('position',get(0,'screensize')); set(gca,'cLim',[0 1]); hold on
set(gca,'xLim',xGr([1 end]),'yLim',[zGr(end) zGr(1)])
patch(xSand,zSand,[1 1 0.6]);
plot(xSat,zSat,'b');
patch(xCanL,zCanL,'w');
patch(xCanR,zCanR,'w');
%patch(xCanM,zCanM,'w');
patch(xWatL,zWatL,'b');
patch(xWatR,zWatR,'b');
%patch(xWatM,zWatM,'b');
patch(xWel ,zWel   ,'b');


if film
    aviobj=avifile(basename,'compression','Cinepak',...
        'fps',5,...
        'videoname','Pennink, series 5, March 1905');
end
hold on;

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,zm,XS(C(1).values),concs,'g');
[c,hh]=contour (xm,zm,XS(H(1).values),chead,'r');
[c,hp]=contour (xGr(2:end-1),squeeze(Z),Psi,cpsi,'k');

title(sprintf('time = %g min',C(1).time));
xlabel('x [m]');                 % set x-label
ylabel('z [m]');                 % set y-label
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=H(i).period;
    tstp  =H(i).tstp;                   % choose a period 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{ifrf});
    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


%% Discharge through boundary conditions
IL=cellIndex(IL(:,3),IL(:,2),IL(:,1),Nx,Ny,Nz);
IR=cellIndex(IR(:,3),IR(:,2),IR(:,1),Nx,Ny,Nz);
IM=cellIndex(IM(:,3),IM(:,2),IM(:,1),Nx,Ny,Nz);
IW=cellIndex(IW(:,3),IW(:,2),IW(:,1),Nx,Ny,Nz);
IS=cellIndex(IS(:,3),IS(:,2),IS(:,1),Nx,Ny,Nz);

fprintf('%10s %10s %10s %10s %10s %10s %10s\n','Period','tstep','QL','QM','QR','QS','Qtot');
ccmpm2lph=1e-6*60; % m3/min to L / h
for i=1:length(B) 
    QL=sum(B(i).term{ighb}(IL))*ccmpm2lph;
    QM=sum(B(i).term{ighb}(IM))*ccmpm2lph;
    QR=sum(B(i).term{ighb}(IR))*ccmpm2lph;
    QS=sum(B(i).term{ighb}(IS))*ccmpm2lph;
    fprintf('%10d %10d %10g %10g %10g %10g %10g\n',B(i).period,B(i).tstp,QL,QM,QR,QS,sum([QL QM QR QS]));
end


%% Interface computation
%  The head in the milk is set in the input. It serves to keep the
%  interface at the desired elevation. However, the correct milk-head is
%  known only afterwards and iteration is cumbersum. So we compute here the
%  desired milk-head given the freshwater head and the saltwater head at
%  two strategic positions in the model (below the right ditch at z=50 cm,
%  which should be a good position for the freshwater head without much
%  influence of the converging flow near te ditch and still fresh, and the
%  pint near the milk supply. If the correct head is chosen, the interface
%  should be stable and little milk flow should occur and so this saltwater
%  head should be more or less that in the milk reservoir).
rhof=1000; rhos=1029; zIface=23;  % must match what has been specified in model

[ix,iy,iz]=xyzindex([95 95]',[1 1]',[50 5]',xGr,yGr,zGr(:)); % point positions
I=cellIndex(ix,iy,iz,Nx,Ny,Nz); % position indices in data array

phif=NaN(size(H));  % freshwater head in first  point
phis=NaN(size(H));  % saltwater  head in second point 
for i=1:length(H)
    phif(i)=H(i).values(I(1));
    phis(i)=H(i).values(I(2));
end
figure; hold on

subplot(2,1,1);  % plot phif, phis and phis desired to maintain zIface
plot([H.totim],phif,'b',...
     [H.totim],phis,'r',...
     [H.totim],(rhos-rhof)/rhos*zIface+rhof/rhos.*phif,'g');
legend('phi_f(50)','phi_s(5)','phi_s desired');
xlabel time; ylabel head; title('fresh and saltwater head and satwater head need to keep interface at zIface');

subplot(2,1,2); % plot computed interface form heads and desired zIface
plot([H.totim],(rhos*phis-rhof*phif)/(rhos-rhof),'r',...
     [H.totim],zIface*ones(1,length(H)),'g');
legend('computed interface elevation','desired elevation zIface')
xlabel time; ylabel head; title('computed interface given the heads in the model');
