function [varargout] = eval_Ipl(x,g,~,model,gauss)
% function [Ipl,dIpl_dx] = eval_Ipl(x,g,v,model,whichgrid)
% #codegen
coder.extrinsic('warning')

if nargout>1
    der = true; % compute derivative
else
    der = false;
end

% 1/(8*pi^3*mu0) * 1/(2*Phib) * F*g2*g3/rho * dpsi/drho
mu0 = pi*4e-7;
Phib = eval_Phib([],g,[],model,[]);
F    = eval_F([],g,[],model,gauss);
g23or = eval_g23or([],g,[],model,gauss);
psip = eval_psip(x,[],[],model,gauss);

k = 1./(16*pi^3*mu0*Phib);
tt = bsxfun(@times,k,F.*g23or);

Ipl = tt .* psip;

if ~gauss && (psip(end) == 0)
    warning('Ip eval called with diverted plasma. Call FF if you see this')
    psipp = eval_psipp(x,[],[],model,gauss);
    g23orp = eval_g23orp([],g,[],model,gauss);
    Ipl(end) = k*F(end)*g23orp(end)./psipp(end);
else
    psipp = psip; % dummy
    g23orp = psip; % dummy
end

if der
     [~,dpsip_dx] = eval_psip(x,[],[],model,gauss);
     dIpl_dx = bsxfun(@times,tt,dpsip_dx);

     if ~gauss && (psip(end)==0)
         % special case for l'Hopital
         [~,dpsipp_dx] = eval_psipp(x,[],[],model,gauss);
         dIpl_dx(end,:) = -k*F(end)*g23orp(end)./(psipp(end).^2).*(dpsipp_dx(end,:));
     end
    varargout{2} = dIpl_dx;
end

varargout{1} = Ipl;

return
