Objects such as spacecraft and large particles in the coma of an active comet are perturbed by gas drag. These perturbations are not currently well understood. Previous work by Moretto and McMahon (2020) explored the effects of the radial component of coma drag on the mean elements by expanding that force as a Fourier series in argument of latitude and averaging over Gauss' planetary equations. This takes advantage of the fact that radial perturbations do not affect the orbit plane and that the coma density varies as the inverse of the square of the radial distance. For a more general approach, the drag force from the coma can be represented using spherical harmonic expansions for the radial, in-track, and cross-track components of the perturbing acceleration. Using Kaula's inclination function (Kaula, 1960), the spherical harmonics expansions can be reduced to a Fourier series in the true anomaly and averaged using the same procedure. The resulting differential equations yield analytical insights into the evolution of an orbit over many orbit periods and can be applied to mission design or the motion of natural grains. The derivation of the differential equations for the mean elements is shown and numerically validated. Application to comets is shown by expanding a simplified flat plate drag model using spherical harmonics.