close all

mf_settings

load('name.mat') % loads name.mat, which only contains the variable "basename"                                
load(basename);  % having retrieved the baasename value, we can load
                 
load('adapt_settings.mat') % Load settings added by Yannick

H=readDat([basename,'','.hds']);  % use readDAT to read the heads  output file
B=readBud([basename,'','.bgt']);  B=maskHC(B,IBOUND);


[NROW,NCOL,NLAY]=size(IBOUND); % get the size of the model

[xGr,yGr,xm,ym,DELX,DELY,NX,NY]=modelsize(xGr,yGr); % xGr and yGr were saved in mf_adapt
% Read unformatted budget file and mask noflow cells if they exist

if strcmpi( FlowDirection, 'horizontal'),    
    A = (Z(1,1,1)-Z(1,1,2)) * (yGr(1)-yGr(end));
    dl = abs( xGr(end)-xGr(1) );
elseif strcmpi( FlowDirection, 'vertical')   
    A = (Z(1,1,1)-Z(1,1,2)) * (xGr(end)-xGr(1));
    dl = abs( yGr(end)-yGr(1) );
end
% hcellsize = dl / size(MAP, 2);

%%
% figure hydraulic head
figure('visible','off');
imagesc(H.values);
axis image;
h = colorbar;
ylabel(h, 'Hydraulic head (m)');
set( gca, 'XTick', [])
set( gca, 'YTick', [])
% title('Hydraulic head')



% Scale bar
scalebar
% bar_pixels = scale_bar_length / cellsize;
% for idx = 0.2*bar_pixels: 0.2*bar_pixels:bar_pixels,
%     rectangle('Position', [size(MAP,2)-bar_pixels - 50 + idx, 5, round(0.1*bar_pixels),10], 'FaceColor','w');
%     rectangle('Position', [size(MAP,2)-bar_pixels - 50 + idx + 0.1*bar_pixels, 5, round(0.1*bar_pixels),10], 'FaceColor','k');   
% end
% text( size(MAP,2)- 30 - round(0.6*bar_pixels), 30, barlabel, 'FontWeight', 'bold', 'Color','k');

% Load in geological map as base image
[BASE, cmapbase] = imread( baseimagefile ); 

figure('Visible','off');

if strcmpi( FlowDirection, 'horizontal'),
%   Plot geological map as base map
    imagecm(BASE, cmapbase);
    hold on
    
%   Create space above and below map  
    axis( [-9.5, size(BASE,2) + 10.5, 0.5, size(BASE,1) + 0.5 ]);
%   Draw rectangles to illustrate boundary conditions:    
    rectangle('Position', [ -9.5, 0.5, 10, size(BASE,1)], 'FaceColor', 'k' , 'EdgeColor', 'none');
    rectangle('Position', [ size(BASE,2)-0.5, 0.5, 11, size(BASE,1)], 'FaceColor', 'k' , 'EdgeColor', 'none');
%   Plot Streamlines  
    streamlinesh( B, 50, streamlinestyle);    
elseif strcmpi( FlowDirection, 'vertical')
%   Plot geological map as base map
    imagecm(BASE, cmapbase);
    hold on
    
%   Create space above and below map  
    axis( [0.5, size(BASE,2)+0.5, -9.5, size(BASE,1) + 10.5]);
%   Draw rectangles to illustrate boundary conditions:    
    rectangle('Position', [1, -9.5, size(BASE,2),10], 'FaceColor', 'k' , 'EdgeColor', 'none');
    rectangle('Position', [1, size(BASE,1), size(BASE,2),11], 'FaceColor', [ 0.5, 0.5, 0.5] , 'EdgeColor', 'none');
    rectangle('Position', [min( inflowrange), size(BASE,1), max(inflowrange) - min( inflowrange) ,11], 'FaceColor','k', 'EdgeColor', 'none');
%   Plot Streamlines  
    streamlinesv( B, 50, streamlinestyle);    
end

set( gca, 'XTick', [])
set( gca, 'YTick', [])

% xlabel( sprintf( '%1.1f m', dl) );
% ylabel( sprintf( '%1.1f m', A) );
scalebar






%% In and outflow through constant head cells
dh = STRTHD(1, 1) - STRTHD(end, end);

CH=B.term{strmatch('CONSTANTHEAD',B.label)};
QCHin =sum(CH(CH(:)>0));
QCHout=sum(CH(CH(:)<0));

bulkconductivity = QCHin / ( A * dh/dl) ;
bulkpermeability = bulkconductivity / (2.75 .* 0.001 .* (1/3.28)) ;

mass_balance_error = abs( 100*( abs(QCHin)-abs(QCHout))/abs(QCHin) );

fprintf(['Total inflow  across constant head cells    =  %g m^3/d\n',...
         'Total outflow across constant head cells    =  %g m^3/d\n',...
         '                                               -----------------\n',...
         'mass balance error                          =  %g %% \n',...
         'Bulk permeability                           =  %g mD \n'],...
         QCHin,QCHout, mass_balance_error, bulkpermeability);



%% If mass balance error is bigger then 1, indicate it on the images

if mass_balance_error > 1,
    figs = findall(0, 'type', 'figure');
    for figidx = 1: length( figs),
        figure( figs(figidx) );
        ax = axis;
        plot( ax( [ 1, 2 ] ), ax( [ 3, 4 ]), 'r-', ax( [1,2] ), ax( [4, 3]), 'r-' );
        text( ax( 1 ), ( 0.9 * ax( 4 )) , [ 'mass balance error = ', num2str( mass_balance_error ) ] );
    end
    
%     return;
end
     


%% Analyze flow per lithological unit
% 
QX = B.term{strmatch('FLOWRIGHTFACE', B.label)};
QY = B.term{strmatch('FLOWFRONTFACE', B.label)};

clear QXperlithology QYperlithology

ctr = 1;
for idx = imagevaluesused,
    QXperlithology( ctr ) =  sumall( abs(QX( MAP == idx)));
    QXperlithologypercell( ctr ) =  sumall( abs(QX( MAP == idx))) / sumall( MAP == idx);
    ctr = ctr+1;
end

ctr = 1;
for idx = imagevaluesused,
    QYperlithology( ctr ) =  sumall( abs(QY( MAP == idx)));
    QYperlithologypercell( ctr) =  sumall( abs(QY( MAP == idx))) / sumall( MAP == idx);
    ctr = ctr + 1;
end

% calculate as percentage.
QXperlithology = 100 * QXperlithology / sum(QXperlithology);
QYperlithology = 100 * QYperlithology / sum(QYperlithology);

figure('Visible','off');
bar( [ QXperlithology; QYperlithology]')
set(gca, 'XTickLabel', legenditems)
title('Flow per lithological unit')
legend( 'Flow in x direction', 'Flow in y direction')
ylabel('Percentage of total flow');

figure('Visible','off');    
bar( [ QXperlithologypercell; QYperlithologypercell]')
set(gca, 'XTickLabel', legenditems)
title('Flow per cell per lithological unit')
legend( 'Flow in x direction', 'Flow in y direction')
ylabel('{m}^{3} / day');


%% Alteration streamlines map

if altavailable,
    
    figure('Visible','off');
    [ALTIM, cmapalt] = imread(altimagefile) ;
    imagecm( ALTIM, cmapalt);
    hold on
    clear ALTIM, cmapalt;
    set( gca, 'XTick', [])
    set( gca, 'YTick', [])
    axis image

    hold on


    if strcmpi( FlowDirection, 'horizontal'),        
    %   Create space above and below map  
        axis( [-9.5, size(BASE,2) + 10.5, 0.5, size(BASE,1) + 0.5 ]);
    %   Draw rectangles to illustrate boundary conditions:    
        rectangle('Position', [ -9.5, 0.5, 10, size(BASE,1)], 'FaceColor', 'k' , 'EdgeColor', 'none');
        rectangle('Position', [ size(BASE,2)-0.5, 0.5, 11, size(BASE,1)], 'FaceColor', 'k' , 'EdgeColor', 'none');
    %   Plot Streamlines  
        streamlinesh( B, 50, streamlinestyle);    
    elseif strcmpi( FlowDirection, 'vertical')
    %   Create space above and below map  
        axis( [0.5, size(BASE,2)+0.5, -9.5, size(BASE,1) + 10.5]);
    %   Draw rectangles to illustrate boundary conditions:    
        rectangle('Position', [1, -9.5, size(BASE,2),10], 'FaceColor', 'k' , 'EdgeColor', 'none');
        rectangle('Position', [1, size(BASE,1), size(BASE,2),11], 'FaceColor', [ 0.5, 0.5, 0.5] , 'EdgeColor', 'none');
        rectangle('Position', [min( inflowrange), size(BASE,1), max(inflowrange) - min( inflowrange) ,11], 'FaceColor','k', 'EdgeColor', 'none');
    %   Plot Streamlines  
        streamlinesv( B, 50, streamlinestyle);    
    end    
    
    scalebar

    % Scale bar
%     for idx = 0.2*bar_pixels: 0.2*bar_pixels:bar_pixels,
%         rectangle('Position', [size(MAP,2)-bar_pixels - 50 + idx, 15, round(0.1*bar_pixels),10], 'FaceColor','w');
%         rectangle('Position', [size(MAP,2)-bar_pixels - 50 + idx + 0.1*bar_pixels, 15, round(0.1*bar_pixels),10], 'FaceColor','k');
% 
%     end
%     text( size(MAP,2)- 30 - round(0.55*bar_pixels), 32, barlabel, 'FontWeight', 'bold', 'Color','k');    
    
end

%%
screensize = get(0,'ScreenSize');
op = get( 1, 'OuterPosition');

if altavailable, 
    figs = 5;
else
    figs = 4;
end

for idx = 1 : figs,     
    set( idx, 'OuterPosition', screensize);
    set( idx, 'PaperPositionMode', 'auto');
end

print(1 ,'-depsc2', 'Head.eps');
print(2 ,'-depsc', 'Streamlines.eps');
print(3 ,'-depsc2', 'FlowPerUnit.eps');
print(4 ,'-depsc2', 'FlowPerUnitpercell.eps');
if altavailable,
    print(5 ,'-depsc', 'AlterationStreamlines.eps');
end

system( 'convert -density 300 Head.eps Head_eps.png');
system( 'convert -density 300 Streamlines.eps Streamlines_eps.png');
system( 'convert -density 300 FlowPerUnit.eps FlowPerUnit_eps.png');
system( 'convert -density 300 FlowPerUnitpercell.eps FlowPerUnitpercell_eps.png');
if altavailable,
    system( 'convert -density 300 AlterationStreamlines.eps AlterationStreamlines_eps.png');
end

save( ['mf_analyze_results_', datestr(now, 'yyyy_mm_dd HH_MM_SS')] );
save( ['bulkperm_', datestr(now, 'yyyy_mm_dd HH_MM_SS'), '.txt' ], 'bulkpermeability', '-ascii');
    
% saf
clear [op, screensize ]