%% Analyzing output of polder1 model
% TO 091226

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

load pol  % saved from mf_adapt: R RP Iout Iin hIn hOut

[xGr,yGr,zGr,xm,ym,zm,dx,dy,dz,Nx,Ny,Nz]=modelsize3(xGr,yGr,zGr);
[DX,DY]=meshgrid(dx,dy);


%% load the unformatted files with the heads, the concentration and the budget terms
H=readDat([basename,'.hds']); H=maskHC(H,IBOUND);   % read heads and mask inactive cells with NaN
C=readMT3D('MT3D001.UCN');    C=maskHC(C,ICBUND);
B=readBud([basename,'.bgt']);

t=[H.totim]; Nt=length(t);  % check that length B equals length C and length H

%% remove intermediate concentration times
tEPS=1e-3;
tc=[C.time];
Itc=NaN(size(t));
for it=1:length(t)
    Itc(it)=find(tc>=t(it)-tEPS,1,'first');
end
C=C(Itc);


%% Indices of top layer inside and outside polder
Iin =Iin( Iin <=Nx*Ny);
Iout=Iout(Iout<=Nx*Ny);

%% Total discharge and total fresh discharge


iGHB=strmatch('HEADDEPBOUNDS',B(1).label);

% Dit moet nog uitgewerkt worden
% Seepage flow entire polder
Qpol =NaN(Nt,1);
Qsalt=NaN(Nt,1);
Lpol =NaN(Nt,1);
cpol =NaN(Nt,1);
Asalt=NaN(Nt,1);

% Xsec through polder at ym=0
iy=find(ym<100,1,'first');
IX=find(abs(xm)<RP);          % inside polder
qcross=NaN(Nt,length(IX));
ccross=NaN(Nt,length(IX));

cMax = max(C(1).values(:));
A    =sum(DX(Iin).*DY(Iin)); fprintf('Polder area = %.0f ha\n',A/1e4);

for it=1:Nt
    Qpol(it)=-sum(B(it).term{iGHB}(Iin));
    Lpol(it)=-sum(B(it).term{iGHB}(Iin).*C(it).values(Iin));
    cpol(it)= Lpol(it)/Qpol(it);
    Qsalt(it)=cpol(it)/cMax*Qpol(it);
    I=find(C(it).values(Iin)>cMax/2);
    Asalt(it)=sum(DX(Iin(I)).*DY(Iin(I)));    
    qcross(it,:)=-B(it).term{iGHB}(iy,IX,1)./(dy(iy)*dx(IX));
    ccross(it,:)= C(it).values(iy,IX,1);  % mg/L
end

%% Visualize results

%% Total flow and fresh flow out of polder
% Graph of total polder discharge, dicharge-load and discharge
% concentration versus simulation time.
figure

dpy=365.25;   % days per year, to show in years instead of days
ty=t/dpy;     % time in years

subplot(4,1,1); hold on;
plot(ty, Qpol*dpy/1e6,'b');  % discharge in million m3/year
plot(ty,Qsalt*dpy/1e6,'r');  % saltwater discharge  m3/year
legend('Qtotal','Qsalt');
ylabel('Mm^3/y'); title('Total polder and total saltwater discharge');
grid on

subplot(4,1,2); hold on;
plot(ty,Qsalt/Qpol,'b');  % Salt discharge fraction m3/year
plot(ty,Asalt/A   ,'r');     % Polder area fraction with c>0.5CMax
set(gca,'ylim',[0 1]);
legend('Qsalt/Q','Asalt/A');
ylabel('fraction'); title('Relative salt discharge and relative salt area');
grid on

subplot(4,1,3);
plot(ty,Lpol*dpy/1000); % dicharge chloride in tons/year
ylabel('Ton/y');  title('Total polder Cl^- discharge');
grid on

subplot(4,1,4);
plot(ty,cpol*1000); % discharged chloride concentration in mg/L
ylabel('mg/L');  title('Average polder discharge Cl concentration');
grid on
xlabel('time [y]');

%% Concentration and flux and concentration distribution over time
% Show dynamic graph of discharge flux acros poler at ym=0
figure
subplot(2,1,1); hold on;  % disharge flux
h=line(xm(IX),qcross(it,:)); ylim=[0 max(qcross(:))]; set(gca,'ylim',ylim); grid on
xlabel('x [m]'); ylabel('seepage [mm/d]');
for it=1:Nt,     set(h,'ydata',qcross(it,:)); set(gca,'ylim',ylim); drawnow; end

subplot(2,1,2); hold on  % discharge concentration
h=line(xm(IX),ccross(it,:)); ylim=[0 max(ccross(:))]; set(gca,'ylim',ylim); grid on
xlabel('x [m]'); ylabel('Cl^- conc g/L');
for it=1:Nt,     set(h,'ydata',ccross(it,:)); set(gca,'ylim',ylim); drawnow; end


%% Contouring the cross section throuth ym=0 over time (movie)
figure
cvals=0.5:1:9.5; % contour values in g/m3.

% setup first contour to get handles to graphic objects
it=1;  
[c,h]=contourf(xm,squeeze(zm),squeeze(permute(C(it).values(iy,:,:),[3,2,1])),cvals); hold on
ht=title(''); xlabel('Distance from polder center [m]'); ylabel('Depth [m]');
colorbar; title('gCl^-/L');

% show over time, by replacing data values to those of next time step
for it=1:Nt
    set(h,'ZData',squeeze(permute(C(it).values(iy,:,:),[3,2,1])));
    set(ht,'String',sprintf('concentration contours t=%.0f y colors in gCl^-/L',ty(it)));
    drawnow; % update figure
    if rem(it,20)==0,
        fprintf('Pauzing ... press any key');
        pause;
        fprintf('\n');
    end
end
    

%% Threedimensional movie of development of salt distribution
% it's a bit complex showing chloride isosrface and contour slices yielding
% %D contourlines in layers.

figure; view(10,10); grid on; hold on

xlabel('x [m]'); ylabel('y [m]'); zlabel('z [m]');
set(gca,'xlim',xGr([1 end]),'ylim',yGr([end 1]),'zlim',zGr([end 1]))

set(gca,'DataAspectRatio',[10,10,1]);

Nt   =length(C);
cpatch=5; cvals=0.5:1:9.5; caxis([0 10]);
s=[]; for i=1:length(cvals), s=[s sprintf(' %g',cvals(i))]; end
F(Nt)=getframe(gcf);
for it=1:Nt
    if any(exist('h1','var')) && any(ishandle(h1)), delete(h1); end
    if any(exist('h2','var')) && any(ishandle(h2)), delete(h2); end
    values=C(it).values;
    h1=patch(isosurface(xm,ym,zm,values,cpatch));
    set(h1,'FaceColor','r','EdgeColor','none');
    h2=contourslice(xm,ym,zm,values,xm(1),ym(end),zm,cvals);
    title(sprintf('%s Patch=%.0f; Cont=%s, t=%2d y',basename,cpatch,s,round(C(it).time/365)));   
    F(it) = getframe(gcf);
end

%% Make avi file
avifname=[basename,'.avi'];
movie2avi(F,avifname,'Quality',100,'Compression','none');
fprintf('\n\nAvi file <<%s>> has been generated for future playback\n\n\n',avifname);

%% Show movie if desired
%h2 = figure;
%movie(h2,F)
