  A method for 3D interpolation between hard spheres is described. The function
to be interpolated could be the charge density between atoms in condensed
matter. Its electrostatic potential is found analytically, and so are various
integrals. Periodicity is not required. The interpolation functions are
localized structure-adapted linear combinations of spherical waves, socalled
unitary spherical waves (USWs), centered at the spheres where they have
cubic-harmonic character. Input to the interpolation are the coefficients in
the cubic-harmonic expansion of the target function at and slightly outside the
spheres; specifically, the values and 3 first radial derivatives. To fit this,
we use USWs with 4 negative energies. Each interpolation function is actually a
linear combination of these 4 sets of USWs with the following properties: (1)
It is centered at a specific sphere where it has a specific cubic-harmonic
character and radial derivative. (2) Its value and first 3 radial derivatives
vanish at all other spheres and for all other cubic-harmonics. It is therefore
highly localized. Explicit expressions are given for these value-and-derivative
(v&d) functions and their Coulomb potentials in terms of the USWs and their
structure matrix, including those with zero energy, as well as for integrals
over the interstitial of the v&d functions and of their products. Use of point-
and space-group symmetries can significantly reduce matrix sizes and the number
of v&d functions. For open structures, additional constraints are needed to
pinpoint the interpolated function deep in the interstitial. This requires the
USW set at one extra energy. As examples, we consider a constant density and
the valence-electron densities in zinc-blende structured Si, ZnSe, and CuBr.
