A mean-field model for the analysis of electrostatic interactions in liquid crystal phases is proposed. The model is based on a realistic description of the probe molecule and a continuum approximation for its environment. Electrostatic interactions, treated with the reaction field model, are superimposed to short-range interactions, which are accounted for by a shape method. The model is used to predict orientational order parameters of some representative solutes, characterized by different shapes and charge distributions, dissolved in nematics with different dielectric permittivities.