%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-10-26 14:56:26
 % @LastEditTime: 2021-10-26 16:09:49
%--------------------------
function H = Heaviside_simply(phi, radius)
    num_all=[1:length(phi(:))]'; 
    num1=find(phi>radius); 
    H(num1)=1; 
    num2=find(phi<-radius); 
    H(num2)=0; 
    num3=setdiff(num_all,[num1;num2]); 
    H(num3) = (radius - phi(num3)).^2 / 4 .* (2*radius + phi(num3)) / radius^3;
    % H(num3)=3*(1-alpha)/4*(phi(num3)/epsilon-phi(num3).^3/(3*(epsilon)^3))+ (1+alpha)/2; 
end