A new computer implementation of a method of modeling the interaction of naturally-occuring electromagnetic waves with three-dimensional Earth structures has been developed. The method is derived from the theory developed by Lord Rayleigh (1896) to describe the scattering of acoustic waves from rough surfaces. Specific application is made to magnetotelluric modeling. The Earth models considered are defined by a set of three-dimensional smoothly-varying interfaces that divide space into regions of constant electromagnetic properties. Adjacent interfaces may coincide over part or all of their areas, allowing the representation of closed bodies. A scattering solution is found by imposing boundary conditions derived from Maxwell's equations at points on the interfaces and electric field divergence conditions within the regions. The solution consists of the superposition of plane waves of varying orientations and amplitudes. The implementation has been tested by comparisons to an existing two-dimensional Rayleigh scattering program. The magnetotelluric response of a simple three-dimensional topographic feature is presented.