function [g23orp] = eval_g23orp(~,g,~,model,gauss)
% function [g23orp] = eval_g23orp(x,g,v,model,whichgrid)
% computes d/drho(g2*g3/rho) = ddrho(gg2*g3*rho./H^2)
% = d/drho(gg2).*rho.*g3./H^2 
%    + d./drho(g3).*rho.*gg2./H^2 
%    + gg2.*g3./H.^2
%    - rho.*gg2.*g3./H.^2.*d/drho(H^2)
% #codegen


g3 = eval_g3([],g,[],model,gauss);
g3p = eval_g3p([],g,[],model,gauss);
gg2 = eval_gg2([],g,[],model,gauss);
H  = eval_H([],g,[],model,gauss);
H2 = (H.*H);

if ~gauss
    rho = model.rgrid.rho;
    gg2p = model.geom.LamGp*g(model.geom.ind_gg2);
    H2p = 2*H.*(model.geom.LamGp*g(model.geom.ind_H));
    %dH2/drho  = 2*H*dH/drho
else
    rho = model.rgrid.rhogauss;
    gg2p = model.geom.LamGpgauss*g(model.geom.ind_gg2);
    H2p = 2*H.*(model.geom.LamGpgauss*g(model.geom.ind_H));
end

iggpsV2 = 1./H2;
g23orp = (rho.*gg2.*g3p + rho.*g3.*gg2p + gg2.*g3 ...
  - rho.*gg2.*g3.*H2p.*iggpsV2).*iggpsV2;

debug_plot = false;
if debug_plot
eval_debugging_plot(mfilename,'rhotor',[],model,gauss, g3, 'g3', H, 'H', gg2, 'gg2', g3p, 'g3p', ...
    rho.*gg2.*g3p, 'rho.*gg2.*g3p', rho.*g3.*gg2p, 'rho.*g3.*gg2p', gg2.*g3, 'gg2.*g3',rho.*gg2.*g3.*H2p.*iggpsV2,  'rho.*gg2.*g3.*H2p.*iggpsV2');
end



return