function sigConc_p1=get_sigConc_p1(sigConc, actConcUse, xPos, yPos, p)
%{
Solves the PDE for the external signal dynamics
Called by main_Dicty.m

*Input*
sigConc .. array of current signal concentration
actConcUse .. Nx1 vector indicating which cells are active and secrete
signal
xPos, yPos .. current cell positions
p .. parameters structure


*Output*
sigConc_p1 .. array of signal concentration at the next timestep (p1 stands for "+1")

Written 12/2021 by Angelika Manhart
%}


if p.sigSpreadRad==0 % all the signal of an active cell is released onto one grid point

    % find positions of cells on grid
    [binX, binY]=get_cellPosOnGrid(xPos, yPos, p);
    
    % map the concentration of activator onto the grid (Note: not well-written, could be made faster and
    % included in the function above!)
    
    ind=sub2ind(size(sigConc),binY, binX);
    actConcOnGrid=zeros(size(sigConc));
    
    for i=1:p.N
        actConcOnGrid(ind(i))=actConcOnGrid(ind(i))+actConcUse(i);
    end

    actConcOnGrid=actConcOnGrid/(p.dx*p.dy);
    
else 
    %{
    The signal of an active cell is spread over a box of side length 2*p.sigSpreadRad around the cells
    
    Note: This is currently only writte for no-flux BCs, for periodic BCs,
    the code would have to be adapted here!    
    %}
    
    if p.optBC==0
       disp('You are using periodic BCs, the signal spreading code in get_sigConc_p1 should be adapted!')
       pause         
    end
    
    %{
     the following defines the spreading function. For numerical efficiency, we use
     the product of 2 1D bumps, as oppsed to 1 2D bump
    %}
    
    intVal1D=0.443994; % normalization constant
    
    % define bump function
    bump1D=@(r) (1/(intVal1D*p.sigSpreadRad))*exp(-1./(1-((r/p.sigSpreadRad).^2))).*(abs(r)<p.sigSpreadRad);

    actConcOnGrid=zeros(size(sigConc)); % initialize
    nrSpreadRad=floor(p.sigSpreadRad/p.dx); % spread radius in grid points
    
    % loop over active cells only
    for iCell=(find(actConcUse>0))' 
        
        s=[(xPos(iCell)+p.Lx/2)/p.dx (yPos(iCell)+p.Ly/2)/p.dy];
        
        ss1=ceil(s(1)); ss2=ceil(s(2)); % find lower bordering grid indices
        
        i1=(ss1-nrSpreadRad+1):(ss1+nrSpreadRad); i2=(ss2-nrSpreadRad+1):(ss2+nrSpreadRad); % define affected grid indices
        i1(i1<1)=[]; i1(i1>p.nX)=[]; % delete those outside of domain
        i2(i2<1)=[]; i2(i2>p.nY)=[]; % delete those outside of domain
        
        r1=xPos(iCell)-p.xVals(i1); r2=yPos(iCell)-p.yVals(i2); % calculate distances to cell
        
        bb=bump1D(abs(r1))'*bump1D(abs(r2)); % evaluate bump function
        
        actConcOnGrid(i2,i1)=actConcOnGrid(i2,i1)+bb'; % update added signal to grid
    end
    
end

% define added signal concentration (source term in signal PDE)
addedSigConc=p.sigSecret.*actConcOnGrid;

%% Solve PDE

switch p.optSigDisc
    
    case 0 % fully explicit

        % A_{i,j+1}=circishift(A,-1,2) , affects x-coord
        % A_{i,j-1}=circishift(A,1,2)  , affects x-coord
        % A_{i+1,j}=circishift(A,-1,1) , affects y-coord
        % A_{i-1,j1}=circishift(A,1,1) , affects y-coord
        
        % simple finite difference disc of diffusion term            
        diffTerm=(circshift(sigConc,1,2)+circshift(sigConc,-1,2)-2*sigConc)/p.dx/p.dx+(circshift(sigConc,1,1)+circshift(sigConc,-1,1)-2*sigConc)/p.dy/p.dy;    % diffusion in y

        %  stability condition needs to be fulfilled!
        if (1-2*p.dt*p.sigDiff*(1/p.dx^2+1/p.dy^2))<0
            disp('explicit scheme not stable!')
        end
        
        sigConc_p1=sigConc+p.dt*addedSigConc+...
                        p.dt*p.sigDiff*diffTerm-p.dt*p.sigDecay*sigConc;

        if p.optBC==1 % no-flux BCs
            
            % Left/Right
            sigConc_p1(2:p.nY-1,1)=sigConc_p1(2:p.nY-1,2); 
            sigConc_p1(2:p.nY-1,end)=sigConc_p1(2:p.nY-1,end-1); 
            
            % Upper/Lower
            sigConc_p1(1,2:p.nX-1)=sigConc_p1(2, 2:p.nX-1); 
            sigConc_p1(end,2:p.nX-1)=sigConc_p1(end-1, 2:p.nX-1); 
            
            % Corners - assume diagonally directed normal vectors
            sigConc_p1(1,1)=(p.dy*sigConc_p1(1,2)+p.dx*sigConc_p1(2,1))/(p.dx+p.dy);
            sigConc_p1(1,end)=(p.dy*sigConc_p1(1,end-1)+p.dx*sigConc_p1(2,end))/(p.dx+p.dy);
            sigConc_p1(end,end)=(p.dy*sigConc_p1(end,end-1)+p.dx*sigConc_p1(end-1,end))/(p.dx+p.dy);
            sigConc_p1(end,1)=(p.dy*sigConc_p1(end,2)+p.dx*sigConc_p1(end-1,1))/(p.dx+p.dy);
            
        end
                    
    case 1 % full implicit (slower per timestep, but allows larger timesteps)
        
        % prepare coefficients
        ax=p.dt*p.sigDiff/(p.dx^2); ay=p.dt*p.sigDiff/(p.dy^2);
        c=1+p.dt*p.sigDecay+2*ax+2*ay;
        
        % set up matrix
        NN=p.nX*p.nY;
        
        A = sparse(NN,NN);
                
        % set up RHS
        b=reshape((sigConc+p.dt*addedSigConc)',NN,1);
        
        if p.optBC==0 % PBC
            
            % set-up nX times nX matrices
            aa1=c*eye(p.nX);
            aa1=aa1+diag(-ax*ones(p.nX-1,1),1)+diag(-ax*ones(p.nX-1,1),-1);
            aa1(1,end)=-ax;
            aa1(end,1)=-ax;
            
            aa2=-ay*eye(p.nX);
            
            % populate matrix
            for ii=1:p.nY
                A(p.nX*(ii-1)+1:p.nX*(ii-1)+p.nX,p.nX*(ii-1)+1:p.nX*(ii-1)+p.nX)=aa1;
                if ii<=p.nY-1
                    A(p.nX*(ii-1)+1+p.nX:p.nX*(ii-1)+2*p.nX, p.nX*(ii-1)+1:p.nX*(ii-1)+p.nX)=aa2;
                    A(p.nX*(ii-1)+1:p.nX*(ii-1)+p.nX, p.nX*(ii-1)+1+p.nX:p.nX*(ii-1)+2*p.nX)=aa2;
                end
            end
            
            A(1:p.nX,NN-p.nX+1:NN)=aa2;
            A(NN-p.nX+1:NN,1:p.nX)=aa2;
            
        elseif p.optBC==1 % no-flux BC
        
            % set-up nX times nX matrices
            aa1=c*eye(p.nX);
            aa1=aa1+diag(-ax*ones(p.nX-1,1),1)+diag(-ax*ones(p.nX-1,1),-1);
            
            aa2=-ay*eye(p.nX);
            
            for ii=1:p.nY
                A(p.nX*(ii-1)+1:p.nX*(ii-1)+p.nX,p.nX*(ii-1)+1:p.nX*(ii-1)+p.nX)=aa1;
                if ii<=p.nY-1
                    A(p.nX*(ii-1)+1+p.nX:p.nX*(ii-1)+2*p.nX, p.nX*(ii-1)+1:p.nX*(ii-1)+p.nX)=aa2;
                    A(p.nX*(ii-1)+1:p.nX*(ii-1)+p.nX, p.nX*(ii-1)+1+p.nX:p.nX*(ii-1)+2*p.nX)=aa2;
                end
            end

            % Left/Right
            rowNr=p.nX+1:p.nX:(p.nY-1)*p.nX;
            colNr=p.nX+1:p.nX:(p.nY-1)*p.nX;
            val1=1; val2=-1;            
            ind=sub2ind(size(A),rowNr,colNr);
            A(rowNr,:)=0; b(rowNr)=0;             
            A(ind)=val1;
            A(ind+NN)=val2;
         
            rowNr=2*p.nX:p.nX:p.nY*p.nX-p.nX;
            colNr=2*p.nX:p.nX:p.nY*p.nX-p.nX;
            
            val1=1; val2=-1;
            ind=sub2ind(size(A),rowNr,colNr);
            A(rowNr,:)=0; b(rowNr)=0;
            A(ind)=val1;
            A(ind-NN)=val2;
            
            % Up/Down
            rowNr=2:p.nX-1;
            colNr=2:p.nX-1;
            val1=1; val2=-1;
            ind=sub2ind(size(A),rowNr,colNr);
            A(rowNr,:)=0; b(rowNr)=0;
            A(ind)=val1;
            A(ind+p.nX*NN)=val2;
            
            rowNr=(p.nY-1)*p.nX+2:NN-1;
            colNr=(p.nY-1)*p.nX+2:NN-1;
            val1=1; val2=-1;
            ind=sub2ind(size(A),rowNr,colNr);
            A(rowNr,:)=0; b(rowNr)=0;
            A(ind)=val1;
            A(ind-p.nX*NN)=val2;
            
            
            % corners
            val1=1;
            val2=p.dx/(p.dx+p.dy);
            val3=p.dy/(p.dx+p.dy);
            
            A(1,:)=0; b(1)=0;
            A(1,1)=val1;
            A(1,2)=-val3;
            A(1,p.nX+1)=-val2;
            
            A(p.nX,:)=0; b(p.nX)=0;
            A(p.nX, p.nX)=1;
            A(p.nX, p.nX-1)=-val3;
            A(p.nX, 2*p.nX)=-val2;
            
            A(NN,:)=0; b(NN)=0;
            A(NN,NN)=1;
            A(NN,NN-1)=-val3;
            A(NN,NN-p.nX)=-val2;
            
            A(NN-p.nX+1,:)=0; b(NN-p.nX+1)=0;
            A(NN-p.nX+1,NN-p.nX+1)=1;
            A(NN-p.nX+1,NN-p.nX+2)=-val3;
            A(NN-p.nX+1,NN-2*p.nX+1)=-val2;
            
        end
        
        A=sparse(A);
        
        % SOLVE
        sigConcLin_p1=A\b;
        
        sigConc_p1=reshape(sigConcLin_p1, p.nX,p.nY)';
      
end

end

%% ACCESSORY FUNCTIONS

function [binX, binY]=get_cellPosOnGrid(xPos, yPos, p)

%{
Function mapping cell position onto a rectangular grid

*Input*
xPos, yPos .. current cell positions
p .. parameter structure (contains information about grid)

*Output*
cellPosOnGrid .. array that contains the number of cells assigned to each
grid point
binX, binY .. describes which cells are assigned to which grid point


Written 12/2021 by Angelika Manhart
%}

%% find near which grid point the particles lies (Note: consider rewriting using weights, i.e. cells are assigned to several grid points with different weights)
        
        if p.optBC==0 % for periodic BCs
        
            % edges for histcounts such that grid points are the bin centers
            Xedges=[p.xVals p.Lx/2]-p.dx/2;
            Yedges=[p.yVals p.Ly/2]-p.dy/2;

            % shift those particles that would be "outside" the grid, due to the
            % periodicity
            xPosUse=xPos; xPosUse(xPosUse>max(p.xVals)+p.dx/2)=xPosUse(xPosUse>max(p.xVals)+p.dx/2)-p.Lx;
            yPosUse=yPos; yPosUse(yPosUse>max(p.yVals)+p.dy/2)=yPosUse(yPosUse>max(p.yVals)+p.dy/2)-p.Ly;

        elseif p.optBC==1 % for no-flux BCs
        
            Xedges=[p.xVals p.Lx/2+p.dx]-p.dx/2;
            Yedges=[p.yVals p.Ly/2+p.dy]-p.dy/2;
        
            xPosUse=xPos; yPosUse=yPos;
            
        end
        
        % assign particles to grid points
        [~, ~, ~, binX, binY]=histcounts2(xPosUse,yPosUse,Xedges, Yedges);
       
        
end