I have a three dimensional real function that I am attempting to represent as a sum of real spherical harmonics ($Y_{l,m}$ instead of $Y_l^m$) multiplied by the distance from origin raised to some power. Essentially, I am trying to write $$f(r,\theta,\phi) = \sum_l \sum_m f_{l,m} \times r^lY_{l,m}(\theta,\phi)$$ where $$f_{l,m} = \int_\Omega d\Omega ~ f'(\theta,\phi) Y_{l,m}(\theta,\phi)$$ with $f'(\theta,\phi) = f(R,\theta,\phi)\frac{1}{C}$ and $C$ is a scalar that I'm choosing (i.e. I evaluate $f$ at some constant radius $R$ and scale it by some factor to account for the radius actually being nonunity).
A test that I'm trying to running is fitting $f = z^2$ to this form. I get agreement on the sphere of radius that I evaluate this function on, however off of that sphere I am getting a really poor fit in 3D (e.g. a DC value shows up).
In retrospect, this is sensible to me, since the Spherical Harmonics only form a basis for the unit sphere. Additionally, since $z^2$ (and the functions I actually will want to fit to this form) doesn't scale uniformly under varying the radius, the coefficients will depend heavily on the radius that I evaluate $f$ at. Finally perhaps the most obvious reason against this approach just by inspection is that you can't construct strictly $z^2/r^2$ out of the spherical harmonics.
My question is twofold: first, is my intuition at why this is a poor idea/not working correct? And second, is there a nice way to salvage the general idea of what I'm going for? I really would like to have terms of the form $x^iy^jz^k$ to project onto for the work that follows this step. The alternative is to just take $x^iy^jz^k$ or the spherical harmonics and do a LLS fit to the data, but I would like to avoid that if possible