Electroencephalography (EEG) is a key modality to monitor brain activity with high temporal resolution. EEG makes use of an array of electrodes to measure the electrical potential on the scalp. While most traditional EEG analyses have looked at EEG rhythms in different frequency bands, another important application of EEG is source imaging; i.e., map back the measured scalp potential to the underlying source distribution. However, source localization is an ill-posed problem. Most approaches consider a grid of neurobiologically relevant dipoles with fixed positions and directions of their moments, which makes for many more unknown intensities than the number of measures. To make the solution unique, one needs to add regularization strategies such as smoothness or sparsity. The alternative, which is the approach that we will follow in this work, is to impose a sparse and parametric source model with few unknown parameters, but which include both dipole positions and moments thus rendering the problem solution highly non-linear. The proposed method, for which we coin the term "analytic sensing", is based on two main working principles. First, based on the divergence theorem, the sensing principle relates the boundary potential to a volumetric information about the sources. This principle makes use of a mathematical test function ("analytic sensor") that needs to be a homogeneous solution of the Poisson equation, which is the governing equation of the quasi-static electromagnetic setting. The analytic sensor can be defined for different geometries and conductivity profiles of the domain of interest. We derive closed-form expressions for 3D multi-layer spherical models that are often used as a head model in EEG. Recent advances in signal processing known as "annihilation filter" or "finite rate of innovation", have extended Prony's method to recover sparse sources in a robust way, typically a stream of Diracs. We propose to apply the annihilation principle by imposing a particular choice of multiple analytic sensors; i.e., the virtual measurements obtained with these functions can be annihilated by a filter that allows us to recover the dipoles' positions in a non-iterative way. The dipoles' moments can subsequently be retrieved by solving a linear system of equations. While the application of the sensing principle is intrinsically 2D or 3D, the annihilation principle gives only access to the orthographic projection of the source distribution. Two approaches have been developed. First, using coordinate transformations, one can obtain multiple projections of the sources and recombine them to reconstruct the full 3D information. Second, using a second set of analytic sensors, it is possible to retrieve the missing coordinate by solving another linear system of equations. Successful application of the method requires careful implementation of both principles. For the sensing principle, we project the boundary measurements on a set of spherical harmonics with proper regularization to cover parts of the boundary that are not measured. We also show how the annihilation step can be implemented to be robust to noise and numerically stable. We demonstrate the precision and robustness of the method by both experimental results and theoretical Cramèr-Rao bounds. Finally, we show effective source localization for real-world experimental EEG data; i.e., we identify the underlying sources for several time instants of a visual evoked potential.