function [xPos_p1, yPos_p1, cellDir_p1, cellSpeed_p1, chemSens_p1]=...
    get_cellMove_p1(xPos, yPos, cellDir, cellSpeed, chemSens, sigConc, sigConcPerCell, repForce, p)
%{
This function updates the cell position, direction, speed and chemotacic
sensitivity.
    
*Input*
Current positions, directions, speed and sensitivity:
- xPos, yPos, cellDir, cellSpeed, chemSens    
- sigConc .. signal concentration on grid
- sigConcPerCell .. signal concentration per cell
- repForce .. repulsive force     
- p .. parameters structure        
    
*Output*
positions, directions, speed and sensitivity as the next time step (p1 stands for "+1")    
        
Written 12/2021 by Angelika Manhart   
%}
    
    
% gradient of signal concentration
dxSigConc=(circshift(sigConc,-1,2)-circshift(sigConc,1,2))/(2*p.dx);
dySigConc=(circshift(sigConc,-1,1)-circshift(sigConc,1,1))/(2*p.dy);

[x,y] = meshgrid(p.xVals,p.yVals);

% update grid for interpolation

if p.optBC==0 % period BCs
    
    % account for periodicity
    dxSigConc2=[dxSigConc dxSigConc(:,end)];
    dxSigConc2=[dxSigConc2; dxSigConc2(end,:)];
    
    dySigConc2=[dySigConc dySigConc(:,end)];
    dySigConc2=[dySigConc2; dySigConc2(end,:)];
    
    x2=[x x(:,end)+p.dx];
    x2=[x2; x2(end,:)];
    
    y2=[y y(:,end)];
    y2=[y2; y2(end,:)+p.dy];
    
elseif p.optBC==1 % no flux BC
    
    dxSigConc(:,1)=0; dxSigConc(:,end)=0;
    dySigConc(1,:)=0; dySigConc(end,:)=0;
    
    dxSigConc2=dxSigConc; dySigConc2=dySigConc;
    x2=x; y2=y;
end
 
% interpolate
dxSigConcPerCell=interp2(x2,y2,dxSigConc2, xPos, yPos);
dySigConcPerCell=interp2(x2,y2,dySigConc2, xPos, yPos);

% calculate norm of gradient and direction
normGrad=sqrt(dxSigConcPerCell.^2+dySigConcPerCell.^2);
dirGrad=mod(atan2(dySigConcPerCell,dxSigConcPerCell),2*pi);

%% Update cell speed, direction and sensitivity according to equations
%{
Equations are given in: RHS_modChemSens, RHS_modCellSpeed, RHS_modCellDir
Here they are solved on [0 dt].
This could be replaced by a one-step solver to make faster
%}

% Update chemotactic sensitivity

y0=chemSens';
tspan=[0; p.dt];
                
[~,y] = ode45(@(t,y)RHS_modChemSens(t,y,sigConcPerCell,p),tspan,y0);
chemSens_p1=y(end,:)';

% Update cell speed

y0=cellSpeed';
tspan=[0; p.dt];
                
[~,y] = ode45(@(t,y)RHS_modCellSpeed(t,y,chemSens, normGrad,p),tspan,y0);
cellSpeed_p1=y(end,:)';


% Update cell direction

y0=cellDir';
tspan=[0; p.dt];
                
[~,y] = ode45(@(t,y)RHS_modCellDir(t,y,chemSens, normGrad, dirGrad,p),tspan,y0);
cellDir_p1=y(end,:)';


%% Update cell position

xPos_p1=xPos+p.dt*repForce(:,1)+p.dt*cellSpeed.*cos(cellDir)+sqrt(2*p.movNoise*p.dt)*randn(p.N,1);
yPos_p1=yPos+p.dt*repForce(:,2)+p.dt*cellSpeed.*sin(cellDir)+sqrt(2*p.movNoise*p.dt)*randn(p.N,1);
    
end

%% ACCESSORY FUNCTIONS

function dydt=RHS_modChemSens(t,y,sigConcPerCell, p)

N=length(sigConcPerCell);

chemSens=y(1:N);

dChem=(-p.chemDesens*sigConcPerCell.*chemSens+1-chemSens)/p.chemTimeScale;

dydt=dChem;

end

function dydt=RHS_modCellSpeed(t,y,chemSens, normGrad, p)

N=length(chemSens);

cellSpeed=y(1:N);

dSpeed=p.speedModelP1*chemSens.*normGrad-p.speedModelP2*cellSpeed;

dydt=dSpeed;

end

function dydt=RHS_modCellDir(t,y,chemSens, normGrad, dirGrad, p)

N=length(chemSens);

cellDir=y(1:N);

dDir=p.turnRate*chemSens.*normGrad.*sin(dirGrad-cellDir);

dydt=dDir;

end
    