function phi = probability(phi_w, t_low, t_high, w_s, w_t)
si = size(phi_w);
    function p = prob_pass(phi_w, threshold)
        p = phi_w;
        pw = 0;
        for x=1:si(1)
            for y=1:si(2)
                for z=1:si(3)
                    for t=1:si(4)
                        grad_sum = 0;
                        val = phi_w(x, y, z, t);
                        tot = 0;
                        if 1<x
                           grad_sum = grad_sum +  w_s*(phi_w(x-1, y, z, t) - val);
                           tot = tot + w_s;
                        end
                        if x<si(1)-1
                           grad_sum = grad_sum + w_s*(phi_w(x+1, y, z, t) - val);
                           tot = tot + w_s;
                        end
                        if 1<y
                           grad_sum = grad_sum + w_s*(phi_w(x, y-1, z, t) - val);
                           tot = tot + w_s;
                        end
                        if y<si(2)-1
                           grad_sum = grad_sum + w_s*(phi_w(x, y+1, z, t) - val);
                           tot = tot + w_s;
                        end
                        if 1<z
                           grad_sum = grad_sum + w_s*(phi_w(x, y, z-1, t) - val);
                           tot = tot +w_s;
                        end
                        if z<si(3)-1
                           grad_sum = grad_sum + w_s*(phi_w(x, y, z+1, t) - val);
                           tot =tot +w_s;
                        end
                        if 1<t
                           grad_sum = grad_sum + w_t*(phi_w(x, y, z, t-1) - val);
                           tot = tot +w_t;
                        end
                        if t<si(4)-1
                           grad_sum = grad_sum + w_t*(phi_w(x, y, z, t+1) - val);
                           tot = tot +w_t;
                        end
                        prob = grad_sum /(2*pi*tot);
                        if prob > threshold
                           p(x, y, z, t) = val + 2*pi;
                           pw = pw +1;
                        end
                        if prob < -threshold
                           p(x, y, z, t) = val - 2*pi;
                           pw = pw +1;
                        end
                    end
                end
            end
        end
        pw;
    end
phi = prob_pass(phi_w, t_low);
for i =1:9
    phi = prob_pass(phi, t_low);
end
for i=1:10
    phi = prob_pass(phi, t_high);
end
end