Modelling and optimization of non-linear multilayer all-optical devices including a grating coupler yield numerous prerequisite calculations in the linear regime. Indeed, it is of prime importance to know the influence of various physical and geometrical parameters (thickness and refractive index of each layer; grating groove depth, period and shape; angular spectrum of the input beam) on the diffraction properties. In this paper, it is shown how the well known (but often not well used) analytical Rayleigh-Fourier (RF) method can be used in many cases of practical interest, well above the validity limit generally admitted for a grating of sinusoidal shape (i.e. up to at least a groove depth equal to the grating period). Numerical problems due to the well known 'exponential overflow' due to evanescent waves are solved and so, up to 25 diffraction orders can be currently taken into account (up to 31 in some cases), whatever the number and thickness of substrate layers lying under the guiding grooved layer are. The use of equivalence rules for non-sinusoidal gratings is described. As a straightforward application of the plane wave calculation, the case of a finite beam is treated by taking advantage of its angular spectrum (a fast Fourier transform is used). The theoretical results obtained using the RF method are checked using the rigorous coupled wave method of Chateau and Hugonin. Finally, on a test example, dealing with a limited beam, we show that the diffraction efficiencies, calculated by the method we have implemented, are in perfect coincidence with those obtained experimentally.