function [a,b,at] = RD_spots2(cells_x,cells_y,cells_dists,...
    s,beta_init,beta_rand,Da,Db,max_iter,max_dist,change_thresh)
%
% [a,b] = RD_spots2(x,y,d,0.05,12,0.1,0.25,0.0625,20000,2);

[n_side,n_side] = size(cells_x);
a = 4*ones(n_side,n_side);
b = 4*ones(n_side,n_side);
beta = beta_init*ones(n_side,n_side);
beta = beta + beta_rand*(0.5-rand(n_side,n_side));
neighbors = cells_dists<max_dist&cells_dists>0;
trace_index = 1;
at(:,:,trace_index) = a;

done = 0;
iter = 0;

w1 = waitbar(0,'Spots');
while done==0
    waitbar(iter/max_iter);
    iter = iter + 1;
    a_change = 0;
    b_change = 0;
    for x = 1:n_side
        for y = 1:n_side
            index = RD_coords2index(x,y,n_side);
            nei_indexes = find(neighbors(index,:));
            if ~isempty(nei_indexes)
                a_c = a(index);
                b_c = b(index);
                num_nei = length(nei_indexes);
                nei_sum_a = 0;
                nei_sum_b = 0;
                for n = 1:num_nei
                    [x_ind,y_ind] = RD_index2coords(nei_indexes(n),n_side);
                    nei_sum_a = nei_sum_a + a(x_ind,y_ind);
                    nei_sum_b = nei_sum_b + b(x_ind,y_ind);
                end
                Lap_a = nei_sum_a - num_nei*a_c;
                Lap_b = nei_sum_b - num_nei*b_c;
                da = s*(16-a_c*b_c) + Da*Lap_a;
                db = s*(a_c*b_c-b_c-beta(x,y)) + Db*Lap_b;
                a_next = a_c + da;
                b_next = b_c + db;
                if a_next>0
                    a_change = a_change + abs(a_c-a_next);
                    a(x,y) = a_next;
                end
                if b_next>0
                    b_change = b_change + abs(b_c-b_next);
                    b(x,y) = b_next;
                end
                if (a_change<change_thresh)&(b_change<change_thresh)
                    done = 1;
                end
            end
        end
    end
    trace_index = trace_index + 1;
    at(:,:,trace_index) = a;
    if iter>max_iter
        done = 1;
    end
end
close(w1);
