function [pathout NodeP Nodeu Nodev] =...
    Solver_pathplan_NB(app,SetupScript,PlotScript) 
%function [pathout NodeP Nodeu Nodev] =...
%    Solver_pathplan(app,SetupScript,PlotScript) 
%
% This is the Narrowband LS solver for analytically specified flowfields.
% Pay special attention to app.NBlen and app.reinitLS
%
%INPUTS:
%       LS_MASK:    Interior point mask that defines forbidden regions,
%                   where the levelset is not allowed to go (but the 
%                   velocity field is)
%           app:	The application data structured. At a minimum it should 
%                   contain the following fields (which may be added through
%                   the SetupScript).
%          app.NBlen This is the length or width of the narrow band. It is
%                    specified as a fraction of (dx^2+dy^2)^0.5. So, for
%                    NBlen = 1, the narrow band will be about 1 grid cell
%                    on either side of the level set. For NB = 1000 (i.e.
%                    large), the entire field is calculated at each
%                    timestep. The default value is large. 
%          app.reinitLS:   When using the Narrow band (i.e. small
%                          app.NBlen), the level set within the narrow band
%                          needs to be reinitialized every now and then for
%                          an accurate calculation. This part of the
%                          process is, unfortunately, expensive...
%          app.Nx:         Number of interior nodes in x-direction (integer
%                           number)
%          app.Ny:         Number of interior nodes in y-direction (integer
%                           number)
%          app.dt:         Size of timestep (TotalTime/Nt)
%          app.U:          Function handle for u-velocity f(x,y,t)
%          app.V:          Function handle for v-velocity f(x,y,t)
%          app.pathin:     Structure containing information about the
%                           path-planning section
%          app.pathin.p:   Starting positions of assets p = [x;y]
%          app.pathin.mag: Magnitude of velocity that vehicle has
%          app.pathin.v:   Prescribed vector velocities of assets relative to the
%                     asset coordinate system. That is, the coordinate
%                     system is fixed in space, not taking into account the
%                     flowfield velocity. v = [vx;vy;t];
%       PlotIntrvl: Number of time integration intervals between of calling
%                   "plotscript."
%       PlotScript: String input of filename used to plot specific problem.
%                   Can also be used to save outputs.
%                       Note: Variables u, v, and P, will be available for
%                       plotting.
%       SetupScript: String input of filename used to setup the problem.
%           This file needs to correctly create the following variables:
%           dx:     Size of control volume in x-direction
%           dy:     Size of control volume in y-direction
%           Nt:         Number of timesteps (integer number)
%           bcsDrho:Density Dirichlet boundary conditions
%           bcsNrho:Density Neumann boundary conditions
%           NodeP:  Pressure node-numbering matrix (see NodePad.m)
%           Nodeu:  U-velocity node-numbering matrix (see NodePad.m)
%           Nodev:  V-velocity node-numbering matrix (see NodePad.m)
%           idsP:   Interior Pressure node id's (see NodePad.m)
%           idsu:   Interior U-velocity node id's (see NodePad.m)
%           idsv:   Interior V-velocity node id's (see NodePad.m)
%           XU:     Matrix of x-coordinates corresponding to U-velocity CV
%                  centers. XU(idsu) used as 'x' input to Fu. (see meshgrid)
%           XV:     Matrix of x-coordinates corresponding to V-velocity CV
%                  centers. YU(idsu) used as 'y' input to Fv. (see meshgrid)
%           YU:     Matrix of y-coordinates corresponding to U-velocity CV
%                  centers. XV(idsv) used as 'x' input to Fu. (see meshgrid)
%           YV:     Matrix of y-coordinates corresponding to V-velocity CV
%                  centers. YV(idsv) used as 'y' input to Fv. (see meshgrid)
%
% OUTPUTS:
%           pathout:Information about the path
%           PATHOUT:        Structure containing level set information
%   `       PATHOUT.MAG    Maximum magnitude relative velocity of vehicle
%            PATHOUT.N      Number of assets (vehicles)
%            PATHOUT.C      Cell array points defining level set curves (curves
%                            created using contourc.
%            PATHOUT.P      Points tracing current asset path (used for
%                            comparison and plots)
%            PATHOUT.T      Time at which each curve in pathout.c was generated
%            PATHOUT.V      Velocity vectors of flowfield at points in
%                            pathout.c
%           NodeP:  Pressure node-numbering matrix (see NodePad.m)
%           Nodeu:  U-velocity node-numbering matrix (see NodePad.m)
%           Nodev:  V-velocity node-numbering matrix (see NodePad.m)

%
%   Uses Matlab functions: 
%   See also: mk_DiffusionOperator, FixDiffuvDbcs, Grad2D, Div2D
%
%   Written by: Matt Ueckermann, Konur Yigit, Pierre Lermusiaux, Pat Haley


%Run setup script
eval(SetupScript)
Node2P = NodeP(2:end-1,2:end-1); % Convenience
if ~exist('bcsDrho','var'),bcsDrho=zeros(size(bcsDP));end
if ~exist('bcsNrho','var'),bcsNrho=zeros(size(bcsNP));end
if ~exist('rho','var'),rho=zeros(size(P));end
if ~isfield(app,'U')
    display('***Need to specify app.U to use this solver')
    pathout =[];
    return; 
end
if ~isfield(app,'V')
    display('***Need to specify app.V to use this solver')
    pathout =[];
    return; 
end
if ~exist('LS_MASK','var'),LS_MASK = ones(size(XP));end
if ~isfield(app,'NBlen')
    NBlen = max(app.Nx, app.Ny) * 100;
else
    NBlen = app.NBlen;
end
if ~isfield(app,'pathin') 
    app.pathin.p = []; 
    app.pathin.mag = 0; 
end
pathout.mag = app.pathin.mag;

%Find number of boundary conditions
Nbcs = length(bcsDrho) + length(bcsNrho) + length(bcsNR);

% Get interior ids
pid = Nbcs+1:length(rho);

reinitLS=0;

%Check if re-initialization was set. If not, set it to default value of 5
if ~isfield(app,'reinitLS')
    app.reinitLS = app.T/app.dt * 100;
end
pathcount=0;            % Counter for number of level sets stored  

indic = 0;

%% Solving system
T=0;k=0;
eval(PlotScript);
display('Paused at initial condition. Hit any key to continue');drawnow;pause
for k=1:Nt
    T = app.dt * k;
    %Get the velocity at the midpoint of the timestep
    u = app.U(XXU, YYU, T-app.dt/2);
    v = app.V(XXV, YYV, T-app.dt/2);
    
    if T>=0.1
      %  keyboard
    end
    
    % ---------------------------
    % MOTION IN NORMAL DIRECTION
    % ---------------------------
    % Get RHS contribution due to level set normal advection
    
    [RHSq NBids] = Narrowband_LS_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt, ...
                                         rho, u, v, NodeP, Nodeu, Nodev, ...
                                         idsP, Nbcs, length(bcsDrho), ...
                                         Nbcs, Nbcs, NBlen);                              
    RHSq = - app.pathin.mag * RHSq;
    
    
    if RHSq(56)~=0
        indic=1
        break
    end
    %Advance the levelset normal to itself by half a timestep
    rho(pid) = rho(pid) + app.dt/2 * RHSq .* LS_MASK(idsP);
    
    % --------------------------------
    % ADVECTION DUE TO VELOCITY FIELD
    % --------------------------------
    % Get RHS contribution due to motion by velocity field for a timestep
    
    RHSq = Narrowband_Advect_LS_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt,...
                                        rho, u, v, NodeP, Nodeu, Nodev,...
                                        idsP, Nbcs, length(bcsDrho),...
                                        Nbcs, Nbcs, NBlen);
    
    rho(pid) = rho(pid) + app.dt*RHSq.*LS_MASK(idsP);
    
    % ---------------------------
    % MOTION IN NORMAL DIRECTION
    % ---------------------------
    % Get RHS contribution due to level set normal advection
    
    [RHSq NBids] = Narrowband_LS_TVD_mex(app.Nx, app.Ny, dx, dy, app.dt,...
                                         rho, u, v, NodeP, Nodeu, Nodev,...
                                         idsP, Nbcs, length(bcsDrho),...
                                         Nbcs, Nbcs, NBlen);
    RHSq = - app.pathin.mag * RHSq;
    
    % Advance the levelset normal to itself by half a timestep
    rho(pid) = rho(pid) + app.dt/2 * RHSq .* LS_MASK(idsP);
    % =====================================================================
   
    % ---------------------------------------------------------------%
    % Extract, store and plot the zero level set at every time step: %
    % ---------------------------------------------------------------%
    
    figure(1)
   
    pathcount = pathcount+1;     % Update pathcount to store level sets
   % keyboard                          
    pathout.C{pathcount}=contour(XP,YP,rho(NodeP(2:app.Ny+1,2:app.Nx+1))...
                                ,[0.0 0.0],'k','linewidth',2);
    hold on
    plot(1,1,'o','markersize',15,'markeredgecolor','k','markerfacecolor','k')
    plot(1.8,1.8,'p','markersize',15,'markeredgecolor','k','markerfacecolor','k')
    xlabel('x')
    ylabel('y')
    tfin(k+1)=toc;
    title(sprintf('Time = %f, Bandwidth = %g, RS = %g, Runtime = %f',...
                  T,app.NBlen,app.reinitLS,tfin(k+1)))
    
    axis equal
    xlim([0 2.5])
    ylim([0 2.5])
    grid on
    hold off
    % =====================================================================
    
    %Evolve the vehicle trajectory
    for i = 1:pathout.N
        try
            vrel = [interp1(app.pathin.v(:,3,i),app.pathin.v(:,1,i), T),...
                interp1(app.pathin.v(:,3,i),app.pathin.v(:,2,i), T)];
        end
        APath(i,:) = APath(i,:) + app.dt / 2 *vrel;
        try
            vrel = [interp2(XU,YU,u(Nodeu(2:end-1,2:end-1)),APath(i,1),APath(i,2)),...
                interp2(XV,YV,v(Nodev(2:end-1,2:end-1)),APath(i,1),APath(i,2))];
        end
        APath(i,:) = APath(i,:) + app.dt * vrel;
        try
            vrel = [interp1(app.pathin.v(:,3,i),app.pathin.v(:,1,i), T),...
                interp1(app.pathin.v(:,3,i),app.pathin.v(:,2,i), T)];
        end
        APath(i,:) = APath(i,:) + app.dt / 2 * vrel;
    end
    for i = 1:pathout.N
        pathout.p(pathcount,:,i) = [APath(i,:) T];
    end
    
  pathout.t(pathcount) = T;
    
 % ========================================================================
   
% -----------------------------------------------
% Number of pieces in the zero level set contour:
% -----------------------------------------------
%   num = Number of discrete points in first contour
    numS = 1;                                 % Indicator variable     
    numC = 0;                                 % Number of contour pieces
    
    while numS < size(pathout.C{pathcount},2) % While indicator does not 
                                              % reach the end of array
        num = pathout.C{pathcount}(2,numS);   % Update 'num' each time a 
                                              % contour piece is accounted 
        numS = num+numS+1;                    % Update indicator position  
        numC=numC+1;                          % Update number of contours  
    end
% -------------------------------------------------------------------------    
    numS=1;                             
    CC=[];  % Stores points of the zero level set (same as pathout.C except
            % it does not have the breaks in between two pieces of the
            % contour.
    for i = 1:numC
        num = pathout.C{pathcount}(2,numS);
        ids = numS+1:numS+num;
        CC = [CC,pathout.C{pathcount}(1:2,ids)];
        numS = 1+numS+num;
    end
    
    pathout.v{pathcount}=full(double([interp2(XU,YU,u(Nodeu(2:end-1,2:end-1)),...
                              CC(1,:),CC(2,:));
                              interp2(XV,YV,v(Nodev(2:end-1,2:end-1)),...
                              CC(1,:),CC(2,:))]));
    % Flowfield velocities at zero level set position at each time
    
    c=sqrt(sum(pathout.v{pathcount}.^2));
    pathout.c{pathcount}=c;   % Magnitude of flow velocity
    % =====================================================================
    % --------------------------------------------------
    % Stop evolving level sets when they reach the goal:
    % --------------------------------------------------
   
    tmptrglev = interp2(XP,YP,rho(NodeP(2:end-1,2:end-1)),...
                        app.targets(:,1),app.targets(:,2));
    % Variable indicating the level set function at each of the
    % goal locations at every timestep.
    
    if k == 1
        tmptrglevindic = zeros(size(tmptrglev));
    end
    for ii = 1:length(tmptrglev)
        if tmptrglev(ii) < 0 && tmptrglevindic(ii) == 0
            pathout.crosscount(ii) = pathcount; 
            
            % pathout.crosscount indicates, for each goal location, when
            % the zero level set has crossed it. 
            
            tmptrglevindic(ii) = 1;
            
            % Variable that indicates whether a goal point has been reached
            % (value = 1) or not (value = 0). 
        end
    end
    
    if tmptrglev < 0   % Forward computation stops when all the entries of
                       % tmptrglev are negative.
        k = Nt;
        save NB_NO_F
        break
    end        
    
% =========================================================================    
% ------------------%
% REINITIALIZATION: %
% ------------------%   
   if T >=0.9
     %  keyboard
   end
    % Need to re-init one grid-cell larger than narrow band
    % Get IDs of points in narrow band, then IDs of points grid lengths
    % away from the narrow band points.
    
%     NB_IDS = find(NBids);
%     NB_IDS((NB_IDS <= 2*app.Ny) | (NB_IDS >= app.Ny * app.Nx - 2*app.Ny)) = [];
%     NB2ids = false(size(NBids));
%     
%     
%     NB2ids([NB_IDS + 1; NB_IDS-1; NB_IDS+app.Ny; NB_IDS-app.Ny;
%           NB_IDS + 2; NB_IDS-2; NB_IDS+2*app.Ny; NB_IDS-2*app.Ny ]) = true;
%     
%     NB2ids(NB_IDS) = false;
%     NB2_IDS = find(NB2ids);
%     
%     % Need to prepare the outside cells at every timestep
%     % if NBlen < 2
%     D = rho(NodeP(2:end-1,2:end-1));
%     numS = 1;
%     for i = 1:numC
%         num = pathout.C{pathcount}(2,numS);
%         D(NB2ids) = min(dpoly([XP(NB2ids) YP(NB2ids)],...
%             pathout.C{pathcount}(:,numS+1:numS+num)'),D(NB2ids));
%         numS = numS+num+1;
%     end
%     rho(pid) = D(idsP);
%     %reinitLS=0;
%     % end
%     
%     NBids([NB_IDS + 1; NB_IDS-1; NB_IDS+app.Ny; NB_IDS-app.Ny;
%           NB_IDS + 2; NB_IDS-2; NB_IDS+2*app.Ny; NB_IDS-2*app.Ny ]) = true;
%     reinitLS=reinitLS+1;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % DRASTIC CHANGES TO NB _ TAPOVAN!!
    
    % First Step in the two step update process for adding points inside
    % the NB. 
    
    % Manually remove a layer of two grid cells from the left, right, top 
    % and bottom edges of domain from NB. Also remove points which lie
    % inside obstacles. 
   
    NB_IDS = find(NBids);
    NB_IDS((NB_IDS <= 2*app.Ny) | (NB_IDS >= app.Ny * app.Nx - 2*app.Ny)) = [];
    NB_IDS((mod(NB_IDS,app.Ny)<=2) | (mod(NB_IDS,app.Ny)>=app.Ny-1)) = [];
    NB_IDS(Node2P(NB_IDS)<=Nbcs) = []; % Inside Obstacles.
    NBids = false(size(NBids));
    NBids(NB_IDS) = true; 
   % Preliminary Check: Checks if points around a NB point are in the NB or
   % not. If not, the current point is removed from the band.
    NB_IDS = find(NBids);
    % IF POINTS ARE ISOLATED, REMOVE THEM:
%     ind = 1;
%     l = length(NB_IDS);
%     while ind <= l
%         dumv=NB_IDS(ind);
%         if ((NBids(dumv+1)*NBids(dumv+2)==0) &&...
%                 (NBids(dumv-1)*NBids(dumv-2)==0) && ...
%                 (NBids(dumv+app.Ny)*NBids(dumv+2*app.Ny)==0) &&...
%                 (NBids(dumv-app.Ny)*NBids(dumv-2*app.Ny)==0))
%             disp('Cleaning ...')
%             NBids(dumv) = 0;
%             NB_IDS(NB_IDS == dumv) = [];
%             int = 0;
%             l = length(NB_IDS);
%         end
%         ind = ind+1;
%     end
    
    NB_IDS = find(NBids);   
    NB1ids = false(size(NBids)); % Initialize band of NB1 ids. 
    NB1ids([NB_IDS + 1; NB_IDS-1; NB_IDS+app.Ny; NB_IDS-app.Ny]) = true;
    % Setting NB1ids one grid cell around NBids to 1. 
    
    % Remove the points of NB1ids which fall inside obstacles or within two
    % grid cells away from the open boundaries.
    NB1_IDS = find(NB1ids);
    NB1_IDS((NB1_IDS <= 2*app.Ny) | (NB1_IDS >= app.Ny * app.Nx - 2*app.Ny)) = [];
    NB1_IDS((mod(NB1_IDS,app.Ny)<=2) | (mod(NB1_IDS,app.Ny)>=app.Ny-1)) = [];
    NB1_IDS(Node2P(NB1_IDS)<=Nbcs) = []; % Inside Obstacles.
    NB1ids = false(size(NBids)); % Initialize band of NB1 ids. 
    NB1ids(NB1_IDS) = true;  % Update NB1ids based on NB1_IDS.
    
    NB1ids(NB_IDS) = false;  % Set NB1ids of all NBids to zero to include 
                             % only feasible and incremental points. 
    
    NB1_IDS = find(NB1ids);  % Find updated IDS on NB1 points.
    
    % Now, remove NB1 points which are isolated. 
    ind = 1; 
    l = length(NB1_IDS);
    while ind <=l
        dumv=NB1_IDS(ind);
        if ((NBids(dumv+1)*NBids(dumv+2)==0) &&...
                (NBids(dumv-1)*NBids(dumv-2)==0) && ...
                (NBids(dumv+app.Ny)*NBids(dumv+2*app.Ny)==0) &&...
                (NBids(dumv-app.Ny)*NBids(dumv-2*app.Ny)==0))
           % disp('Cleaning ...')
            NB1ids(dumv) = 0;    % Updated
            NB1_IDS(NB1_IDS==dumv) = [];
            ind = 0;
            l = length(NB1_IDS);
        end
        ind = ind+1;
    end 
   NB1_IDS = find(NB1ids);       % Final NB1_IDS 
   % Now, NB1_IDS and NB1ids are prepared.
    
    for ind = 1:length(NB1_IDS)
        dumv = NB1_IDS(ind);
        if (dumv+2*app.Ny<=Node2P(end))   % For safety near boundaries
            rhor1 = rho(Node2P(dumv+app.Ny));
            rhor2 = rho(Node2P(dumv+2*app.Ny));
        else
            rhor1 = 0;rhor2 = 0;
        end
        rhou1 = rho(Node2P(dumv-1));
        rhou2 = rho(Node2P(dumv-2));
        if (dumv-2*app.Ny > 0)            % For safety near boundaries   
            rhol1 = rho(Node2P(dumv-app.Ny));
            rhol2 = rho(Node2P(dumv-2*app.Ny));
        else
            rhol1 = 0;rhol2 = 0;
        end
        rhod1 = rho(Node2P(dumv+1));
        rhod2 = rho(Node2P(dumv+2));
        dum =  ((2*rhor1 - rhor2)*NBids(min(dumv+app.Ny,length(NBids)))*NBids(min(dumv+2*app.Ny,length(NBids))) + ...
                (2*rhod1 - rhod2)*NBids(dumv+1)     *NBids(dumv+2) + ...
                (2*rhou1 - rhou2)*NBids(dumv-1)     *NBids(dumv-2) + ...
                (2*rhol1 - rhol2)*NBids(max(dumv-app.Ny,1))*NBids(max(dumv-2*app.Ny,1)));
         
        if dum > 0
            rho(Node2P(dumv)) = dum/(NBids(min(dumv+app.Ny,length(NBids)))*NBids(min(dumv+2*app.Ny,length(NBids))) + ...
                NBids(dumv+1)     *NBids(dumv+2) + ...
                NBids(dumv-1)     *NBids(dumv-2) + ...
                NBids(max(dumv-app.Ny,1))*NBids(max(dumv-2*app.Ny,1)));
        end
        if rho(Node2P(dumv))<0
        disp('You are doomed ... ');
        end
   
    end 
    NBids(NB1_IDS) = true;  % Finally, add NB1 points to current NB.
    NB_IDS = find(NBids);
%     % IF POINTS ARE ISOLATED, REMOVE THEM:
%     for ind = 1:length(NB_IDS)
%         dumv=NB_IDS(ind);
%         if ((NBids(dumv+1)*NBids(dumv+2)==0) &&...
%                 (NBids(dumv-1)*NBids(dumv-2)==0) && ...
%                 (NBids(dumv+app.Ny)*NBids(dumv+2*app.Ny)==0) &&...
%                 (NBids(dumv-app.Ny)*NBids(dumv-2*app.Ny)==0))
%             disp('Cleaning ...')
%             NBids(dumv) = 0;
%         end
%     end    
    % End of First Step
    
    
    % Second Step in the two step update process for adding points inside
    % the NB. 
    
    % Manually remove a layer of two grid cells from the left, right, top 
    % and bottom edges of domain from NB. Also remove points which lie
    % inside obstacles. 
    NB_IDS = find(NBids);
    NB_IDS((NB_IDS <= 2*app.Ny) | (NB_IDS >= app.Ny * app.Nx - 2*app.Ny)) = [];
    NB_IDS((mod(NB_IDS,app.Ny)<=2) | (mod(NB_IDS,app.Ny)>=app.Ny-1)) = [];
    NB_IDS(Node2P(NB_IDS)<=Nbcs) = []; % Inside Obstacles.
    NBids = false(size(NBids));
    NBids(NB_IDS) = true; 
   % Preliminary Check: Checks if points around a NB point are in the NB or
   % not. If not, the current point is removed from the band.
    NB_IDS = find(NBids);
%    IF POINTS ARE ISOLATED, REMOVE THEM:
%     ind = 1;
%     l = length(NB1_IDS);
%     while ind <= l
%         dumv=NB_IDS(ind);
%         if ((NBids(dumv+1)*NBids(dumv+2)==0) &&...
%                 (NBids(dumv-1)*NBids(dumv-2)==0) && ...
%                 (NBids(dumv+app.Ny)*NBids(dumv+2*app.Ny)==0) &&...
%                 (NBids(dumv-app.Ny)*NBids(dumv-2*app.Ny)==0))
%            % disp('Cleaning ...')
%             NBids(dumv) = 0;
%             NB_IDS(NB_IDS == dumv) = [];
%             int = 0;
%             l = length(NB_IDS);
%         end
%         ind = ind + 1;
%     end
% %     
    NB_IDS = find(NBids);   
    NB1ids = false(size(NBids)); % Initialize band of NB1 ids. 
    NB1ids([NB_IDS + 1; NB_IDS-1; NB_IDS+app.Ny; NB_IDS-app.Ny]) = true;
    % Setting NB1ids one grid cell around NBids to 1. 
    
    % Remove the points of NB1ids which fall inside obstacles or within two
    % grid cells away from the open boundaries.
    NB1_IDS = find(NB1ids);
    NB1_IDS((NB1_IDS <= 2*app.Ny) | (NB1_IDS >= app.Ny * app.Nx - 2*app.Ny)) = [];
    NB1_IDS((mod(NB1_IDS,app.Ny)<=2) | (mod(NB1_IDS,app.Ny)>=app.Ny-1)) = [];
    NB1_IDS(Node2P(NB1_IDS)<=Nbcs) = []; % Inside Obstacles.
    NB1ids = false(size(NBids)); % Initialize band of NB1 ids.
    NB1ids(NB1_IDS) = true;  % Update NB1ids based on NB1_IDS.
    
    NB1ids(NB_IDS) = false;  % Set NB1ids of all NBids to zero to include 
                             % only feasible and incremental points. 
    
    NB1_IDS = find(NB1ids);  % Find updated IDS on NB1 points.
    
    % Now, remove NB1 points which are isolated. 
    ind = 1; 
    l = length(NB1_IDS);
    while ind <= l
        dumv=NB1_IDS(ind);
        if ((NBids(dumv+1)*NBids(dumv+2)==0) &&...
                (NBids(dumv-1)*NBids(dumv-2)==0) && ...
                (NBids(dumv+app.Ny)*NBids(dumv+2*app.Ny)==0) &&...
                (NBids(dumv-app.Ny)*NBids(dumv-2*app.Ny)==0))
            %disp('Cleaning ...')
            NB1ids(dumv) = 0;    % Updated 
            NB1_IDS(NB1_IDS == dumv) = [];
            ind = 0;
            l = length(NB1_IDS);
        end
        ind = ind+1;
    end 
   NB1_IDS = find(NB1ids);       % Final NB1_IDS 
   % Now, NB1_IDS and NB1ids are prepared.
    
    for ind = 1:length(NB1_IDS)
        dumv = NB1_IDS(ind);
        if (dumv+2*app.Ny<=Node2P(end))   % For safety near boundaries
            rhor1 = rho(Node2P(dumv+app.Ny));
            rhor2 = rho(Node2P(dumv+2*app.Ny));
        else
            rhor1 = 0;rhor2 = 0;
        end
        rhou1 = rho(Node2P(dumv-1));
        rhou2 = rho(Node2P(dumv-2));
        if (dumv-2*app.Ny > 0)            % For safety near boundaries   
            rhol1 = rho(Node2P(dumv-app.Ny));
            rhol2 = rho(Node2P(dumv-2*app.Ny));
        else
            rhol1 = 0;rhol2 = 0;
        end
        rhod1 = rho(Node2P(dumv+1));
        rhod2 = rho(Node2P(dumv+2));
       
       dum = ((2*rhor1 - rhor2)*NBids(min(dumv+app.Ny,length(NBids)))*NBids(min(dumv+2*app.Ny,length(NBids))) + ...
                            (2*rhod1 - rhod2)*NBids(dumv+1)     *NBids(dumv+2) + ...
                            (2*rhou1 - rhou2)*NBids(dumv-1)     *NBids(dumv-2) + ...
                            (2*rhol1 - rhol2)*NBids(max(dumv-app.Ny,1))*NBids(max(dumv-2*app.Ny,1)));
        if dum > 0
            rho(Node2P(dumv)) = dum/(NBids(min(dumv+app.Ny,length(NBids)))*NBids(min(dumv+2*app.Ny,length(NBids))) + ...
                NBids(dumv+1)     *NBids(dumv+2) + ...
                NBids(dumv-1)     *NBids(dumv-2) + ...
                NBids(max(dumv-app.Ny,1))*NBids(max(dumv-2*app.Ny,1)));
        else
          %  NB1_IDS(NB1_IDS == dumv) = [];
          %  NB1ids(dumv) = 0;
        end
        if rho(Node2P(dumv))<0
        disp('You are doomed ... ');
        end
    end 
        
    NBids(NB1_IDS) = true; 
   % Finally, add NB1 points to current NB.%    % Preliminary Check: Checks if points around a NB point are in the NB or
   % not. If not, the current point is removed from the band.
    NB_IDS = find(NBids);
%     % IF POINTS ARE ISOLATED, REMOVE THEM:
%     for ind = 1:length(NB_IDS)
%         dumv=NB_IDS(ind);
%         if ((NBids(dumv+1)*NBids(dumv+2)==0) &&...
%                 (NBids(dumv-1)*NBids(dumv-2)==0) && ...
%                 (NBids(dumv+app.Ny)*NBids(dumv+2*app.Ny)==0) &&...
%                 (NBids(dumv-app.Ny)*NBids(dumv-2*app.Ny)==0))
%             disp('Cleaning ...')
%             NBids(dumv) = 0;
%         end
%     end
    % End of Second Step
      
%     NB1_IDS(Node2P(NB1_IDS+1)<=Nbcs) = []; % Inside Obstacles.
%     NB1_IDS(Node2P(NB1_IDS-1)<=Nbcs) = []; % Inside Obstacles.
%     NB1_IDS(Node2P(NB1_IDS+app.Ny)<=Nbcs) = []; % Inside Obstacles.
%     NB1_IDS(Node2P(NB1_IDS-app.Ny)<=Nbcs) = []; % Inside Obstacles.
   
    
   
       
    
%     NB_IDS(Node2P(NB_IDS)<=Nbcs) = []; % Inside Obstacles.
%     NB_IDS(Node2P(NB_IDS+1)<=Nbcs) = []; % Inside Obstacles.
%     NB_IDS(Node2P(NB_IDS-1)<=Nbcs) = []; % Inside Obstacles.
%     NB_IDS(Node2P(NB_IDS+app.Ny)<=Nbcs) = []; % Inside Obstacles.
%     NB_IDS(Node2P(NB_IDS-app.Ny)<=Nbcs) = []; % Inside Obstacles.
    
            
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     
    
    if k>0 && (reinitLS>=app.reinitLS || k * app.dt * app.pathin.mag < norm([dx, dy].^2))
        D = rho(NodeP(2:end-1,2:end-1));
        numS = 1;
        for i = 1:numC
            num = pathout.C{pathcount}(2,numS);
            D(NBids) = min(dpoly([XP(NBids) YP(NBids)],...
                pathout.C{pathcount}(:,numS+1:numS+num)'),D(NBids));
            numS = numS+num+1;
        end
        rho(pid) = D(idsP);
        reinitLS=0;
    end
    
    %Plot solution
    if ~mod(k,app.PlotIntrvl(1))
        eval(PlotScript);drawnow
    end
   
end
if (~mod(k,app.PlotIntrvl(1))) == 0
    eval(PlotScript);drawnow
end
