function [g1] = eval_g1(~,g,~,model,gauss)
% function [g1] = eval_g1(x,g,v,model,whichgrid)
% g1 = <(nabla V)^2>
% #codegen

% g1 = gg1.*rho/ H2

H = eval_H([],g,[],model,gauss);

if ~gauss
    r = model.rgrid.rho;
    gg1 = model.geom.LamG*g(model.geom.ind_gg1,:);
else
    r = model.rgrid.rhogauss;
    gg1 = (model.geom.LamGgauss*g(model.geom.ind_gg1,:));
end

diverted = (H(end,:)==0);
g1 = bsxfun(@times,r,gg1)./(H.*H);
g1(end,diverted) = inf;

return