function [c,f,s] = AVShearPDE(y,t,u,dudy)
%This function contains parameters and equations for a quasi-2D channel
%flow. The particles are three dimensional, however, any derivative in the
%z-direction (perpendicular to the image plane) is neglected.

%Parameters (can be changed freely)
%Maximum solid volume fraction
alphaMax=0.74;
%density ratio
D=1;
%Reynolds number
ReF=1e-3;
%Particle geometry
epsilon=16;%aspect ratio
volumeP=1e-5;
%active force
activeF=500;
%Pressure gradient
dpdx=-1e6;
dpdy=0;
dpdz=0;
gradP=[dpdx;dpdy;dpdz];

%Do not change anything below this comment!
%volume fraction
%The 13th variable in the system is not the solid phase volume fraction.
%Instead we apply a transformation with a logistic function
%alphaS=alphaMax/(1+exp(-gAlpha)) and solve for gAlpha.
%Subsequently, to obtain the correct alphaS we need to retransform
%u13=gAlpha.
alphaS=alphaMax/(1+exp(-u(13)));
dalphaSdy=alphaMax*exp(-u(13))*dudy(13)/(1+exp(-u(13)))^2;
alphaF=1-alphaS;
dalphaFdy=-dalphaSdy;

%Particle length
a=(volumeP*epsilon^2*6/pi)^(1/3);

%Channel height
H=1;

%Mixture Reynolds number
muInt=2.5+0.123*(epsilon-1)^(0.925);
ReM=ReF*(1-alphaS/alphaMax)^(muInt*alphaMax)*(1-heaviside(alphaS-alphaMax));

%Fluid velocity
velF=[u(1);u(2);u(3)];
dvelFdy=[dudy(1);dudy(2);dudy(3)];

%Solid translational velocity
velS=[u(4);u(5);u(6)];
dvelSdy=[dudy(4);dudy(5);dudy(6)];

%Relative velocity
w=[u(1)-u(4);u(2)-u(5);u(3)-u(6)];
absW=sqrt(w(1)^2+w(2)^2+w(3)^2);

%Solid rotational velocity
rotS=[u(7);u(8);u(9)];
drotSdy=[dudy(7);dudy(8);dudy(9)];

%Orientation
oriS=[u(10);u(11);u(12)];
doriSdy=[dudy(10);dudy(11);dudy(12)];
eS=sqrt(u(10)*u(10)+u(11)*u(11)+u(12)*u(12))+1e-15;
eSCorr=sqrt(1/(1+exp(-u(14))));
oriS2=eSCorr*oriS/eS;
%%
%Sphericity
Ap=4*pi*((2*(a^2/epsilon)^(8/5)+(a^2/epsilon^2)^(8/5))/3)^(5/8);
Sphericity=pi^(1/3)*(6*volumeP)^(2/3)/Ap;
sqrtSpher=sqrt(Sphericity);
cosAngleIncidence=u(10)*w(1)+u(11)*w(2)+u(12)*w(3)/(absW*eSCorr+1e-15);
%Cross sectional spericity
aPerp=sqrt((a*cosAngleIncidence)^2+(a*sqrt(1-cosAngleIncidence^2)/epsilon)^2);
sqrtSpherPerp=sqrt(a^2/(epsilon*aPerp^2));

%%
%Moment of inertia
%Principle axes
ThetaA=D/5*(2*(a/epsilon)^2);
ThetaB=D/5*(a^2+(a/epsilon)^2);
ThetaC=D/5*(a^2+(a/epsilon)^2);
Theta=[ThetaA 0 0; 0 ThetaB 0; 0 0 ThetaC];

%%
%Forces
%Active force
active=alphaS*activeF*oriS2;

%Drag force
cD_Ac=5/6;%correction factor due to active boundary condition
cD_H=1+2.33*epsilon^(2/5)*a/(epsilon*H);%correction factor due to the confinement
CD=pi*a*cD_Ac*cD_H*(1/(sqrtSpherPerp)+2/(sqrtSpher))/(ReM);
drag=alphaS*CD*w; 

%Saffman force
rotF=[dudy(3);0;-dudy(1)];
absRotF=sqrt(rotF(1)^2+rotF(2)^2+rotF(3)^2);
saffman=1.615*alphaS*sqrt(a/(ReM+1e-15))*(absRotF)^(-0.5)*cross(w,rotF);

%Magnus force    
diffRot=0.5*rotF-rotS;
magnus=pi*a^3*alphaS*cross(diffRot,w)/8;

%Circulation lift force (airfoil)
cE=0.5*cos(2*eSCorr*pi/sqrt(2))*(1-heaviside(eSCorr-sqrt(2)/2))+0.5+0.5*cos(2*pi*eSCorr/(2-sqrt(2))-2*pi/(2-sqrt(2)))*heaviside(eSCorr-sqrt(2)/2);
CL=CD*0.92*sin(0.647*acos(cosAngleIncidence))*sin(0.647*acos(cosAngleIncidence)+0.351*pi)*cos(0.647*acos(cosAngleIncidence)+0.177*pi);
lift=alphaS*CL*sign(dot(w,oriS2))*cE*pi*a^2*cross(cross(w,oriS2),w)/8;

%Pressure force
pressureForce=-alphaS*volumeP*gradP;

%Sum of all forces  
allForces=(active+drag+saffman+magnus+lift+pressureForce)/volumeP;
%%
%Torques
cDT=0.5;%Correction factor due to the active boundary condition
CR=pi*a^3*cDT/(ReM*sqrtSpher^(3/2));
dragTorque=CR*(0.5*rotF-rotS)*alphaS;

%Pitching torque, caused by the hydrodynamic forces
centerOfPressure=oriS2*sign(dot(w,oriS2))*(1-(sqrt(1-cosAngleIncidence^2))^3)*(1-exp(1-epsilon))*a/4;
pitchingTorque=cross(centerOfPressure,allForces*volumeP);

allTorques=(dragTorque+pitchingTorque)/volumeP;

%%
%Diffusion coefficient
muT=1/ReM;    
dmuTdy=-muInt*dalphaSdy*(1-alphaS/alphaMax)^(alphaMax*muInt-1);
%%
%Active attraction
lambda=(volumeP/alphaS)^(1/3);%averaged distance between two particles
cKac=(a/2)^2;
Kac=11*D/((lambda^2+cKac)^2)*(activeF*ReM*a/(32*pi))^2;
dKacdy=11*D/((lambda^2+cKac)^2)*(activeF*a/(32*pi))^2*dalphaSdy*(ReM^2*4*lambda^5/(3*volumeP*(cKac+lambda^2))-2*ReM*muInt*(1-alphaS/alphaMax)^(muInt*alphaMax-1));
dalphaSKacdy=Kac*dalphaSdy+alphaS*dKacdy;
hu=(y)/a;
ho=(1-y)/a;
deltaH1=1/(hu^2+cKac)+1/(hu^2+cKac);
deltaH2=1/(hu^2+cKac)^2+1/(hu^2+cKac)^2;
deltaH3=1/(hu^2+cKac)^3+1/(hu^2+cKac)^3;
deltaH6=1/(hu^3+cKac^(3/2))^2+1/(ho^3+cKac^(3/2))^2;
wallNormal=[0;1;0];
cosKappa=dot(oriS2,wallNormal)/eSCorr;
sinKappa=sqrt(1-cosKappa^2);
zeta=(epsilon^2-1)/(epsilon^2+1);

KW=-D*(3*activeF*ReM*a/(64*pi))^2/2*deltaH2*(1-3*cosKappa^2)^2;
dKWdy=-D/2*(3*activeF*a/(64*pi))^2*(muInt*dalphaSdy*deltaH2*(1-alphaS/alphaMax)^(muInt*alphaMax-1)*2*ReM*(1-3*cosKappa^2)^2+ReM^2*(hu+ho)/a*deltaH3*(1-3*cosKappa^2)^2+ReM^2*deltaH2*2*(1-3*cosKappa^2)*(6*cosKappa)*dot(wallNormal,doriSdy));

Wac=[ThetaA;ThetaB;ThetaC]/((lambda^3+cKac^(3/2))^2)*(activeF*ReM*a/(32*pi))^2*(3-4*zeta)^2;
WW=[ThetaA;ThetaB;ThetaC]*deltaH6/2*(3*activeF*a*ReM/(64*pi))^2*cosKappa^2*sinKappa^2*(1+zeta*(1+cosKappa^2)/2)^2;

KEW=3*activeF*a*ReM/(64*pi)*deltaH1*(3*cosKappa^2-1)*cosKappa;
dKEWdy=3*activeF*a/(64*pi)*(-muInt*dalphaSdy*(1-alphaS/alphaMax)^(muInt*alphaMax-1)*(3*cosKappa^2-1)*cosKappa*deltaH1-ReM*2*(hu+ho)*a*deltaH2*(3*cosKappa^2-1)*cosKappa-ReM*deltaH1*(dot(wallNormal,doriSdy)*(3*cosKappa^2-1)+cosKappa*6*cosKappa*dot(wallNormal,doriSdy)));
%%
%gAlpha
%Convection
convAlphaT1=velS(2)*dudy(13);
convAlphaT2=(exp(u(13))+1)*dvelSdy(2);

%%
%Fluid velocity
%Convection
convVelF=alphaF*velF(2)*dvelFdy;

%Diffusion
diffVelF=1/ReM*(alphaF*dvelFdy+velF*dalphaFdy);

%Production
prodVelF1=alphaF*gradP;
prodVelF2=-2*dalphaSKacdy/3;

%%
%Solid translational velocity
%Convection
convVelS=D*alphaS*velS(2)*dvelSdy;

%Diffusion
diffVelS=D*muT*alphaS*dvelSdy;

%Production due to active stress
prodVelS=2*dalphaSKacdy/3+(2*KW*dalphaSdy/3-2*alphaS*dKWdy);

%%
%Solid rotational velocity
%Convection
convRotS=Theta*alphaS*velS(2)*drotSdy;

%Diffusion
diffRotS=Theta*muT*alphaS*drotSdy;

%Production
prodRotS1=alphaS*cross(rotS,Theta*rotS);

%Production due to active stress
prodRotS2=2*(WW+Wac)*alphaS/3;
%%
%Orientation
%Convection
convOriS=alphaS*velS(2)*doriSdy;

%Diffusion
muE=muT;
diffOriS=muE*alphaS*doriSdy;

%Production
prodOriS1=alphaS*cross(rotS,oriS);

%Production due to collisions
prodOriS2=dmuTdy*doriSdy(1);
prodOriS3=2/3*(KEW*dalphaSdy+alphaS*dKEWdy);

%Orientation vector length
%Convection
convOriLength1=alphaS*exp(u(14))/(2*(1+exp(u(14)))^2)*velS(2)*dudy(14);

%Diffusion
diffOriLength=alphaS*muE*(dot(oriS,doriSdy)+oriS(2)*doriSdy(2));
notDiffOriLength=alphaS*muE*(dot(doriSdy,doriSdy)+doriSdy(2)*doriSdy(2));
%%
%System
%Temporal derivatives
c=[alphaF;alphaF;alphaF;D*alphaS;D*alphaS;D*alphaS;Theta(1,1)*alphaS;Theta(2,2)*alphaS;Theta(3,3)*alphaS;alphaS;alphaS;alphaS;1;alphaS*exp(u(14))/(2*(1+exp(u(14)))^2)];

%Diffusion
f=[diffVelF(1);
   2*diffVelF(2);
   0;
   diffVelS(1);
   2*diffVelS(2);
   0;
   0;
   0;
   diffRotS(3);
   diffOriS(1);
   2*diffOriS(2);
   0;
   0;
   diffOriLength
   ];

%Convection and production
s=[-allForces(1)-convVelF(1)-prodVelF1(1);
   -allForces(2)-convVelF(2)-prodVelF1(2)-prodVelF2;
    0;
    allForces(1)-convVelS(1);
    allForces(2)-convVelS(2)-prodVelS;
    0;
    0;
    0;
    allTorques(3)-convRotS(3)-prodRotS1(3)-prodRotS2(3);
    prodOriS1(1)-convOriS(1)+prodOriS2;
    prodOriS1(2)-convOriS(2)-prodOriS3;
    0;
   -convAlphaT1-convAlphaT2;
   -u(11)*prodOriS3+u(10)*prodOriS2-convOriLength1-notDiffOriLength
   ];
end

