function visualize_v7(fin,step_in,scales,labelflag,contourflag,boreholes)
%visualize_v7(fin,step_in,scales,labelflag,contourflag,boreholes)
%Visualization tool for output from runs using network_percolation_v2 using
%a network produced by hewitt_grid_v2
%The code produces a multipanel figure with rows of panels showing, in that
%sequence, size of S_R conduits, S_K conduits, effective pressure N with an
%option to superimpose surface and bed elevation contours, cavity storage
%(if activated in network_percolation_v2 through the Kamb_storage flag),
%and connectivity of S_R and S_K conduits along network edges, with red
%dots showing moulin locations. The code assumes qfun as supplied with the
%code package was used to define water input.
%Each row of these panels corresponds to a snapshot in time. Below the
%snapshots panels, the code places a plot showing water input time series
%for all moulins with red marker lines indicating the snapshots.
%To ensure short run times and minimize memory use, this tool projects the
%network on a pixellated grid so the matlab routines pcolor and contour can
%be used. The Hewitt (2013, EPSL) network is a tiling of individual "cells"
%of the form
%     x --- x
%     | \ / |
%     |  x  |
%     | / \ |
%     x --- x
%where x denotes a node and lines denote network edges. visualize_v6
%renders quantities defined at nodes in the following way
%     x avh x 
%    avv x avv
%     x avh x
%Where x is the value at a node itself (either those at the corners or
%centers of the cell depicted above. avh and avv denote algebraic means
%taken over nodal values that are adjacent horizontally (avh) or vertically
%(avv) to the positions of the "av" symbol. For effective pressure, for
%instance, these would correspond to the averages Pe_{ij} computed on the
%corresponding horizontal or vertical network edges. In other words, values
%are shown at nodes and (to fill gaps between them) on horizontal and
%vertical network edges, averaged appropriately. The diagonal edges are not
%rendered for quantities defined at nodes.
%Quantities defined at nodes are rendered as
%   max -- max
%    |  \/  |
%    |  /\  |
%   max -- max
%where max denotes the maximum value of the required quantity taken over
%all edges that connect to the corresponding node (including edges outside
%of the unit cell shown): in other words, horizontal and vertical edges
%take up two pixels each, diagonal edges take up one pixel each, the
%corners of the unit cell are filled with the maximum value over adjacent
%edges, and the node at the center of the unit cell defined above is not
%rendered when plotting quantities defined on network edges
%
%Input arguments:
%   fin:    output structure from a run of network_percolation, assumes the
%           water supply rate function qfun has been used.
%   step_in:list of integer indices of time steps in the solution to be
%           rendered as snapshots
%   scales: scale parameters for displaying output; for use for instance
%           when the solution was rendered in dimensionless terms. scales
%           is a structure with fields
%               N, S_K, S_R, Phi_0, Q, x, t: effective pressures, conduit
%               sizes, hydraulic potential values, water supply rates,
%               positions in fin.parameters.grid and times will be
%               multiplied by these scales. SUpply a value of "1" for each
%               if you do not want your output to be scaled
%               max: a substructure with fields N, S_R, S_K, Q: this bounds
%               the colour axis /  plot axes above; lower bound is zero for
%               all. If no substructure "max" is supplied as a field, the
%               code will assess maximum values in the solution to
%               determine the bounds.
%               ii_start: optional starting index, sets t = 0 to be the iith
%               time stamp in fin.t; useful for cutting initial transients
%               from the display.
%           If scales is not supplied (fewer than three arguments or an
%           empty array [] is supplied), default values are scales
%           appropriate to the parameters structure in default_params.mat
%           (where all parameter values are dimensionless)
%   labelflag: ensures panels are labelled if set to true, labelling scheme
%           is a letter for each row of columns, with numbers
%           distinguishing columns, and a separate letter label for the
%           water input plot at the bottom of the figure. Assumes the
%           vector step_in does not have more than four entries
%   contourflag: if set to true, superimposes bed and surface contours on
%           the effective pressure plots. If set to true, fin.parameters
%           needs to have n_nodes-by-1 fields Base and Surface (which are
%           not required to run network_percolation_v5) giving bed and
%           surface elevations at all node positions in the network. If
%           these are not present, the code will create a surface and bed
%           elevation from the hydraulic potential fin.parameters.Phi_0
%           assuming a specific cross-glacier surface profile as assumed
%           for the parameter structure in default_params.mat
%   boreholes: list of node indices for which effective pressure is plotted
%           against in a separate figure (optional, no additional figure is
%           created if this argument is not supplied). This can be used to
%           create time series from synthetic "borehole arrays".


%determine problem type
Kamb_storage = fin.parameters.Kamb_storage;

%set up number of panels
plotparams.n_rows = length(step_in)+1;
if Kamb_storage
    plotparams.n_cols = 5;
else
    plotparams.n_cols = 4;
end
plotparams.row_offset = .35;
plotparams.row_offset2 = .075;
plotparams.row_offset3 = .025;
plotparams.col_offset = .1;
plotparams.col_tickoffset = .025;
plotparams.row_tickoffset = .025;
plotparams.row_tickoffset2 = .1;
plotparams.row_cbaroffset =-.9/plotparams.n_rows*.8;

%generate plots
figure, hold on
f0 = gcf;
set(f0,'Units','centimeters')
set(f0,'Position',[10 10 8 3+4*plotparams.n_rows])
fontsize = 10;


%set up split colourmap
n_colours = 128;
cmap = colormap([channelcolour(n_colours); pressurecolour(n_colours)]);

%scales
if nargin < 3 || isempty(scales)
    scales.N = 0.7634;
    scales.S_K = 0.0332;
    scales.S_R = 0.0332;
    scales.x = 0.4681;
    scales.t = 3.7790;
    scales.Phi_0 =0.7634;
    scales.Q = 0.0231;
    scales.max.N = 1.5000;
    scales.max.S_R = 8;
    scales.max.S_K = 0.0100;
    scales.max.Q = 4;
end
if nargin >= 3
    if isfield(scales,'max')
        S_R_scale = scales.max.S_R;
        S_K_scale = scales.max.S_K;
        N_scale = scales.max.N;
        if Kamb_storage
            V_K_scale = scales.max.V_K;  
        end
    else
        S_R_scale = scales.S_R*max(max(fin.S_R));
        S_K_scale = scales.S_K*max(max(fin.S_K));
        N_scale = scales.N*max(max(fin.N));
        if Kamb_storage
            V_K_scale = scales.V_K*max(max(fin.V_K));  
        end
    end
else
    S_R_scale = max(max(fin.S_R));
    S_K_scale = max(max(fin.S_K));
    N_scale = max(max(fin.N));
    if Kamb_storage
        V_K_scale = max(max(fin.V_K));  
    end
end

%start time
if nargin < 3 || ~isfield(scales,'ii_start')
    ii_start = 1;
else
    ii_start = scales.ii_start;
end

%labels
if nargin < 4
    labelflag = false;
end

%contours
if nargin < 5
    contourflag = true;
end
if contourflag
    contourint = 100;
end

%boreholes
if nargin >= 6 && isempty(boreholes)
    boreholes =6860:5:6885;
end


%extract parameters
Phi_0 = fin.parameters.Phi_0;   %reduced hydraulic potential rho_i s + (rho_w-rho_i) b at each node (n_nodes-by-1 vector) 
S_P_R = fin.parameters.S_P_R;   %percolation cut-off size for R conduit
S_P_K = fin.parameters.S_P_K;   %percolation cut-off size for K conduit
if nargin >= 3
    Phi_0 = Phi_0*scales.Phi_0;
    S_P_R = S_P_R*scales.S_R;
    S_P_K = S_P_K*scales.S_K;
end

%load geometry
X = fin.parameters.grid.y;
Y = -fin.parameters.grid.x;
if nargin >= 3
    X = X*scales.x;
    Y = Y*scales.x;
end
Y = Y-min(Y);

%load surface and bed, if present
if isfield(fin.parameters,'Base')
    Base = fin.parameters.Base;
    Surface = fin.parameters.Surface;
end

%extract moulin parameters
moulins = 1:length(fin.parameters.moulins.base);
moulins = moulins(fin.parameters.moulins.base>0);
base =  fin.parameters.moulins.base(moulins);
envelope = fin.parameters.moulins.envelope(moulins,:);
omega_envelope = fin.parameters.moulins.omega_envelope(moulins,:);
phase_envelope = fin.parameters.moulins.phase_envelope(moulins,:);
amplitude = fin.parameters.moulins.amplitude(moulins,:);
omega_moulin = fin.parameters.moulins.omega_moulin(moulins,:);
phase_moulin = fin.parameters.moulins.phase_moulin(moulins,:);
Q = zeros(length(moulins),length(fin.t));
for ii=1:length(fin.t)
    t=fin.t(ii); 
    Q(:,ii) = base + amplitude(:,1).*(1+envelope(:,1).*cos(omega_envelope(:,1)*t+phase_envelope(:,1))).*cos(omega_moulin(:,1)*t+phase_moulin(:,1)) ...
        + amplitude(:,2).*(1+envelope(:,2).*cos(omega_envelope(:,2)*t+phase_envelope(:,2))).*cos(omega_moulin(:,2)*t+phase_moulin(:,2)); 
end
t = fin.t-min(fin.t);
if nargin >= 3
    t = scales.t*t;
    Q = scales.Q*Q;
end
Xmoulin = X(moulins);
Ymoulin = Y(moulins);

if nargin >= 3 && isfield(scales,'max')
    Q_scale = scales.max.Q;
else
    Q_scale = 1.25*max(max(Q));
end

%labels
labelstring = {'a','b','c','d','e','f'};
labelx = min(X) + 0.05*(max(X)-min(X));
labely = min(Y) + 0.95*(max(Y)-min(Y));

%edges
up_node = fin.parameters.grid.up_node;
down_node = fin.parameters.grid.down_node;
X_up = X(up_node); X_down = X(down_node);
Y_up = Y(up_node); Y_down = Y(down_node);
X_s = (X_up+X_down)/2;
Y_s = (Y_up+Y_down)/2;
ii_s = 1:fin.parameters.grid.n_edges;
%node spacing
dxp = max(abs(X_up-X_down))/2;
dyp = max(abs(Y_up-Y_down))/2;
nxp = round((max(X)-min(X))/dxp)+1;
nyp = round((max(Y)-min(Y))/dyp)+1;
%distinguish diagonally from vertically / horizontally connected nodes
X_square = X; Y_square = Y;
ii_square  = 1:fin.parameters.grid.n_nodes;
nxsquare = (nxp+1)/2; nysquare=(nyp-1)/2;
nxcross = nxp-nxsquare; nycross=nyp-nysquare;
[Y_square,indsquare] = sort(Y_square); X_square = X_square(indsquare); ii_square = ii_square(indsquare);
indsquare2 =  kron(nxcross:nxcross:(nycross-1)*nxcross,ones(1,nxsquare))+kron(0:nxsquare:(nycross-2)*nxsquare,ones(1,nxsquare)) + repmat(1:nxsquare,1,nysquare);
X_cross = X_square; Y_cross = Y_square; ii_cross = ii_square;
X_cross(indsquare2) = []; Y_cross(indsquare2) = []; ii_cross(indsquare2) = [];
X_cross = reshape(X_cross,nxcross,[]); Y_cross = reshape(Y_cross,nxcross,[]); ii_cross = reshape(ii_cross,nxcross,[]);
X_square = X_square(indsquare2); Y_square = Y_square(indsquare2); ii_square = ii_square(indsquare2);
X_square = reshape(X_square,nxsquare,[]); Y_square = reshape(Y_square,nxsquare,[]); ii_square = reshape(ii_square,nxsquare,[]);
[X_square,indsquare3] = sort(X_square,1);  
indsquare3 = indsquare3+kron(0:length(X_square(1,:))-1,length(X_square(:,1))*ones(length(X_square(:,1)),1));
Y_square = Y_square(indsquare3); ii_square = ii_square(indsquare3);
[X_cross,indcross] = sort(X_cross,1); 
indcross = indcross+kron(0:length(X_cross(1,:))-1,length(X_cross(:,1))*ones(length(X_cross(:,1)),1));
Y_cross = Y_cross(indcross); ii_cross = ii_cross(indcross);
%distinguish vertical, horizontal, diagonal edges
X_vert = X_s(X_s==X_up);
Y_vert = Y_s(X_s==X_up);
ii_vert = ii_s(X_s==X_up);
X_hor = X_s(Y_s==Y_up);
Y_hor = Y_s(Y_s==Y_up);
ii_hor = ii_s(Y_s==Y_up);
%ii_hor(Y_s==Y_up) = true;
X_diag = X_s(X_s~=X_up&Y_s~=Y_up);
Y_diag = Y_s(X_s~=X_up&Y_s~=Y_up);
ii_diag = ii_s(X_s~=X_up&Y_s~=Y_up);
%ii_diag(X_s~=X_up&Y_s~=Y_up) = true;
%grid parameters
dxs = max(abs(X_up-X_down))/3;
dys = max(abs(Y_up-Y_down))/3;
nxs = (nxp-1)*3/2;
nys = (nyp-1)*3/2;
%rearrange indexing of vertical, horizontal and diagonal edges 
nvert = nxs/3;
nhor = nxs/3-1;
ndiag = nxs/3*2;
[Y_vert,indvert] = sort(Y_vert); X_vert = X_vert(indvert); ii_vert = ii_vert(indvert);
X_vert = reshape(X_vert,nvert,[]); Y_vert = reshape(Y_vert,nvert,[]); ii_vert = reshape(ii_vert,nvert,[]);
[X_vert,indvert2] = sort(X_vert,1); 
indvert2 = indvert2+kron(0:length(X_vert(1,:))-1,length(X_vert(:,1))*ones(length(X_vert(:,1)),1));
Y_vert = Y_vert(indvert2); ii_vert = ii_vert(indvert2);
[Y_hor,indhor] = sort(Y_hor); X_hor = X_hor(indhor); ii_hor = ii_hor(indhor);
X_hor = reshape(X_hor,nhor,[]); Y_hor = reshape(Y_hor,nhor,[]); ii_hor = reshape(ii_hor,nhor,[]);
[X_hor,indhor2] = sort(X_hor,1);
indhor2 = indhor2+kron(0:length(X_hor(1,:))-1,length(X_hor(:,1))*ones(length(X_hor(:,1)),1));
Y_hor = Y_hor(indhor2); ii_hor = ii_hor(indhor2);
[Y_diag,inddiag] = sort(Y_diag); X_diag = X_diag(inddiag); ii_diag = ii_diag(inddiag);
X_diag = reshape(X_diag,ndiag,[]); Y_diag = reshape(Y_diag,ndiag,[]); ii_diag = reshape(ii_diag,ndiag,[]);
[X_diag,inddiag2] = sort(X_diag,1);
inddiag2 = inddiag2+kron(0:length(X_diag(1,:))-1,length(X_diag(:,1))*ones(length(X_diag(:,1)),1));
Y_diag = Y_diag(inddiag2); ii_diag = ii_diag(inddiag2);

%create mesh
[Xp,Yp] = ndgrid(min(X)+(0:nxp-1)*dxp,min(Y)+(0:nyp-1)*dyp);
[Xs,Ys] = ndgrid(min(X)+(0:nxs-1)*dxs,min(Y)+(0:nys-1)*dys);

%handle array for panels
if Kamb_storage
    hp = zeros(length(step_in),5);
else
    hp = zeros(length(step_in),4);
end

%initialize data arrays
Np = nan(size(Xp));
Phi_0p = nan(size(Xp));
SRs = nan(size(Xs));
SKs = nan(size(Xs));
SPRs = nan(size(Xs));
SPKs = nan(size(Xs));
if Kamb_storage, VKs = nan(size(Xs)); end

%allocation arrays
jjvert11= 2+(0:nvert-1)*3;
jjvert12 = 1:length(SRs(1,:));
jjvert12(3:3:end) = [];
jjvert12(2:2:end) = [];

jjvert21= 2+(0:nvert-1)*3;
jjvert22 = 1:length(SRs(1,:));
jjvert22(3:3:end) = [];
jjvert22(1:2:end) = [];

jjhor12= 3+(0:nhor-1)*3;
jjhor11 = 1:length(SRs(:,1));
jjhor11(2:3:end) = [];
jjhor11(2:2:end) = [];
jjhor11(1) = [];

jjhor21 = 1:length(SRs(:,1));
jjhor21(2:3:end) = [];
jjhor21(1:2:end) = [];
jjhor21(end) = [];
jjhor22 = jjhor12;

jjdiag1 = 1:length(SRs(:,1));
jjdiag1(2:3:end) = [];
jjdiag2 = 1:length(SRs(1,:));
jjdiag2(3:3:end) = [];

%hole filling arrays - edges
jj_holes1 = 2+(0:nhor)*3;
jj_holes2 = 1:length(SRs(1,:));
jj_holes2(1:3:end) = [];
jj_holes2(1:2:end) = [];

jj_holes1_plus = jj_holes1+1;
jj_holes1_plus(jj_holes1_plus>length(SRs(:,1))) = jj_holes1(jj_holes1_plus>length(SRs(:,1)))-1;

jj_holes1_minus = jj_holes1-1;
jj_holes1_minus(jj_holes1_minus<1) = jj_holes1(jj_holes1_minus<1)+1;

jj_holes2_plus = jj_holes2+1;
jj_holes2_plus(jj_holes2_plus>length(SRs(1,:))) = jj_holes2(jj_holes2_plus>length(SRs(1,:)))-1;

jj_holes2_minus = jj_holes2-1;
jj_holes2_minus(jj_holes1_minus<1) = jj_holes2(jj_holes2_minus<1)+1;

%hole filling arrays - nodesfout.rho_i
jj_holes11 = 2:2:nxp;
jj_holes12 = 2:2:nyp;

jj_holes11_plus = jj_holes11+1;
jj_holes11_plus(jj_holes11_plus>length(Np(:,1))) = jj_holes11(jj_holes11_plus>length(Np(:,1)))-1;

jj_holes11_minus = jj_holes11-1;
jj_holes11_minus(jj_holes11_minus<1) = jj_holes11(jj_holes11_minus<1)+1;


jj_holes21 = 1:2:nxp;
jj_holes22 = 1:2:nyp;

hold on
jj_holes22_plus = jj_holes22+1;
jj_holes22_plus(jj_holes22_plus>length(Np(1,:))) = jj_holes22(jj_holes22_plus>length(Np(1,:)))-1;

jj_holes22_minus = jj_holes22-1;
jj_holes22_minus(jj_holes22_minus<1) = jj_holes2(jj_holes22_minus<1)+1;

%Static data allocation - edge values
SPRs(jjvert11,jjvert12) = S_P_R(ii_vert);
SPRs(jjvert21,jjvert22) = S_P_R(ii_vert);
SPRs(jjhor11,jjhor12) = S_P_R(ii_hor);
SPRs(jjdiag1,jjdiag2) = S_P_R(ii_diag);

SPKs(jjvert11,jjvert12) = S_P_K(ii_vert);
SPKs(jjvert21,jjvert22) = S_P_K(ii_vert);
SPKs(jjhor11,jjhor12) = S_P_K(ii_hor);
SPKs(jjdiag1,jjdiag2) = S_P_K(ii_diag);

%fill holes - static data
SPRs(jj_holes1,jj_holes2) = max(SPRs(jj_holes1_plus,jj_holes2),SPRs(jj_holes1_minus,jj_holes2));
SPRs(jj_holes1,jj_holes2) = max(SPRs(jj_holes1,jj_holes2),SPRs(jj_holes1,jj_holes2_plus));
SPRs(jj_holes1,jj_holes2) = max(SPRs(jj_holes1,jj_holes2),SPRs(jj_holes1,jj_holes2_minus));
SPRs(jj_holes1,jj_holes2) = max(SPRs(jj_holes1,jj_holes2),SPRs(jj_holes1_minus,jj_holes2_minus));
SPRs(jj_holes1,jj_holes2) = max(SPRs(jj_holes1,jj_holes2),SPRs(jj_holes1_minus,jj_holes2_plus));
SPRs(jj_holes1,jj_holes2) = max(SPRs(jj_holes1,jj_holes2),SPRs(jj_holes1_plus,jj_holes2_minus));
SPRs(jj_holes1,jj_holes2) = max(SPRs(jj_holes1,jj_holes2),SPRs(jj_holes1_plus,jj_holes2_plus));

SPKs(jj_holes1,jj_holes2) = max(SPKs(jj_holes1_plus,jj_holes2),SPKs(jj_holes1_minus,jj_holes2));
SPKs(jj_holes1,jj_holes2) = max(SPKs(jj_holes1,jj_holes2),SPKs(jj_holes1,jj_holes2_plus));
SPKs(jj_holes1,jj_holes2) = max(SPKs(jj_holes1,jj_holes2),SPKs(jj_holes1,jj_holes2_minus));
SPKs(jj_holes1,jj_holes2) = max(SPKs(jj_holes1,jj_holes2),SPKs(jj_holes1_minus,jj_holes2_minus));
SPKs(jj_holes1,jj_holes2) = max(SPKs(jj_holes1,jj_holes2),SPKs(jj_holes1_minus,jj_holes2_plus));
SPKs(jj_holes1,jj_holes2) = max(SPKs(jj_holes1,jj_holes2),SPKs(jj_holes1_plus,jj_holes2_minus));
SPKs(jj_holes1,jj_holes2) = max(SPKs(jj_holes1,jj_holes2),SPKs(jj_holes1_plus,jj_holes2_plus));

%static data allocation - node values
Phi_0p(1:2:end,2:2:end) = Phi_0(ii_square); Phi_0p(2:2:end,1:2:end) = Phi_0(ii_cross);

%fill holes - nodes
Phi_0p(jj_holes11,jj_holes12) = .5*(Phi_0p(jj_holes11_plus,jj_holes12)+Phi_0p(jj_holes11_minus,jj_holes12));
Phi_0p(jj_holes21,jj_holes22) = .5*(Phi_0p(jj_holes21,jj_holes22_minus)+Phi_0p(jj_holes21,jj_holes22_plus));


%Surface = 0.1700*Yp*1000;
if contourflag
    if ~isfield(fin.parameters,'Base')
        Surfacep = 10^6/9.8/1000*Phi_0p(1,:);
        Surfacep = repmat(Surfacep,length(Phi_0p(1,:)),1);
        Surfacep = Surfacep - 50*Xp.*(1-Xp).*Yp/.5^2/5;
        Basep = (10^6*Phi_0p/9.8/910 - Surfacep)*910/(1000-910);
    else
        Surfacep(1:2:end,2:2:end) = Surface(ii_square); Surfacep(2:2:end,1:2:end) = Surface(ii_cross);
        Surfacep(jj_holes11,jj_holes12) = .5*(Surface(jj_holes11_plus,jj_holes12)+Surface(jj_holes11_minus,jj_holes12));
        Surfacep(jj_holes21,jj_holes22) = .5*(Surface(jj_holes21,jj_holes22_minus)+Surface(jj_holes21,jj_holes22_plus));
        Basep(1:2:end,2:2:end) = Base(ii_square); Basep(2:2:end,1:2:end) = Base(ii_cross);
        Basep(jj_holes11,jj_holes12) = .5*(Base(jj_holes11_plus,jj_holes12)+Base(jj_holes11_minus,jj_holes12));
        Basep(jj_holes21,jj_holes22) = .5*(Base(jj_holes21,jj_holes22_minus)+Base(jj_holes21,jj_holes22_plus));
    end
end

%set limits
if nargin >= 3 && isfield(scales,'lims')
    plotparams.lims = lims;
else
    plotparams.lims.x = [min(Xs)-dxs max(Xs)+dxs]; plotparams.lims.y = [min(Ys)-dys max(Ys)+dys];
end
plotparams.equal = true;


%prepare figure
clf

for ii = 1:length(step_in) 
    
step = step_in(ii);

%load data
N = fin.N(:,step);
if nargin >= 3
    N = N*scales.N;
end
S_R = fin.S_R(:,step);
if nargin >= 3
    S_R = S_R*scales.S_R;
end
S_K = fin.S_K(:,step);
if nargin >= 3
    S_K = S_K*scales.S_K;
end
if Kamb_storage
    V_K = fin.V_K(:,step);
    if nargin >= 3
        V_K = V_K*scales.V_K;
    end
end

%allocate values
Np(1:2:end,2:2:end) = N(ii_square); Np(2:2:end,1:2:end) = N(ii_cross);
SRs(jjvert11,jjvert12) = S_R(ii_vert);
SKs(jjvert11,jjvert12) = S_K(ii_vert);
if Kamb_storage, VKs(jjvert11,jjvert12) = V_K(ii_vert); end
SRs(jjvert21,jjvert22) = S_R(ii_vert);
SKs(jjvert21,jjvert22) = S_K(ii_vert);
if Kamb_storage, VKs(jjvert21,jjvert22) = V_K(ii_vert); end
SRs(jjhor11,jjhor12) = S_R(ii_hor);
SKs(jjhor11,jjhor12) = S_K(ii_hor);
if Kamb_storage, VKs(jjhor11,jjhor12) = V_Kfout.rho_i(ii_hor); end
SRs(jjhor21,jjhor22) = S_R(ii_hor);
SKs(jjhor21,jjhor22) = S_K(ii_hor);
SPRs(jjhor21,jjhor22) = S_P_R(ii_hor);
SPKs(jjhor21,jjhor22) = S_P_K(ii_hor);
if Kamb_storage, VKs(jjhor21,jjhor22) = V_K(ii_hor); end
SRs(jjdiag1,jjdiag2) = S_R(ii_diag);
SKs(jjdiag1,jjdiag2) = S_K(ii_diag);
if Kamb_storage, VKs(jjdiag1,jjdiag2) = V_K(ii_diag); end

%fill holes - edge values
SRs(jj_holes1,jj_holes2) = max(SRs(jj_holes1_plus,jj_holes2),SRs(jj_holes1_minus,jj_holes2));
SRs(jj_holes1,jj_holes2) = max(SRs(jj_holes1,jj_holes2),SRs(jj_holes1,jj_holes2_plus));
SRs(jj_holes1,jj_holes2) = max(SRs(jj_holes1,jj_holes2),SRs(jj_holes1,jj_holes2_minus));
SRs(jj_holes1,jj_holes2) = max(SRs(jj_holes1,jj_holes2),SRs(jj_holes1_minus,jj_holes2_minus));
SRs(jj_holes1,jj_holes2) = max(SRs(jj_holes1,jj_holes2),SRs(jj_holes1_minus,jj_holes2_plus));
SRs(jj_holes1,jj_holes2) = max(SRs(jj_holes1,jj_holes2),SRs(jj_holes1_plus,jj_holes2_minus));
SRs(jj_holes1,jj_holes2) = max(SRs(jj_holes1,jj_holes2),SRs(jj_holes1_plus,jj_holes2_plus));

SKs(jj_holes1,jj_holes2) = max(SKs(jj_holes1_plus,jj_holes2),SKs(jj_holes1_minus,jj_holes2));
SKs(jj_holes1,jj_holes2) = max(SKs(jj_holes1,jj_holes2),SKs(jj_holes1,jj_holes2_plus));
SKs(jj_holes1,jj_holes2) = max(SKs(jj_holes1,jj_holes2),SKs(jj_holes1,jj_holes2_minus));
SKs(jj_holes1,jj_holes2) = max(SKs(jj_holes1,jj_holes2),SKs(jj_holes1_minus,jj_holes2_minus));
SKs(jj_holes1,jj_holes2) = max(SKs(jj_holes1,jj_holes2),SKs(jj_holes1_minus,jj_holes2_plus));
SKs(jj_holes1,jj_holes2) = max(SKs(jj_holes1,jj_holes2),SKs(jj_holes1_plus,jj_holes2_minus));
SKs(jj_holes1,jj_holes2) = max(SKs(jj_holes1,jj_holes2),SKs(jj_holes1_plus,jj_holes2_plus));

%fill holes - nodes
Np(jj_holes11,jj_holes12) = .5*(Np(jj_holes11_plus,jj_holes12)+Np(jj_holes11_minus,jj_holes12));
Np(jj_holes21,jj_holes22) = 0.5*(Np(jj_holes21,jj_holes22_minus)+Np(jj_holes21,jj_holes22_plus));

%Cap SKs as needed
SRs = min(SRs,S_R_scale*(1-1/n_colours));
SKs = min(SKs,S_K_scale*(1-1/n_colours));

%same for Np
Np = max(Np,-N_scale/n_colours);

hp(ii,1)=plotpanel_v2(plotparams.n_cols*(ii-1)+1,plotparams);
pcolor(Xs,Ys,SRs/S_R_scale); shading flat
colormap(cmap)
caxis([0 2])
if labelflag
    text(labelx,labely,strcat('(',labelstring{ii},num2str(1),')'));
end
set(hp(ii,1),'YTick',[0 1 2 3 4 5])
set(hp(ii,1),'XTick',[])
if ii == length(step_in)
    xlabel(hp(ii,1),'$x$ (km)','Interpreter','Latex')
set(hp(ii,1),'XTick',[0 .5 1])
end
ylabel(hp(ii,1),'$y$ (km)','Interpreter','Latex')
%axis(h0,'equal')

hp(ii,2)=plotpanel_v2(plotparams.n_cols*(ii-1)+2,plotparams);
pcolor(Xs,Ys,SKs/S_K_scale); shading flat
colormap(cmap)
caxis([0 2])
if labelflag
    text(labelx,labely,strcat('(',labelstring{ii},num2str(2),')'));
end
set(hp(ii,2),'YTick',[])
set(hp(ii,2),'XTick',[])
if ii == length(step_in)
    xlabel(hp(ii,2),'$x$ (km)','Interpreter','Latex')
set(hp(ii,2),'XTick',[0 .5 1])
end
%axis(h2,'equal')

hp(ii,3)=plotpanel_v2(plotparams.n_cols*(ii-1)+3,plotparams);
pcolor(Xp,Yp,1+Np/N_scale); shading flat
colormap(cmap)
caxis([0 2])
if labelflag
    text(labelx,labely,strcat('(',labelstring{ii},num2str(3),')'));
end
set(hp(ii,3),'YTick',[])
set(hp(ii,3),'XTick',[])
%height contours
hold on
if contourflag
   contour(Xp,Yp,Surfacep,contourint:contourint:900,'Color','k') 
   contour(Xp,Yp,Basep,contourint:contourint:850,'Color',[.5 .5 .5]) 
end
if labelflag
    text(labelx,labely,strcat('(',labelstring{ii},num2str(1),')'));
end
if ii == length(step_in)
    xlabel(hp(ii,3),'$x$ (km)','Interpreter','Latex')
    set(hp(ii,3),'XTick',[0 .5 1])
end
%axis(h3,'equal')

if Kamb_storage
    hp(ii,4)=plotpanel_v2(plotparams.n_cols*(ii-1)+5,plotparams);
else
    hp(ii,4)=plotpanel_v2(plotparams.n_cols*(ii-1)+4,plotparams);
end
Kappas = nan(size(SRs));
Kappas(SRs>SPRs|SKs>SPKs) = 1;
pcolor(Xs,Ys,Kappas); shading flat
colormap(cmap)
caxis([0 2])
hold on
for jj = 1:length(moulins)
    plot(Xmoulin(jj),Ymoulin(jj),'o','MarkerFaceColor','r','MarkerEdgeColor','none','MarkerSize',5*Q(jj,step_in(ii))/Q_scale)
end
if labelflag
    if Kamb_storage
        text(labelx,labely,strcat('(',labelstring{ii},num2str(5),')'));
    else
        text(labelx,labely,strcat('(',labelstring{ii},num2str(4),')'));
    end
end
set(hp(ii,4),'XTick',[])
set(hp(ii,4),'YTick',[])
if ii == length(step_in)
    xlabel(hp(ii,4),'$x$ (km)','Interpreter','Latex')
set(hp(ii,4),'XTick',[0 .5 1])
end
if Kamb_storage
    hp(ii,5)=plotpanel_v2(plotparams.n_cols*(ii-1)+4,plotparams);  
    pcolor(Xs,Ys,VKs/V_K_scale); shading flat
    colormap(cmap)
    caxis([0 2])
    text(labelx,labely,strcat('(',labelstring{ii},num2str(4),')'));
    set(hp(ii,5),'XTick',[])
    set(hp(ii,5),'YTick',[])
    if ii == length(step_in)
        set(hp(ii,5),'XTick',[0 .5 1])
        xlabel(hp(ii,5),'$x$ (km)','Interpreter','Latex')
    end
end

end

%colourbars
haux=plotpanel_v2(plotparams.n_cols*length(step_in)+1,plotparams);
colormap(cmap)
caxis([0 2])
h = colorbar('Southoutside');
set(h,'FontSize',fontsize)
set(h,'XLim',[0 .99])
if S_R_scale >= 2
    set(h,'XTick',[0 floor(S_R_scale/2)/S_R_scale]);
    set(h,'XTickLabel',{'0',num2str(floor(S_R_scale/2)),num2str(2*floor(S_R_scale/2))})
elseif S_R_scale >= 1
    set(h,'XTick',[0 floor(S_R_scale)/S_R_scale]);
    set(h,'XTickLabel',{'0',num2str(floor(S_R_scale))});
end
xlabel(h,'$S_R$ (m$^2$)','Interpreter','Latex');
axis(haux,'off');

haux=plotpanel_v2(plotparams.n_cols*length(step_in)+2,plotparams);
colormap(cmap)
caxis([0 2])
h = colorbar('Southoutside');
set(h,'FontSize',fontsize)
set(h,'XLim',[0 .99])
pwr = floor(log10(S_K_scale));
if pwr >= 0
    if S_K_scale >= 2
        set(h,'XTick',[0 floor(S_K_scale/2)/S_K_scale]);
        set(h,'XTickLabel',{'0',num2str(floor(S_K_scale/2))});
    elseif S_K_scale >= 1
        set(h,'XTick',[0 floor(S_K_scale)/S_K_scale]);
        set(h,'XTickLabel',{'0',num2str(floor(S_K_scale))});
    end
else
    if S_K_scale/10^pwr > 1
        set(h,'XTick',[0 10^pwr*floor(S_K_scale/2/10^pwr)/S_K_scale]);
        set(h,'XTickLabel',{'0',num2str(floor(S_K_scale/2))});
    else
        set(h,'XTick',[0 .5])
        set(h,'XTickLabel',{'0',num2str(10^pwr/2)})
    end
end
xlabel(h,'$S_K$ (m$^2$)','Interpreter','Latex')
axis(haux,'off');

haux=plotpanel_v2(plotparams.n_cols*length(step_in)+3,plotparams);
colormap(cmap)
caxis([0 2])
h = colorbar('Southoutside');
set(h,'XLim',[1+1/n_colours 2])
if N_scale >= 2
    set(h,'XTick',1+[1/n_colours floor(N_scale/2)/N_scale 2*floor(N_scale/2)/(N_scale)])
    set(h,'XTickLabel',{'0',num2str(floor(N_scale/2)),num2str(2*floor(N_scale/2))})
elseif N_scale >= 1
    set(h,'XTick',1 + [1/n_colours floor(N_scale)/N_scale 2*floor(N_scale)/N_scale])
    set(h,'XTickLabel',{'0',num2str(floor(N_scale)),num2str(2*floor(N_scale))})
end
set(h,'FontSize',fontsize)
xlabel(h,'$N$ (MPa)','Interpreter','Latex')
axis(haux,'off');

if Kamb_storage
    haux=plotpanel_v2(plotparams.n_cols*length(step_in)+4,plotparams);
    colormap(cmap)
    caxis([0 2])
    h = colorbar('Southoutside');
    set(h,'FontSize',fontsize)
    set(h,'YLim',[0 1])
    if V_K_scale >= 2
        set(h,'YTick',[0 floor(V_K_scale/2)/V_K_scale 2*floor(V_K_scale/2)/V_K_scale])
        set(h,'YTickLabel',{'0',num2str(floor(V_K_scale/2)),num2str(2*floor(V_K_scale/2))})
    elseif V_K_scale >= 1
        set(h,'YTick',[0 floor(V_K_scale)/V_K_scale 2*floor(V_K_scale)/V_K_scale])
        set(h,'YTickLabel',{'0',num2str(floor(V_K_scale)),num2str(2*floor(V_K_scale))})
    end
    set(h,'FontSize',fontsize)
    xlabel(h,'$V_K$ (m$^{3}$)','Interpreter','Latex')
    axis(haux,'off');
end

boxwidth = 1 - plotparams.col_offset;
plotwidth = boxwidth - plotparams.col_tickoffset;
boxheight2 = plotparams.row_offset - plotparams.row_offset2 - plotparams.row_offset3;
plotheight = boxheight2 - plotparams.row_tickoffset2;
haux = subplot('position',[plotparams.col_offset plotparams.row_offset2 plotwidth plotheight]);
t_min = t(ii_start);
plot(t-t_min,Q,'k'); box on
xlabel(haux,'$t$ (days)','Interpreter','Latex')
ylabel(haux,'$Q$ (m$^3$ s$^{-1}$)','Interpreter','Latex')
xlim([min(t) max(t)-t_min]);
ylim([0 Q_scale])
hold on
for ii=1:length(step_in)
        plot([t(step_in(ii)) t(step_in(ii))]-t_min,[0 Q_scale],'r','LineWidth',1)
end
if labelflag
    text(min(t)+0.025*(max(t)-min(t)),0.9*Q_scale,strcat('(',labelstring{length(step_in)+1},')'));
end

if nargin >= 6
    figure, hold on
    for ii = 1:length(boreholes)
        plot(t-t_min,N_scale*fin.N(boreholes(ii),:),'k')
    end
    xlabel('$t$ (days)','Interpreter','Latex')
    ylabel('$N$ (MPa)','Interpreter','Latex')
    xlim([min(t) max(t)-t_min]);
    ylim([0 N_scale])
    box on
end

end

function z_out = channelcolour(ncolours)

if nargin == 0, ncolours = 64; end

shade = linspace(0,1,ncolours).';
shade = shade.^(1/3);

dark_min = 0;
dark_max = 2.5;
saturate = 0.7;

beta = -(dark_max-dark_min)/2;
delta = dark_min/2;
alpha = saturate + delta - 1 - beta;

blue = ones(size(shade));
green = 1-delta - alpha *shade.^2;
red = 1-delta + beta*shade + alpha*shade.^2;

z_out = [red green blue];
end

function z_out = pressurecolour(ncolours)

if nargin == 0, ncolours = 64; end

shade = linspace(0,1,ncolours).';
shade = 1.5*shade.^(5/2);

dark_min = 0;
dark_max = 2.5;
saturate = 0.7;

blue = zeros(size(shade));
green = ones(size(shade));
red = ones(size(shade));

blue(shade < 0.5) = 0.5-shade(shade < 0.5);
red(shade < 0.5) = 0.5+shade(shade < 0.5);
green(shade > 0.5) = 2*(1-shade(shade > 0.5));
green(shade > 1) = 0;
red(shade > 1) = 2-shade(shade > 1);

z_out = [red green blue];
z_out(1,:) = [0 0 1];

end

function hout = plotpanel_v2(ii,parameters,labels,lims,fontsize)
%hout = plotpanel_v2(ii,parameters,labels,lims,fontsize)
%Set up subplots with minimal whitespace. Input variables are
%   parameters: structure, containing the following optional fields
%       n_rows:     number of rows in grid of panels (default 1)
%       n_cols:  number of columns in grid of panels (default 3)
%       col_tickoffset: scaled (to total width) space between columns of
%                   panels (default 0.025)
%       row_tickoffset: scaled (to total height) space between rows of
%                   panels (default 0.025)
%       col_offset: scaled (to total width) space to the left of the columns to
%                   accommodate labels (default 0.1)
%       row_offset: scaled (to total height) space below the tows to
%                   accommodate labels (default 0.1)
%       row_cbaroffset:   additional space for colorbar in lowest row of
%                   panels (default 0)
%       col_cbaroffset:   leaves additional space for colorbar in right-hand column of
%                   panels (default 0)
%       box:        boolean, 'true' forces 'box on' (default false)
%       equal:      boolean, 'true' sets 'axis equal' (default false)
%       hold:       boolean, 'true' sets 'hold on' (default false)
%       log:        creats logscale on both axes
%       semilogx:   creates logscale on x axis
%       semilogy:   creates logscale on y axis
%   labels: structure containing the following optional fields
%       xlab:       x-axis label, displayed only if panel specified through
%                   ii is at the bottom of a column (default '')
%       ylab:       y-axis label, displayed only if panel specified through
%                   ii is at the left-hand end of a row (default '')
%       panel_lab:  panel label, placed top left hand corner (default '')
%   lims: structure containing the following optional fields
%       x:          one-by-two-vector defining x limits
%       y:          one-by-two-vector defining y limits
%   fontsize:       Font size to be used (in pts)
%   ii:   identifies which panel to generate, counting across from top left
%The output variable is
%   hout:           axes handle for axes created

if nargin < 2, parameters.dummy = []; end
if isfield(parameters,'n_rows')
    n_rows = parameters.n_rows;
else
    n_rows = 1;
end
if isfield(parameters,'n_cols')
    n_cols = parameters.n_cols;
else
    n_cols = 3;
end
if isfield(parameters,'col_offset')
    col_offset = parameters.col_offset;
else
    col_offset = .1;
end
if isfield(parameters,'row_offset')
    row_offset = parameters.row_offset;
else
    row_offset = .1;
end
if isfield(parameters,'col_tickoffset')
    col_tickoffset = parameters.col_tickoffset;
else
    col_tickoffset = .025;
end
if isfield(parameters,'row_tickoffset')
    row_tickoffset = parameters.row_tickoffset;
else
    row_tickoffset = .025;
end
if isfield(parameters,'col_cbaroffset')
    col_cbaroffset = parameters.col_cbaroffset;
else
    col_cbaroffset =  0;
end
if isfield(parameters,'row_cbaroffset')
    row_cbaroffset = parameters.row_cbaroffset;
else
    row_cbaroffset = 0;
end
if isfield(parameters,'box')
    boxflag = parameters.box;
else
    boxflag = false;
end
if isfield(parameters,'equal')
    equalflag = parameters.equal;
else
    equalflag = false;
end
if isfield(parameters,'hold')
    holdflag = parameters.hold;
else
    holdflag = false;
end
if isfield(parameters,'log')
    logflag = parameters.log;
else
    logflag = false;
end
if isfield(parameters,'semilogx')
    semilogxflag = parameters.semilogx;
else
    semilogxflag = false;
end
if isfield(parameters,'semilogy')
    semilogyflag = parameters.semilogy;
else
    semilogyflag = false;
end

boxwidth = (1-col_offset-col_cbaroffset)/n_cols;
boxheight = (1-row_offset-row_cbaroffset)/n_rows;

plotwidth = boxwidth-col_tickoffset;
plotheight = boxheight-row_tickoffset;

col = 1+mod(ii-1,n_cols);
row = n_rows-floor((ii-1)/n_cols);

if row<0 && col<0, error('plot outside figure'), end

if col < n_cols && row > 1
    hout = subplot('position',[col_offset+(col-1)*boxwidth row_offset+(row-1)*boxheight+row_cbaroffset plotwidth plotheight]);
elseif row > 1
    hout = subplot('position',[col_offset+(col-1)*boxwidth row_offset+(row-1)*boxheight+row_cbaroffset plotwidth+col_cbaroffset plotheight]);
elseif col < n_cols
    hout = subplot('position',[col_offset+(col-1)*boxwidth row_offset+(row-1)*boxheight plotwidth plotheight+row_cbaroffset]);
else
    hout = subplot('position',[col_offset+(col-1)*boxwidth row_offset+(row-1)*boxheight plotwidth+col_cbaroffset plotheight+row_cbaroffset]);
end
    
    
%set plot type
if logflag
    loglog(NaN,NaN)
elseif semilogxflag
    semilogx(NaN,NaN)
elseif semilogyflag
    semilogy(NaN,NaN)
end

%Font size
if nargin >= 5;
    set(gca,'FontSize',fontsize)
end

%finish labelling etc
if boxflag, box on, end
if equalflag, axis equal, end
if holdflag, hold on, end

if nargin < 3
    xlab  = ''; ylab = ''; panel_lab = '';
else
    if isfield(labels,'xlab')
        xlab = labels.xlab;
    else
        xlab = '';
    end
    if isfield(labels,'ylab')
        ylab = labels.ylab;
    else
        ylab = '';
    end
    if isfield(labels,'panel_lab')
        panel_lab = labels.panel_lab;
    else
        panel_lab = '';
    end
end
if nargin >= 4
    if ~isfield(lims,'x')
        lims.x(1) = -1; lims.x(2) = 1;
    end
    if ~isfield(lims,'y')
        lims.y(1) = -1; lims.y(2) = 1;
    end
end



if col==1
    ylabel(ylab)
else
    set(hout,'YTickLabel',{})
end
if row==1
    xlabel(xlab)
else
    set(hout,'XTickLabel',{})
end
if exist('lims')
    xlim(lims.x);
    ylim(lims.y);
    text(lims.x(1)+.1*(lims.y(2)-lims.y(1)),lims.y(2)-.1*(lims.y(2)-lims.y(1)),panel_lab,'FontSize',fontsize)
end

end