function [f, pf_pw] = part_pf(w, p)

% Returns rhs of particle equation given the current state (w)
% and the parameters of the problem (p)

% Unpack the parameters from p
rhog  = p(1);
a     = p(2);
mu    = p(3);
Omega = p(4);
Wg    = p(5);
mp    = p(6);

% Get particle velocity and position
up = w(1);
vp = w(2);
wp = w(3);
xp = w(4);
yp = w(5);
zp = w(6);

% Zero pf_pw

pf_pw = zeros(6,6);

% Calculate local gas velocity
ug = -Omega*yp;
vg =  Omega*xp;
wg =  Wg;

% Calculate relative velocity
urel      = ug - up;
purel_pyp = -Omega;
purel_pup = -1;

vrel      = vg - vp;
pvrel_pxp = Omega;
pvrel_pvp = -1;

wrel      = wg - wp;
pwrel_pwp = -1;

qrel = sqrt(urel*urel + vrel*vrel + wrel*wrel);

pqrel_purel = urel/qrel;
pqrel_pvrel = vrel/qrel;
pqrel_pwrel = wrel/qrel;

pqrel_pup = pqrel_purel*purel_pup;
pqrel_pvp = pqrel_pvrel*pvrel_pvp;
pqrel_pwp = pqrel_pwrel*pwrel_pwp;

pqrel_pxp = pqrel_pvrel*pvrel_pxp;
pqrel_pyp = pqrel_purel*purel_pyp;

% Calculate drag force
Re        = rhog*qrel*(2*a)/mu;
dRe_dqrel = rhog     *(2*a)/mu;

CD        = 24/Re + 6/(1+sqrt(Re)) + 0.4;
dCD_dRe   = -24*Re^(-2) - 6*(1+sqrt(Re))^(-2)*0.5/sqrt(Re);

Dx        = 0.5*rhog*qrel*pi*a^2*CD*urel;
pDx_pqrel = 0.5*rhog*     pi*a^2*CD*urel + 0.5*rhog*qrel*pi*a^2*dCD_dRe*dRe_dqrel*urel;
pDx_pup   = pDx_pqrel*pqrel_pup + 0.5*rhog*qrel*pi*a^2*CD*purel_pup;
pDx_pvp   = pDx_pqrel*pqrel_pvp;
pDx_pwp   = pDx_pqrel*pqrel_pwp;
pDx_pxp   = pDx_pqrel*pqrel_pxp;
pDx_pyp   = pDx_pqrel*pqrel_pyp + 0.5*rhog*qrel*pi*a^2*CD*purel_pyp;

Dy        = 0.5*rhog*qrel*pi*a^2*CD*vrel;
pDy_pqrel = 0.5*rhog*     pi*a^2*CD*vrel + 0.5*rhog*qrel*pi*a^2*dCD_dRe*dRe_dqrel*vrel;
pDy_pup   = pDy_pqrel*pqrel_pup;
pDy_pvp   = pDy_pqrel*pqrel_pvp + 0.5*rhog*qrel*pi*a^2*CD*pvrel_pvp;
pDy_pwp   = pDy_pqrel*pqrel_pwp;
pDy_pxp   = pDy_pqrel*pqrel_pxp + 0.5*rhog*qrel*pi*a^2*CD*pvrel_pxp;
pDy_pyp   = pDy_pqrel*pqrel_pyp;

Dz        = 0.5*rhog*qrel*pi*a^2*CD*wrel;
pDz_pqrel = 0.5*rhog*     pi*a^2*CD*wrel + 0.5*rhog*qrel*pi*a^2*dCD_dRe*dRe_dqrel*wrel;
pDz_pup   = pDz_pqrel*pqrel_pup;
pDz_pvp   = pDz_pqrel*pqrel_pvp;
pDz_pwp   = pDz_pqrel*pqrel_pwp + 0.5*rhog*qrel*pi*a^2*CD*pwrel_pwp;
pDz_pxp   = pDz_pqrel*pqrel_pxp;
pDz_pyp   = pDz_pqrel*pqrel_pyp;

% Calculate pressure force
Px      = -(4/3)*pi*a^3*rhog*Omega^2*xp;
pPx_pxp = -(4/3)*pi*a^3*rhog*Omega^2;
Py      = -(4/3)*pi*a^3*rhog*Omega^2*yp;
pPy_pyp = -(4/3)*pi*a^3*rhog*Omega^2;
Pz      = 0;

% Sum forces
Fx = Dx + Px;
Fy = Dy + Py;
Fz = Dz + Pz;

% Calculate right-hand sides
f = zeros(6,1);
f(1) = Fx/mp;
f(2) = Fy/mp;
f(3) = Fz/mp;
f(4) = up;
f(5) = vp;
f(6) = wp;

% Assemble derivatives
pf_pw(1,1) = (pDx_pup          )/mp;
pf_pw(1,2) = (pDx_pvp          )/mp;
pf_pw(1,3) = (pDx_pwp          )/mp;
pf_pw(1,4) = (pDx_pxp + pPx_pxp)/mp;
pf_pw(1,5) = (pDx_pyp          )/mp;

pf_pw(2,1) = (pDy_pup          )/mp;
pf_pw(2,2) = (pDy_pvp          )/mp;
pf_pw(2,3) = (pDy_pwp          )/mp;
pf_pw(2,4) = (pDy_pxp          )/mp;
pf_pw(2,5) = (pDy_pyp + pPy_pyp)/mp;

pf_pw(3,1) = (pDz_pup          )/mp;
pf_pw(3,2) = (pDz_pvp          )/mp;
pf_pw(3,3) = (pDz_pwp          )/mp;
pf_pw(3,4) = (pDz_pxp          )/mp;
pf_pw(3,5) = (pDz_pyp          )/mp;

pf_pw(4,1) = 1.0;
pf_pw(5,2) = 1.0;
pf_pw(6,3) = 1.0;