function [g0] = eval_g0(~,g,~,model,gauss)
% function [g0] = eval_g0(x,g,v,model,whichgrid)
% g0 = <|nabla V|>
% #codegen

% g0 = gg0 / H

H = eval_H([],g,[],model,gauss); % (dpsi_dV)^2

if ~gauss
    gg0 = model.geom.LamG*g(model.geom.ind_gg0,:);
else
    gg0 = (model.geom.LamGgauss*g(model.geom.ind_gg0,:));
end

diverted = (H(end,:)==0);
g0 = gg0./H;
g0(end,diverted) = inf;

return