function [iota,diota_dx] = eval_iota(x,g,~,model,gauss)
% function [iota,diota_dx] = eval_iota(x,g,v,model,whichgrid)
% #codegen

if ~gauss;
    ngrid = numel(model.rgrid.rho);
else
    ngrid = numel(model.rgrid.rhogauss);
end

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

nt = size(g,2);
iota = 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.
        diota_dpsihat = (1./(2*Phib)) * model.psi.Lampor;
        iota(:,it) = diota_dpsihat * x(model.psi.xind,it);
        
        if der
            diota_dx(:,model.psi.xind) = diota_dpsihat;
        end
    else
        % all points at rho>0, no need for tricks
        diota_dpsi = 1./(2*Phib) * model.psi.Lamporgauss;
        iota(:,it) = diota_dpsi*x(model.psi.xind,it);
        if der
            diota_dx(:,model.psi.xind) = diota_dpsi;
        end
    end
    
end


return
