function [iotaprime,diotaprime_dx] = eval_diotadrho(x,g,~,model,gaussgrid)
% function [iota,diota_dx] = eval_diotadrho(x,g,model,whichgrid)
if ~gaussgrid;
    gauss = false;
    ngrid = model.rgrid.nrho;
else
    gauss = true;
    ngrid = model.rgrid.nrhogauss;
end

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

nt = size(g,2);
iotaprime = zeros(ngrid,nt);
if nt>1 && der
    error('can not compute derivative for multiple time steps')
end

for it = 1:nt
    Phib = eval_Phib([],g(:,it),[],model,[]);
    if ~gauss
        % avoid singularity at rho=0 by using l'hopital's rule for
        % first rho point.
        diotaprime_dpsihat = (1./(2*Phib)) * model.psi.Lamppr;
        iotaprime(:,it) = diotaprime_dpsihat * x(model.psi.xind,it);
        if der
            diotaprime_dx(:,model.psi.xind) = diotaprime_dpsihat;
        end
    else
        % all points at rho>0, no need for tricks
        diotaprime_dpsi = 1./(2*Phib) * model.psi.Lampprgauss;
        iotaprime(:,it) = diotaprime_dpsi*x(model.psi.xind,it);
        if der
            diotaprime_dx(:,model.psi.xind) = diotaprime_dpsi;
        end
    end
end


return
