We derive the Ornstein-Zernike equation for molecular crystals of axially symmetric particles and apply the Percus-Yevick approximation. The one-particle orientational distribution function rho((1)) (Omega) has a nontrivial dependence on the orientation Omega, in contrast to a liquid, and is needed as an input. Despite some differences, the Ornstein-Zernike equation for molecular crystals has a similar structure as for liquids. We solve both equations numerically for hard ellipsoids of revolution on a simple cubic lattice. Compared to molecular liquids, the orientational correlators in direct and reciprocal space exhibit less structure. However, depending on the lengths a and b of the rotation axis and the perpendicular axes of the ellipsoids, respectively, different behavior is found. For oblate and prolate ellipsoids with b greater, similar 0.35 (in units of the lattice constant), damped oscillations in distinct directions of direct space occur for some of the orientational correlators. They manifest themselves in some of the correlators in reciprocal space as a maximum at the Brillouin zone edge, accompanied by a maximum at the zone center for other correlators. The oscillations indicate alternating orientational fluctuations, while the maxima at the zone center originate from ferrorotational fluctuations. For a less, similar 2.5 and b less, similar 0.35, the oscillations are weaker, leading to no marked maxima at the Brillouin zone edge. For a greater, similar 3.0 and b less, similar 0.35, no oscillations occur any longer. For many of the orientational correlators in reciprocal space, an increase of a at fixed b or vice versa leads to a divergence at the zone center q=0, consistent with the formation of ferrorotational long-range fluctuations, and for some oblate and prolate systems with b less, similar 1.0 a simultaneous tendency to divergence of few other correlators at the zone edge is observed. Comparison of the orientational correlators with those from Monte Carlo simulations shows satisfactory agreement. From these simulations we also obtain a phase boundary in the a-b plane for order-disorder transitions.