function Hess = Hess_DC_u(P,Prob)
% returns the Hessian of the objective function to SNOPT
% J = Sum of the mean of the neural excitations to the power of n_p
% Last modification: 04/10/2008 by Marko Ackermann

n_nod = Prob.user.n_nod;
%model = Prob.user.model; % spedifies model used
n_f =  Prob.user.n_f; % number of DoFs
n_m = Prob.user.n_m; % number of muscles
n_p = Prob.user.n_p; % power in the cost function
c1_J = Prob.user.c1_J; % scaling of cost function
F_isomax = Prob.user.F_isomax; % vector of CE maximal isometric forces
L_ceopt = Prob.user.L_ceopt; % vector of optimal muscle CE lengths

% Allocate memory for the sparse Hessian matrix
Hess = spalloc((2*n_f+3*n_m)*(n_nod-1)+1,(2*n_f+3*n_m)*(n_nod-1)+1,n_m*(n_nod-1));

% Precomputation to decrease computational time
I_V = n_p*(n_p-1)*([L_ceopt L_ceopt].*[F_isomax F_isomax])/((n_nod-1)*[L_ceopt L_ceopt]*[F_isomax F_isomax]');

for k = 1:n_nod-1
    % diagonal 'sub'-submatrix in the (2*n_f+3*n_m)x(2*n_f+3*n_m)-submatrix (k,k)
    Hess((k-1)*(2*n_f+3*n_m)+(2*n_f+2*n_m)+1:k*(2*n_f+3*n_m), (k-1)*(2*n_f+3*n_m)+(2*n_f+2*n_m)+1:k*(2*n_f+3*n_m)) = c1_J*diag(I_V'.*(P((k-1)*(2*n_f+3*n_m)+(2*n_f+2*n_m)+1:k*(2*n_f+3*n_m),1).^(n_p-2))); 
end

% Last column of the Hessian partial derivative with respect to half cycle
% duration - zero