function [jdotB,djdotB_dx] = eval_jdotB(x,g,~,model,gaussgrid)
% function [jdotB,jdotB_dx] = eval_jdotB(x,g,[],model,whichgrid)
% #codegen

if ~gaussgrid;
    gauss = false;
    ngrid = model.rgrid.nrho;
else
    gauss = true;
    ngrid = model.rgrid.nrhogauss;
end

if nargout>1
    der = true; % compute derivative
    djdotB_dx = zeros(ngrid,numel(x)); % init
else
    der = false;
end

mu0 = pi*4e-7;
Phib = eval_Phib([],g,[],model,[]);
k = (2*pi)/(16*pi^3*mu0*Phib);

Vp   = eval_Vp([],g,[],model,gauss);
Vpp  = eval_Vpp([],g,[],model,gauss);
F    = eval_F([],g,[],model,gauss);
g23or = eval_g23or([],g,[],model,gauss); % g2*g3/rho
g23orp = eval_g23orp([],g,[],model,gauss); % d/drho(g2*g3/rho)

kk = k.*F.^2./Vp;

if ~gauss
  
    djdotB_dpsihat = ...
            bsxfun(@times,kk.*g23orp,model.psi.Lamp) + ...
            bsxfun(@times,kk.*g23or,model.psi.Lampp);
    % l'Hopital for 0/0
    djdotB_dpsihat(1,:) = 2*k.*F(1)^2 ./ Vpp(1).*g23orp(1) * model.psi.Lampp(1,:);
else
    % gauss grid
    djdotB_dpsihat = ...
            bsxfun(@times,kk.*g23orp,model.psi.Lampgauss) + ...
            bsxfun(@times,kk.*g23or,model.psi.Lamppgauss);
end

jdotB = djdotB_dpsihat*x(model.psi.xind);
if der
    djdotB_dx(:,model.psi.xind) = djdotB_dpsihat;
end

return
