Starting from WIMP density profiles, in the framework of the Eddington approach, we obtain the energy distribution f ( E ) of dark matter in our vicinity. Assuming a factorizable phase space function, f ( E , L ) = F ( E ) F L ( L ) , we obtain the velocity dispersions and the anisotropy parameter β in terms of the parameters describing the angular momentum dependence. By employing the derived expression f ( E ) we construct axially symmetric WIMP velocity distributions. The obtained distributions automatically have a velocity upper bound, as a consequence of the fact that they are associated with a gravitationally bound system, and are characterized by an anisotropy parameter β . We then show how such velocity distributions can be used in determining the event rates, including modulation, both in the standard as well directional WIMP searches.