  We develop a numerical method for approximating the surface modes of
sphere-like nanoparticles in the quasi-static limit, based on an expansion of
(the angular part of) the potentials into spherical harmonics. Comparisons of
the results obtained in this manner with exact solutions and with a
perturbation ansatz prove that the scheme is accurate if the shape deviations
from a sphere are not too large. The method allows fast calculations for large
numbers of particles, and thus to obtain statistics for nanoparticles with
random shape fluctuations. As an application we present some statistics for the
distribution of resonances, polariziabilities, and dipole axes for particles
with random perturbations.
