function ftrap = calc_ftrap(epsilon,varargin)
% trapped particle fraction expression depending on delta
% Equation by O. Sauter

if nargin>0
    delta = varargin{1};
else
    delta = zeros(size(epsilon)); % zero delta by default
end

assert(all(size(epsilon)==size(delta)),'epsilon, delta must have the same size');

epseff = 0.67 * (1. - 1.4*abs(delta).*delta).*epsilon;
aa=(1. - epsilon)./(1. + epsilon);

assert(all(aa>=0)    ,'(1-epsilon)/(1+epsilon) < 0')
assert(all(epseff>=0),'epseff < 0')

ftrap = 1. - sqrt(aa).*(1. - epseff)./(1. + 2.*sqrt(epseff));
end