We consider the modal analysis of partially coherent submillimeter-wave quasi-optical systems. According to our scheme the cross-spectral density at a plane is expanded as a sum of partially coherent propagating free-space modes. The coherence matrix, the elements of which are determined by evaluating bimodal overlap integrals, completely describes the state of the field at a plane. The coherence matrix can be traced through an optical system to another plane by means of a scattering matrix. Whereas diagonalizing the scattering matrix gives the natural modes of the optical system, diagonalizing the coherence matrix gives the natural modes of the field. We demonstrate the technique by studying the behavior of Gaussian beam telescopes when the field at the source plane is completely incoherent.