Analytic functions with blackbody

See #3077 
Definitely do not mask very small `x` -- this means that the Rayleigh-Jeans limit will be wrong!! I also do not see the reason for `_VERY_LARGE`: numpy does the right thing already!
```
In [7]: with np.errstate(all='ignore'):
    print(np.expm1(1000.))
    print(1./np.expm1(1000.))
   ...:     
inf
0.0
```