function repForce=get_repForce(xPos, yPos, p)
%{
Determines repulsive effect for direct cell-cell interaction

*Input*
xPos, yPos .. current cell positions
p .. parameters structure

*Output*
repForce .. Nx2 array of repulsive force

Written 12/2021 by Angelika Manhart
%}
    repRad=p.repRad;
    Lx=p.Lx; Ly=p.Ly;
 
    if p.optBC==0 % periodic BCs
        
        % copy particles to account for periodic BCs
        delta=0.01/p.Lx;

        % who is near a boundary
        isNearBoundary_L = xPos<=repRad+delta;
        isNearBoundary_R = xPos>=Lx-repRad-delta;
        isNearBoundary_D = yPos<=repRad+delta;
        isNearBoundary_U = yPos>=Ly-repRad-delta;

        % who is near an edge
        isNearEdge_LU=isNearBoundary_L & isNearBoundary_U;
        isNearEdge_LD=isNearBoundary_L & isNearBoundary_D;
        isNearEdge_RU=isNearBoundary_R & isNearBoundary_U;
        isNearEdge_RD=isNearBoundary_R & isNearBoundary_D;

        % extend
        xPosExt=[xPos; xPos(isNearBoundary_L)+Lx; xPos(isNearBoundary_R)-Lx; xPos(isNearBoundary_D); xPos(isNearBoundary_U);...
                xPos(isNearEdge_LU)+Lx; xPos(isNearEdge_LD)+Lx; xPos(isNearEdge_RU)-Lx; xPos(isNearEdge_RD)-Lx];
        yPosExt=[yPos; yPos(isNearBoundary_L); yPos(isNearBoundary_R); yPos(isNearBoundary_D)+Ly; yPos(isNearBoundary_U)-Ly;...
                yPos(isNearEdge_LU)-Ly; yPos(isNearEdge_LD)+Ly; yPos(isNearEdge_RU)-Ly; yPos(isNearEdge_RD)+Ly];
            
    elseif p.optBC==1 % non periodic BCs
        
        xPosExt=xPos;
        yPosExt=yPos;
        
    end
            
    
    % how does the repulsion force depend on distance?              
    distPush=@(d) p.repStrength*(repRad-d)/repRad; % p.repStrength is the maximal pushing strength, which decays to 0
    
    % find neighbors
    neighborID_rep  = rangesearch([xPosExt yPosExt],[xPos yPos],repRad);

    % initialize
    repForce=zeros(p.N,2);
    
    % loop through particles
    for i=1:p.N
        
        % who are my neighbors
        iNeigh_rep=neighborID_rep{i};
        
        % extract vectors
        iVecX=xPos(i)-xPosExt(iNeigh_rep);
        iVecY=yPos(i)-yPosExt(iNeigh_rep);
        
        % normalize them
        iNorm=sqrt(iVecX.^2+iVecY.^2);
        
        % normalized vector
        iVecXnorm=iVecX./iNorm;
        iVecYnorm=iVecY./iNorm;
        
        % get distance dependence force
        iDistPush=distPush(iNorm);
        
        % calculate pushing force on swimmers
        iPushX=nansum(iVecXnorm.*iDistPush);
        iPushY=nansum(iVecYnorm.*iDistPush);
        
        % store forces
        repForce(i,:)=[iPushX iPushY];
        
    end
    
  
        
        
   