%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-10-26 14:56:26
 % @LastEditTime: 2022-02-01 11:56:36
%--------------------------
function H = Heaviside_simply(phi, radius)

    H = (radius - phi).^2 / 4 .* (2*radius + phi) / radius^3;
    H(phi>radius)=1; 
    H(phi<-radius)=1e-3;       
    % H(num3)=3*(1-alpha)/4*(phi(num3)/epsilon-phi(num3).^3/(3*(epsilon)^3))+ (1+alpha)/2; 
end