function [g2p] = eval_g2p(~,g,~,model,gauss)
% function [g2p] = eval_g2p(x,g,v,model,whichgrid)
% d/drho(g2) = d/drho(<(nabla V)^2>)
% #codegen

gg2 = eval_gg2([],g,[],model,gauss);
H = eval_H([],g,[],model,gauss);
H2 = (H.*H);
% g2 = gg2.*r^2 / H2
% g2p = gg2p * r^2/H2 + 2*r*gg2/H2 - H2p * gg2*r^2/H2^2

if ~gauss
    r = model.rgrid.rho;
    gg2p = model.geom.LamGp*g(model.geom.ind_gg2,:);
    H2p = 2*H.*(model.geom.LamGp*g(model.geom.ind_H,:));
else
    r = model.rgrid.rhogauss;
    gg2p = model.geom.LamGpgauss*g(model.geom.ind_gg2,:);
    H2p = 2.*H.*(model.geom.LamGpgauss*g(model.geom.ind_H,:));
end

g2p = bsxfun(@times,r.^2,gg2p)./H2 ...
       + 2*r.*gg2./H2 ...
       - H2p.*(bsxfun(@times,r.^2,gg2)./H2.^2);

diverted = (H(end,:)==0);
g2p(end,diverted) = inf;

return