function F = planckTaper2D(radius,epsilon,X,Y)

R = sqrt(X.^2 + Y.^2);

radius2 = radius*(1-2*epsilon);
F = (1./(exp(-(2*epsilon*radius+eps)*(1./(R-radius) + 1./(R-radius2))) + 1)).*(R<=radius-eps).*(R>radius2) + (R<=radius2);

% F(isnan(F)) = 1;
%   
% 
% F = ones(size(X));
% 
% eps2 = 1-2*eps;
% 
% Sp = sqrt(X.^2 + Y.^2);
% 
% ind = find(Sp>eps2*radius);
% 
% Sp2 = Sp;
% 
% Sp2(ind) = 1;
% 
% ind = find(Sp>=radius);
% 
% Sp2(ind) = 0;
% 
% F(ind) = 0;
% 
% ind = find(Sp <= eps2*radius);
% 
% Sp2(ind) = 0;
% 
% ind = find(Sp2);
% 
% Zplus = -2*eps*radius*(1./(Sp(ind)-radius) + 1./(Sp(ind)-radius*eps2));
% 
% F(ind) = 1./(exp(Zplus) + 1);

end

