function fout = network_timestep_Jacobian_v5(v_in,parameters)
%Jacobian of network_timestep, see the latter for documentation

%unpack parameters
n_nodes = parameters.grid.n_nodes;   %number of nodes
n_edges = parameters.grid.n_edges;   %number of edges
up_node = parameters.grid.up_node;   %list (n_edges-by-1 vector) of 'upstream' node for each edge
down_node = parameters.grid.down_node;   %list (n_edges-by-1 vector) of 'downstream' nodes for each edge
bdy_nodes = parameters.grid.bdy_nodes;   %list (n_bdy_nodes-by-1 vector) of nodes at domain boundary

dt = parameters.dt;         %time step

L = parameters.L;           %list (n_edges-by-1 vector) of edge lengths
n_c = parameters.n_c;       %list (n_edges-by-1 vector) of number of conduits along edge
Phi_0 = parameters.Phi_0;   %reduced hydraulic potential rho_i s + (rho_w-rho_i) b at each node (n_nodes-by-1 vector) 
N_bdy_nodes = parameters.N_bdy_nodes;   %Dirichlet conditions (n_bdy_nodes-by-1 vector) for effective pressure at domain boundary

c_1 = parameters.c_1;       %relates opening rate to Q*Psi
c_2 = parameters.c_2;       %relates closure rate to S*N*n
c_3 = parameters.c_3;       %relates discharge Q to S^alpha*Psi^beta
alpha = parameters.alpha;   %exponent in dependence of discharge Q on cross-section S
beta = parameters.beta;     %exponent in dependence of discharge Q on hydraulic gradient Psi, same convention as in Schoof et al 2012 / Hewitt et al 2012 so 'frozen-time' problenm for N only becomes beta-Laplacian (i.e. p-Laplacian with p=beta)
epsilon = parameters.epsilon;   %regularization parameter for hydraulic gradient
n_Glen = parameters.n_Glen; %Glen's law exponent
uh = parameters.uh;         %list (n_edges-by-1 vector) of cavity opening rates for each network edge
S_0_R = parameters.S_0_R;       %cut-off size for cavity opening
S_0_K  = parameters.S_0_K;  %cut-off size for Kamb cavity opening
V_0_K  = parameters.V_0_K;  %cut-off size for storage cavity opening
S_P_R = parameters.S_P_R;   %percolation cut-off size for R conduit
S_P_K = parameters.S_P_K;   %percolation cut-off size for K conduit
epsilon_P_R = parameters.epsilon_P_R;   %percolation cut-off regularizer for R-conduits;
epsilon_P_K = parameters.epsilon_P_K;   %percolation cut-off regularizer for K-conduits;
T = parameters.T;           %list (n_edges-by-1 vector) of tortuosities for 'cavities' along each network edge
nu = parameters.nu;         %list (n_edges-by-1 vector) of step size ratios for 'channel' along each network edge
gamma_S = parameters.gamma_S;   %contribution of conduit evolution to mass balance
gamma_store = parameters.gamma_store;   %list (n_nodes-by-1 vector) of 'ealstic' water storage capacities at each node (i.e. storage is affine in N with coefficient gamma_store)
k_leak = parameters.k_leak; %scalar leakage permeability
r = parameters.r;           %ice-to-water density ratio

Kamb_storage = parameters.Kamb_storage;     %Boolean that determines whether Kamb-style storage cavities are included at each node
if Kamb_storage
    uh_K = parameters.uh_K; %list (n_nodes-bv-1 vector) of cavity opening rates for storage cavities at each node
end

%unpack input variable v_in
if (Kamb_storage && length(v_in) ~= 2*(n_edges+n_nodes)) || (~Kamb_storage && length(v_in) ~= 2*n_edges+n_nodes), error('input size incompatible with network geometry'), end
                                                        %check input size
S_R = v_in(1:n_edges);
S_K = v_in(n_edges+1:2*n_edges);
N = v_in(2*n_edges+1:2*n_edges+n_nodes);
N(bdy_nodes) = N_bdy_nodes;
if Kamb_storage
    V_K = v_in(2*n_edges+n_nodes+1:2*n_edges+2*n_nodes);  %optional Kamb-style storage cavity at node
end

%compute hydraulic gradient and derivatives with respect to N's along each network edge
Psi = (Phi_0(up_node) - N(up_node) - Phi_0(down_node) + N(down_node))./L;
dPsidN_up = -1./L;
dPsidN_down = 1./L;

%same for regularized absolute values of hydraulic gradient
Psi_abs_R = (Psi.^2+epsilon.^2).^(1/2);
Psi_abs_K = (Psi.^2./T.^2+epsilon.^2).^(1/2);
dPsi_abs_RdPsi = Psi./Psi_abs_R;
dPsi_abs_KdPsi = Psi./(T.^2.*Psi_abs_K);

%compute effective pressure and derivatives with respect to N's along each network edge
N_edge = (N(up_node)+N(down_node))/2;
dN_edgedN_up = 1/2;
dN_edgedN_down = 1/2;

%compute 'channel' flux and derivatives along each edge
Q_R = c_3*(((S_R-S_P_R).^2+epsilon_P_R.^2).^(1/2)-epsilon_P_R).^alpha.*Psi_abs_R.^(beta-2).*Psi;
Q_R(S_R<S_P_R)=0;
%Q_R = c_3*S_R.^alpha.*abs(Psi).^(beta-2).*Psi;
dQ_RdS_R = alpha*c_3*(((S_R-S_P_R).^2+epsilon_P_R.^2).^(1/2)-epsilon_P_R).^(alpha-1).*(S_R-S_P_R)./((S_R-S_P_R).^2+epsilon_P_R.^2).^(1/2)   ...
    .*Psi_abs_R.^(beta-2).*Psi;
dQ_RdS_R(S_R<S_P_R)=0;
dQ_RdPsi = (beta-2).*Q_R.*dPsi_abs_RdPsi./Psi_abs_R ...
    + c_3*(((S_R-S_P_R).^2+epsilon_P_R.^2).^(1/2)-epsilon_P_R).^alpha.*Psi_abs_R.^(beta-2);
dQ_RdPsi(S_R<S_P_R)=0;

%compute 'cavity' flux and derivatives along each edge
Q_K = c_3*(((S_K-S_P_K).^2+epsilon_P_K.^2).^(1/2)-epsilon_P_K).^alpha.*Psi_abs_K.^(beta-2).*Psi./T;
Q_K(S_K<S_P_K)=0;
%Q_K = c_3*S_K.^alpha.*abs(Psi./T).^(beta-2).*Psi./T;
dQ_KdS_K = alpha*c_3*(((S_K-S_P_K).^2+epsilon_P_K.^2).^(1/2)-epsilon_P_K).^(alpha-1).*(S_K-S_P_K)./((S_K-S_P_K).^2+epsilon_P_K.^2).^(1/2)   ...
    .*Psi_abs_K.^(beta-2).*Psi./T;
dQ_KdS_K(S_K<S_P_K)=0;
dQ_KdPsi = (beta-2).*Q_K.*dPsi_abs_KdPsi./Psi_abs_K ...
    + c_3*(((S_K-S_P_K).^2+epsilon_P_K.^2).^(1/2)-epsilon_P_K).^alpha.*Psi_abs_K.^(beta-2)./T;
dQ_KdPsi(S_K<S_P_K)=0;

%compute channel opening rate derivatives
dfout1dS_R = speye(n_edges)/dt + spdiags(-c_1*dQ_RdS_R.*Psi + nu.*uh./S_0_R + c_2*abs(N_edge).^(n_Glen-1).*N_edge,0,n_edges,n_edges);
dfout1dN = sparse(1:n_edges,up_node,-c_1*dQ_RdPsi.*Psi.*dPsidN_up-c_1*Q_R.*dPsidN_up+c_2*n_Glen*S_R.*abs(N_edge).^(n_Glen-1).*dN_edgedN_up,n_edges,n_nodes) ...
    +sparse(1:n_edges,down_node,-c_1*dQ_RdPsi.*Psi.*dPsidN_down-c_1*Q_R.*dPsidN_down+c_2*n_Glen*S_R.*abs(N_edge).^(n_Glen-1).*dN_edgedN_down,n_edges,n_nodes);
%correct for Dirichlet conditions
dfout1dN(:,bdy_nodes) = 0;

%compute cavity opening rate derivatives
dfout2dS_K = speye(n_edges)/dt + spdiags(uh./S_0_K + c_2*abs(N_edge).^(n_Glen-1).*N_edge,0,n_edges,n_edges);
dfout2dN = sparse(1:n_edges,up_node,c_2*n_Glen*S_K.*abs(N_edge).^(n_Glen-1).*dN_edgedN_up,n_edges,n_nodes) ...
    +sparse(1:n_edges,down_node,c_2*n_Glen*S_K.*abs(N_edge).^(n_Glen-1).*dN_edgedN_down,n_edges,n_nodes);
%correct for Dirichlet conditions
dfout2dN(:,bdy_nodes) = 0;

%compute mass conservation derivatives
dfout3dN = -spdiags(gamma_store,0,n_nodes,n_nodes)/dt + sparse(up_node,up_node,(dQ_RdPsi+(n_c-1).*dQ_KdPsi+k_leak).*dPsidN_up,n_nodes,n_nodes) ...
    +  sparse(up_node,down_node,(dQ_RdPsi+(n_c-1).*dQ_KdPsi+k_leak).*dPsidN_down,n_nodes,n_nodes) ...
    - sparse(down_node,up_node, (dQ_RdPsi+(n_c-1).*dQ_KdPsi+k_leak).*dPsidN_up,n_nodes,n_nodes) ...
    - sparse(down_node,down_node,(dQ_RdPsi+(n_c-1).*dQ_KdPsi+k_leak).*dPsidN_down,n_nodes,n_nodes)...
    - sparse(up_node,up_node,r/2*gamma_S(up_node).*L.*(dQ_RdPsi.*Psi+Q_R).*dPsidN_up,n_nodes,n_nodes)...
    - sparse(up_node,down_node,r/2*gamma_S(up_node).*L.*(dQ_RdPsi.*Psi+Q_R).*dPsidN_down,n_nodes,n_nodes)...
    - sparse(down_node,up_node,r/2*gamma_S(down_node).*L.*(dQ_RdPsi.*Psi+Q_R).*dPsidN_up,n_nodes,n_nodes)...
    - sparse(down_node,down_node,r/2*gamma_S(down_node).*L.*(dQ_RdPsi.*Psi+Q_R).*dPsidN_down,n_nodes,n_nodes);
dfout3dS_R = sparse(up_node,1:n_edges,dQ_RdS_R,n_nodes,n_edges) - sparse(down_node,1:n_edges,dQ_RdS_R,n_nodes,n_edges) ...
    + 1/(2*dt)*(sparse(up_node,1:n_edges,L.*gamma_S(up_node),n_nodes,n_edges) + sparse(down_node,1:n_edges,L.*gamma_S(down_node),n_nodes,n_edges))...
    - sparse(up_node,1:n_edges,r/2*gamma_S(up_node).*L.*dQ_RdS_R.*Psi,n_nodes,n_edges) - sparse(down_node,1:n_edges,r/2*gamma_S(down_node).*L.*dQ_RdS_R.*Psi,n_nodes,n_edges);
dfout3dS_K = sparse(up_node,1:n_edges,(n_c-1).*dQ_KdS_K,n_nodes,n_edges) - sparse(down_node,1:n_edges,(n_c-1).*dQ_KdS_K,n_nodes,n_edges) ...
    + 1/(2*dt)*(sparse(up_node,1:n_edges,L.*T.*(n_c-1).*gamma_S(up_node),n_nodes,n_edges) + sparse(down_node,1:n_edges,L.*T.*(n_c-1).*gamma_S(down_node),n_nodes,n_edges));
%correct for Dirichlet conditions
dfout3dN(:,bdy_nodes) = 0;
dfout3dN(bdy_nodes,:) = 0;
dfout3dN = dfout3dN + sparse(bdy_nodes,bdy_nodes,ones(length(bdy_nodes),1),n_nodes,n_nodes);
dfout3dS_R(bdy_nodes,:) = 0;
dfout3dS_K(bdy_nodes,:) = 0;

%correct for Kamb storage, compute Kamb storage change if necessary
if Kamb_storage
    dfout3dV_K = speye(n_nodes)/dt;
    dfout4dV_K = speye(n_nodes)/dt + spdiags(uh_K./V_0_K + c_3*abs(N).^(n_Glen-1).*N,0,n_nodes,n_nodes);
    dfout4dN = spdiags(c_3*n_Glen*V_K.*abs(N).^(n_Glen-1),0,n_nodes,n_nodes);
    %correct for Dirichlet conditions
    dfout3dV_K(bdy_nodes,:) = 0;
    dfout4dN(:,bdy_nodes) = 0;
end

%Assemble final Jacobian
if ~Kamb_storage
    fout = [dfout1dS_R sparse(n_edges,n_edges) dfout1dN; ...
        sparse(n_edges,n_edges) dfout2dS_K dfout2dN; ...
        dfout3dS_R dfout3dS_K dfout3dN];
else
    fout = [dfout1dS_R sparse(n_edges,n_edges) dfout1dN sparse(n_edges,n_nodes); ...
        sparse(n_edges,n_edges) dfout2dS_K dfout2dN sparse(n_edges,n_nodes); ...
        dfout3dS_R dfout3dS_K dfout3dN dfout3dV_K;...
        sparse(n_nodes,2*n_edges) dfout4dN dfout4dV_K];
end

fout = real(fout);