function [trace,pts,pts_v] = BR_Cavity_no_movie_cont(x_min,x_max,y_min,...
    y_max,del_t,max_steps,pts,pts_v)
% BR_Cavity - simulate vortex formation in 2D
% On input:
%     x_min (float): min x value for 2D area
%     x_max (float): max x value for 2D area
%     y_min (float): min y value for 2D area
%     y_max (float): max y value for 2D area
%     del_t (float): time step
%     max_steps (int): number of steps for simulation
% On output:
%     trace (max_steps+1x10 array): each row has trace of the trace
%        particle (selected by variable t_pt below)
%       (:,1): F_long_range x
%       (:,2): F_long_range y
%       (:,3): F_short_range x
%       (:,4): F_short_range y
%       (:,5): F x
%       (:,6): F y
%       (:,7): r x
%       (:,8): r y
%       (:,9): v x
%       (:,10): v y
%     pts (2576x2 array): particle locations
%     vels (2576x2 array): particle velocities
% Method:
%     See Chapter 7, Particle Models
% Call:
%     [t,pts,vels] = BR_Cavity(-6.25,6.25,-6.25,6.25,0.0001,12000);
% Author:
%     T. Henderson
%     UU
%     Summer 2016; revised Spring 2017
%

box = [-6.25,6.25; 6.25,6.25; 6.25,-6.25; -6.25,-6.25];
box = [x_min,y_max; x_max,y_max; x_max,y_min; x_min,y_min];
trace = zeros(max_steps+1,10);
t_pt = 1;

MIN_Y = y_min; %-6.25;
MAX_Y = y_max % 6.25;
MIN_X = x_min; %-6.25;
MAX_X = x_max; %6.25;

MAX_ABS_V = 2000;
INT_MIN = 0.001;
D = 0.35;
V = -10;
V_VEL = [-10,0];
J = 100;
PT_MIN_DIST = 0.001;
LINE_L = [MIN_X,MIN_Y,0;MIN_X,MAX_Y,0];
LINE_U = [MIN_X,MAX_Y,0;MAX_X,MAX_Y,0];
LINE_R = [MAX_X,MIN_Y,0;MAX_X,MAX_Y,0];
LINE_D = [MIN_X,MIN_Y,0;MIN_X,MIN_Y,0];

num_pts = 2576;
pts_last = zeros(num_pts,2);

% restricted set
num_pts2 = 0;
pts2 = [];
for pp = 1:num_pts
    if pts(pp,1)>=MIN_X&pts(pp,1)<=MAX_X&pts(pp,2)>=MIN_Y&pts(pp,2)<=MAX_Y
        num_pts2 = num_pts2 + 1;
        pts2(num_pts2,:) = pts(pp,:);
    end
end

pts_last = pts;
pts = pts2;
pts_last2 = pts;
num_pts = num_pts2;

%pts_v = 0.001*(0.5-rand(num_pts,2));
%pts_v = zeros(num_pts,2);

%num_frames = 0;
%clf
%plot(MAX_X+0.1,MAX_Y+0.1,'w.');
%hold on
%plot(MIN_X-0.1,MIN_Y-0.1,'w.');
%quiver(pts(:,1),pts(:,2),pts_v(:,1),pts_v(:,2));
%%plot(pts(:,1),pts(:,2),'k.');
%num_frames = num_frames + 1;
%text(-1,-7,'Frame 1');
%M(num_frames) = getframe;

m = 1;
g = 980;
G = 0;
H = 100;
pp = 3;
qq = 5;

max_time = del_t*max_steps;
v = 0;
step = 1;
trace(1,:) = ...
    [0,0,0,0,0,0,pts(t_pt,1),pts(t_pt,2),pts_v(t_pt,1),pts_v(t_pt,2)];
wb = waitbar(0,'Cavity Flow No Movie Cont');
for st = 1:max_steps
    t = del_t*st;
    waitbar(t/max_time);
    step = step + 1
    ns = createns(pts,'nsmethod','kdtree');
    idx = rangesearch(ns,pts,0.351);
    new_pts = pts;
    for p = 1:num_pts
        pt = pts(p,:);
        nei = idx{p};
        nei = setdiff(nei,[p]);
        if ~isempty(nei)
            num_nei = length(nei);
            nei_pts = zeros(num_nei,2);
            dists = zeros(num_nei,1);
            for n = 1:num_nei
                nei_pts(n,1) = pts(nei(n),1);
                nei_pts(n,2) = pts(nei(n),2);
                dists(n) = norm(nei_pts(n,:)-pt);
            end
            F_long_range = [0;-g];
            F_short_range = [0;0];
            for n = 1:num_nei
                dn = dists(n);
                if dn<D&dn>0
                    F_short_range = F_short_range + ((G/(dn^pp))...
                        + (H/(dn^qq)))*(pt'-nei_pts(n,:)')/dn;
                end
            end
            F = F_long_range + F_short_range;
            a = F/m;
            pts_v(p,:) = pts_v(p,:) + del_t*a';
            if norm(pts_v(p,:))>MAX_ABS_V
                pts_v(p,:) = MAX_ABS_V*pts_v(p,:)/norm(pts_v(p,:));
            end
            r = pt' + del_t*pts_v(p,:)';
            [new_pt,vel] = BR_bounded_motion(box,pt,pts_v(p,:),del_t);
            r = new_pt;
            pts_v(p,:) = vel;
            if p==t_pt
                trace(step,1) = F_long_range(1);
                trace(step,2) = F_long_range(2);
                trace(step,3) = F_short_range(1);
                trace(step,4) = F_short_range(2);
                trace(step,5) = F(1);
                trace(step,6) = F(2);
                trace(step,7) = r(1);
                trace(step,8) = r(2);
                trace(step,9) = pts_v(p,1);
                trace(step,10) = pts_v(p,2);
            end
            new_pts(p,:) = r';
        end
    end
    pts = new_pts;
%    if mod(step,J)==0
%        pts_v2(:,1) = (pts(:,1)-pts_last2(:,1))/(J*del_t);
%        pts_v2(:,2) = (pts(:,2)-pts_last2(:,2))/(J*del_t);
%        clf
%        plot(MAX_X+0.1,MAX_Y+0.1,'w.');
%        hold on
%        plot(MIN_X-0.1,MIN_Y-0.1,'w.');
%        plot(pts(:,1),pts(:,2),'k.');
%        quiver(pts_last2(:,1),pts_last2(:,2),pts_v2(:,1),pts_v2(:,2));
%        num_frames = num_frames + 1;
%        text((MAX_X-MIN_X)/2,MIN_Y-0.5,['Frame ',num2str(num_frames)]);
%        M(num_frames) = getframe;
%        pts_last2 = pts;
%    end
end
close(wb);
