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=50;
Vp=1e-6;
%active force
A=50;
%Pressure gradient
dpdx=-10000;
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(-alphaT)) and solve for alphaT ("transformed").
%Subsequently, to obtain the correct alpha we need to retransform
%u13=alphaT.
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=(Vp*epsilon^2*6/pi)^(1/3);

%Mixture Reynolds number
ReM=ReF*(1-alphaS/alphaMax)^(5*alphaMax/2)*(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*Vp)^(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
wallNormalU=[0;1;0];
wallTangentialU=[1;0;0];
%Active force
active=alphaS*A*oriS2;

%Drag force
cs=1;
CD=pi*a*(1/(sqrtSpherPerp)+2*cs/(sqrtSpher))/(ReM);
drag=alphaS*CD*w; 

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

%Magnus force    
diffRot=0.5*shearRate-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*Vp*gradP;

%Sum of all forces  
allForces=(active+drag+saffman+magnus+lift+pressureForce)/Vp;
allForces=dot(allForces,wallNormalU)*(1-cos(y*pi)^2)*wallNormalU+dot(allForces,wallTangentialU)*wallTangentialU;
%%
%Torques
CR=pi*a^3*cs/(ReM*sqrtSpher^(3/2));
dragTorque=CR*(0.5*shearRate-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*Vp);

allTorques=(dragTorque+pitchingTorque)/Vp;

%%
%Diffusion coefficients
lambda=(Vp/alphaS)^(1/3);
absV=sqrt(u(4)^2+u(5)^2+u(6)^2);
nu=8*(alphaS)^2*epsilon^2*absV/(pi*a^4)*(epsilon+1);
muT=lambda^2*nu/3;
dmuTdy=8*epsilon^2*(epsilon+1)/(3*pi*a^4)*Vp^(2/3)*(4*alphaS^(1/3)*dalphaSdy*absV/3+alphaS^(4/3)*(u(4)*dudy(4)+u(5)*dudy(5)+u(6)*dudy(6))/absV);
            
%%
%Active gravitation
dalphaSKacdy=-D*(A*ReM*a/(32*pi*Vp^(2/3)))^2/2*2*alphaS*dalphaSdy;
lambdaStrich=a;
hu=(y)/lambdaStrich;
ho=(1-y)/lambdaStrich;
deltaH2=1/hu^4+1/ho^4;
deltaH4=1/hu^4+1/ho^4;
deltaH5=1/hu^5-1/ho^5;
deltaH6=1/hu^6+1/ho^6;
wallNormal=[0;1;0];
cosKappa=dot(oriS2,wallNormal)/eSCorr;
sinKappa=sqrt(1-cosKappa^2);
zeta=(epsilon^2-1)/(2*epsilon^2+1);
KW=-D*(3*A*ReM*a/(64*pi))^2/2*deltaH4*(1-3*cosKappa^2)^2;
dKWdy=-D*(3*A*ReM*a/(64*pi))^2/2*(12*(1-3*cosKappa^2)*cosKappa*dot(wallNormal,doriSdy)*deltaH4+4*(1-3*cosKappa^2)^2/a*deltaH5);

Wac=[ThetaA;ThetaB;ThetaC]/(8*lambda^6)*(A*ReM*a/(8*pi))^2*(3-5*zeta)^2;
WW=[ThetaA;ThetaB;ThetaC]*deltaH6/2*(3*A*a*ReM/(64*pi))^2*cosKappa^2*sinKappa^2*(1+zeta*(1+cosKappa^2))^2;

dKEdy=-A*a*ReM/(96*pi*lambda^2+1e-15)*5/3*alphaS*dalphaSdy;
KEW=3*A*a*ReM/(64*pi*deltaH2)*(3*cosKappa^2-1)*cosKappa;
dKEWdy=3*A*a*ReM/(64*pi)*(-2*lambdaStrich^2*cosKappa*(3*cosKappa^2-1)*(1/y^3+1/(1-y)^3)-sinKappa*dot(wallNormal,doriSdy)*(3*cosKappa^2-1)*(1/hu^2-1/ho^2)+cosKappa*(-6*cosKappa*sinKappa*dot(wallNormal,doriSdy)-1)*(1/hu^2-1/ho^2));
%%
%alphaT
%Convection
convAlphaT1=velS(2)*dudy(13);
convAlphaT2=(exp(u(13))+1)*dvelSdy(2);

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

%Diffusion
diffVelF=1/ReF*(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*(dKEdy+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