The problem of simulating the spectral line shapes of aligned immobile samples arises in solid-state NMR of various biological systems, including integral membrane proteins and peptides, receptor-bound ligands, and macroscopically oriented DNA fibers. An important issue with regard to the extraction of structural information is the correct treatment of the distribution of local symmetry axes relative to the average alignment axis (mosaic spread). Previous formulations have not considered explicitly the three-dimensional uniaxial character of the local axis disorder. Rather, the mosaic spread has been treated simply by convoluting the theoretical line shape function with an effectively two-dimensional distribution of the local symmetry axes. Here a closed-form line shape expression is derived for an axially symmetric distribution of bond orientations, which includes the uniaxial distribution of the local symmetry axis about the average alignment axis. As an illustration, the influences of the bond orientation and the degree of mosaic spread on deuterium (H-2) NMR line shapes are investigated. The closed-form solution in terms of elliptic integrals gives virtually identical results to those of an alternative numerical Monte Carlo line shape simulation method. The derived line shape function yields the correct powder-type limit, and has been tested by simulating a tilt series of H-2 NMR spectra of purple membranes containing bacteriorhodopsin with a specifically deuterated 1R methyl group in the retinal ring. The probability distribution for the bond orientations derived herein can be of potential interest for solid-state NMR spectroscopy of aligned biomolecules involving dipolar, quadrupolar, and chemical shift interactions, such as integral membrane proteins and peptides.