function Z = initial_conditions_4cone(Z,NX,NY,i0,j0,r,h)

for i = 1:NX
    for j = 1:NY
        di = mod(i-i0+NX/2,NX)-NX/2; % difference i-i0 but apply p.b.c. 
        dj = mod(j-j0+NY/2,NY)-NY/2;
        dr = sqrt((di^2)/6+(dj^2)/2);
        if (dr<r)
            Z(i,j) = Z(i,j) + h*(cos(dr/r*pi/2))^2;
        end
    end
end