function varargout = ntm_chie(varargin) 
%#codegen

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin==5;
    stap = varargin{1};
    geop = varargin{2};
    u = varargin{3};
    model  = varargin{4};
    params = varargin{5};
else
    error('must call with 6 inputs');
end

%% Compute output depending on model
ntmfact = ones(model.rgrid.nrhogauss,1); % init
iotagauss = stap.iota;
rhogauss = model.rgrid.rhogauss;

% Loop over all islands
for intm=1:model.ntm.n_ntm
    % Get mode numbers
    m = params.ntm.m(intm);
    n = params.ntm.n(intm);
    % Get width from the state or from the prescribed variables v
    w = stap.wntms(intm);
    % If an island has is present take into account effect on chie
    if (w > 0) && any(iotagauss>=n/m) && any(iotagauss<=n/m) % ensure NTM in plasma
        for imn = model.rgrid.nrhogauss:-1:1; % search from out to in
            if iotagauss(imn)>n/m; irho=imn; break; end
        end
        % interpolate to find rmn
        rhomn = rhogauss(irho) + ...
        (n/m-iotagauss(irho))*diff(rhogauss(irho+(0:1)))./diff(iotagauss(irho+(0:1)));
        % Use Chiefactor as used by S. Nowak for ITM-TF
        wn=w/(model.equi.epsilon*model.equi.R0); % Island width normalization based on minor radius
        ntmfact = ntmfact.*(1+params.ntm.chiefact_A(intm)*exp(-4*(rhogauss-rhomn).^2/(params.ntm.chiefact_w(intm)*wn).^2));
    end
end

% CHIE EFFECT
dntmfact_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);

varargout{1} = ntmfact;
varargout{2} = dntmfact_dx;
return