An interesting improvement would be to use the Schmidt semi-normalised Legendre functions as to avoid the use of factorials, which would make the output much more accurate for high L values.
just one question, let's say i want to find spherical harmonics with m=-1 and l=1, it gives me error statements. im sure it has to do with the modulus sign. if so, why do we put the modulus sign? thanx.