I tried to rotate zonal spherical function, projected onto spherical harmonic basis. Something goes wrong. Then I tried to rotate just one spherical harmonics zonal basis function and I also can't do it, except for band 0 (trivial case, because it is just a sphere) and for band 1.

```
sage: def spherical_point(phi, z, f):
radius = f.subs(phi=phi, z=z)
color = 'blue'
if radius < 0:
color = 'red'
radius = -radius
sine = sqrt(1 - z * z)
return point3d((cos(phi) * sine * radius, sin(phi) * sine * radius, z * radius), color=color, size=1)
sage: var('z, phi', domain='real')
sage: step = 1.0 / 32
sage: band = 2; dir_phi = pi / 4; dir_theta = pi / 4
sage: sh = sum(spherical_harmonic(band, m, dir_theta, dir_phi).real() * spherical_harmonic(band, m, arccos(z), phi).real() * sqrt(4 * pi / (2 * band + 1)) for m in range(-band, band + 1))
sage: C = sum(spherical_point(phi, z, sh) for phi in [0..2*pi, step=2*pi*step] for z in [-1..1, step=step])
sage: C.show(apect_ratio=1, viewer='threejs', axes=True)
```

If `dir_phi`

and `dir_theta`

differs from zero, then I only see a deformed zonal spherical harmonics of the band `band`

.

The only I can suspect is that `spherical_harmonic`

are not normalized properly. Is it true?