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. Input to the interpolation are the coefficients in the cubic-harmonic (lm) expansion of the target function at and slightly outside the spheres; specifically, the values and the 3 first radial derivatives, d. Hence, we use USWs with 4 different non-positive energies. Each interpolation function is actually a linear combination of 4 sets of USWs with the following properties: It is centered at a specific sphere, R, where it has a specific lm and d, and it vanishes at all other spheres, or for all other l, m, or d. It is therefore highly localized. This set of so-called value-and-derivative (v&d) functions depends on the sites R and radii of the hard spheres and can be generated by inversion of 4 structure matrices of dimension N_{R}\(\times\)N_{lm} where N_{R} is of order 100. Explicit expressions are given for the v&d functions and their Coulomb potentials in terms of the USWs, including those with zero energy, as well as for integrals over the interstitial of the v&d functions and of their products. 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 in which cases N_{lm} is merely 4.