Analytic functions with blackbody

See #3077 
@mhvk , I see on the original PR comments that people were generally against automatic unit assignment. So, the middle ground here that I can think of is to set `flux_unit=None` by default, which just returns SI unit, *and* then write two wrappers (`planck_lambda` and `planck_nu`). How does this sound?