There is today a growing need to accurately model the angular scattering response of metasurfaces for optical analog processing applications. However, the current metasurface modeling techniques are not well suited for such a task since they are limited to small angular spectrum transformations, as shall be demonstrated. The goal of this work is to overcome this limitation by improving the modeling accuracy of these techniques and, specifically, to provide a better description of the angular response of metasurfaces. This is achieved by extending the current methods, which are restricted to dipolar responses and weak spatially dispersive effects, so as to include quadrupolar responses and higher-order spatially dispersive components. The accuracy of the newly derived multipolar model is demonstrated by predicting the angular scattering of a dielectric metasurface placed in a vacuum. This results in a modeling accuracy that is substantially better than the standard dipolar model.