@skendrew - just for information, if anyone requests any changes during review, you just need to commit the changes and push to the same branch, and the pull request will update.

Quick comment for now - it looks like the changes in functional_models.py include tabs instead of spaces - could you replace any tab by four spaces?

OK thanks and will do - to both comments. I haven't configured my editor properly yet on my new laptop!!

@kendrew -- Would it be possible to have an option to return f_nu units instead of f_lambda units for the mid-IR through radio users out there? It's also easier to convert f_nu to AB mags.

@aconley - I was going to bring up the same point - I think we can just have two separate functions, planck_nu and planck_lambda. @skendrew, what do you think?

@astrofrog, @skendrew -- on second thought, maybe that's just cluttering up the interface. After all, there is astropy.units.equivalencies.spectral_density. But it would be nice if this were called plank_lambda so users aren't always having to read the docs to remember which one it is.

@skendrew - Thanks for the contribution! I've finally had a chance to review this, and I've left a few comments above - they are all fairly minor! Some are just a question of style or consistency with the rest of Astropy.

The way we usually deal with comments like this is to implement each of them (sometimes one commit per main comment), push the new commits to the branch, then respond to the comments to say if a particular comment has been addressed (even just saying 'Done'). You may have to close on 'show outdated diff' in some cases as the comment may disappear if it is no longer relevant to the code in the pull request.

@aconley - agreed, let's just stick to planck_lambda for now, and we can always add planck_nu later if we want.

@skendrew - minor comments:

1. Would it not be nicer to just call it black_body (or `blackbody)? I have somewhat of an aversion to calling things after people, as it obscures meaning (my favourite example is "Brunt-Väisälä frequency" versus "boyancy frequency"). In this particular case, I think it also is better since the usual symbol for it is B(T).
2. I realise it doesn't really matter, but why is the default output unit in W/m2/micron? (see also below)

@aconley, @astrofrog: I'm not sure I like the suggestion of having *_lambda. With equivalencies, it is easy enough to convert as it is. One suggestion, which would also solve (2) above, is to have the flux always be "per input wavelength/frequency/energy unit", i.e.,

w = 1000. * u.nm
f = 1.4 * u.GHz
e = 1. * u.keV
T = 10000. * u.K
black_body(w, T) # gives F_lambda [W/m2/nm]
black_body(f, T) # gives F_nu [W/m2/GHz]
black_body(e, T) # gives F_E [W/m2/keV]  -> used a lot in X-ray astronomy
black_body(w.to(u.Hz, equivalencies=u.spectral()), T) # F_nu [W/m2/Hz]
black_body(f.to(u.Hz), T) # F_nu [W/m2/Hz]

One advantage of this is that the function will not return frequency units for a wavelength (i.e., length units), which
is not all that physical (e.g., I think it is nicer if one can just do w*black_body(w, T), e*black_body(e,T), and f*black_body(f, T), and always get the correct unit (W/m2). A small disadvantage is that you need to also implement the frequency and energy versions (or do the conversions, though that will be slower).

p.s. As an implementation of the above, one could have the individual black_body_lambda, black_body_frequency and black_body_energy, but those would not convert frequency to wavelength, etc., and then add an upper level black_body function that calls the appropriate one based on the units of w. This has the advantage that if you just want to avoid all the checking, you can use the lower-level function.

@mhvk - Calling it blackbody is fine -- but I can't say I'm worried about anybody likely to be using astropy not immediately knowing the association between Max Planck and the blackbody formula!

For my own personal convenience, I would prefer whatever_lambda and whatever_nu, because the standard in the mid-IR through mm is to work in f_nu as a function of wavelength, and so anybody working there will always have to mucking about with equivalencies. And odd convention, I know, but such is life. At least it isn't as bad as antennae temperature!

I'm with @aconley - I don't think the unit of the input variable should determine whether it's f_nu or f_lambda, and in the mid-IR it is indeed common to use f_nu(lambda).

My suggestion for now is to first implement comments about planck_lambda itself, then once the function is ready, it's easy to simply add a second one called planck_nu and then just adjust the formula. For the purpose of this pull request, let's maybe not worry too much about the underflow/overflow, and we can fix it in a subsequent pull request?

@astrofrog, @aconley - I thought a bit more and agree that the individual functions are useful to have (also for speed-up in modelling), so good to go ahead. As long as we have a *_lambda name, we can always add an overall function later. (And neat that numpy has np.expm1 built-in; fine not to worry about the overflow - it does give the right answer.)

As for the name, yes, of course, most astronomers will understand planck likely is a black body, but I think black_body is even clearer, and probably the first thing one would search for if one wants to know whether astropy includes it.

As a final question: on second thought I also wondered whether one should convert the flux to some specific units at all, whether [W/m2/mu] or anything else; the beauty of quantities is that they represent the same physical quantity (in this case flux with wavelength units) no matter what the precise units. So, perhaps just remove

flux = flux.to(u.Watt / (u.m * u.m * u.um * u.sr)) 

For modelling, this would make it a little faster, as it is unlikely the data being modelled are in exactly those units, so the equivalent of chi2=np.sum((data-model)/sigma)**2) done inside the modelling class will have to do another conversion anyway.

(Also, since you convert wavelengths to [m] anyway, you do automatically get a readable output unit, [J / (m3 s)])

I just wanted to mention that I think a blackbody spectrum model would be a great addition.
A few weeks back I wanted to add it myself (see my old version here), but I didn't know how to handle the issue of units in inputs and outputs.

@nden @embray I think this is the first model where using quantities as parameters and possibly as inputs to __call__ could make sense. What do you think?

I think this should be implemented as a class BlackBody that subclasses Parametric1DModel, not a function?

Now I saw that @eteq already made issue #1464 to discuss if / how models should support quantities.

@plim

Sorry I missed this PR ealier. My comments below would be more useful then.

Coincidentally, STScI Synthetic Photometry package (in development) also has similar Planck blackbody implementations. The default units are, of course, suited for Hubble Space Telescope.

https://github.com/spacetelescope/pysynphot/blob/master/synphot/planck.py
https://github.com/spacetelescope/pysynphot/blob/master/synphot/tests/test_planck.py

The magical flux units are defined in

https://github.com/spacetelescope/pysynphot/blob/master/synphot/units.py

particularly

PHOTLAM = u.def_unit(
    'photlam', u.photon / (u.cm**2 * u.s * u.AA),
    format={'generic': 'PHOTLAM', 'console': 'PHOTLAM'})
PHOTNU = u.def_unit(
    'photnu', u.photon / (u.cm**2 * u.s * u.Hz),
    format={'generic': 'PHOTNU', 'console': 'PHOTNU'})
FLAM = u.def_unit(
    'flam', u.erg / (u.cm**2 * u.s * u.AA),
    format={'generic': 'FLAM', 'console': 'FLAM'})
FNU = u.def_unit(
    'fnu', u.erg / (u.cm**2 * u.s * u.Hz),
    format={'generic': 'FNU', 'console': 'FNU'})

Incidentally, I also implemented a BlackBody1D model that inherits from Parametric1DModel that uses the Planck function. I don't think this is much use until modeling can take Quantities due to the lack of unit flexibility.

class BlackBody1D(Parametric1DModel):
    """Create a blackbody spectrum model with given temperature.

    Parameters
    ----------
    temperature : float
        Blackbody temperature in Kelvin.

    """
    temperature = Parameter('temperature')

    def __init__(self, temperature, **constraints):
        super(BlackBody1D, self).__init__(
            temperature=temperature, **constraints)

    def eval(self, x, temperature):
        """Evaluate the model.

        Parameters
        ----------
        x : number or ndarray
            Wavelengths in Angstrom.

        temperature : number
            Temperature in Kelvin.

        Returns
        -------
        y : number or ndarray
            Blackbody radiation in FLAM.

        """
        wave = u.Quantity(np.ascontiguousarray(x), unit=u.AA)
        t = u.Quantity(temperature, unit=u.K)
        bbflux = planck.bb_photlam(wave, t)
        fluxes = bbflux.to(units.FLAM, equivalencies=u.spectral_density(wave))
        return fluxes.value

    def deriv(self, x, temperature):
        raise NotImplementedError(
            'deriv undefined for {0}.'.format(self.__class__.__name__))

Perhaps some of my work would be useful for this PR as well. I would certainly use the Astropy Planck function when it is available.

I'd like to suggest we revisit our previous decision to implement this as a model. Clearly this would be a useful function regardless of astropy.modeling and it seems silly to hold it back just because we don't support units in modeling yet. Maybe we could have a place in astropy where we can collect commonly used astronomy functions, such as blackbodies, king profiles, lorentzian line profiles, etc. We could always consider having model wrappers later, but why not have these already as directly callable functions? What do others think?

@astrofrog, where would such functions reside in astropy - utils or a new sub-package?

Not sure, but doesn't mean it can't be done :) Not utils though, that's for non astro utils (more infrastructure stuff)

The units thing is something I'm working on almost as we speak. But I agree there's no reason it has to be implemented strictly within the Model framework. In fact I've been trying to design things so that basic functions can have any number of implementations outside the Model framework, and can be easily wrapped (possibly even swapping out implementations where one implementation is somehow better than another for some particular data).

astropy.math? astropy.analytical? astropy.analytic_functions?

Scipy uses the very cryptic scipy.special ;)

I don't really care what it's called, but I've said elsewhere (possibly even in this PR) that it would be good to have a module of general analytic functions and such. These would be used by some of the models in the modeling package in their evaluate methods, but not exclusively either. One other example that has come up before is general rotations.

A quick online search reveals that @eteq implemented some models back in 2010 at https://pythonhosted.org/Astropysics/coremods/models.html -- Any plans to port them to astropy?

@pllim - yes definitely. They're mostly functions that can pretty much be imported wholesale into modeling. The last time I thought about this I was waiting for some of the major modeling changes, but now that I think about it, most of them can probably be moved now, as the auto-conversion of functions to models is already available, even though @embray is still finishing up other major modeling changes.

The only uncertainty is about the Quantity question - probably some of them are much more useful with Quantity support in modeling. But some of them it probably doesn't matter so much...

@eteq Not to get too far off topic, but there are no major changes that would significantly impact porting. Quantity support is on its way but again is not a major factor in implementation.

@embray - I just meant that some of those models would be better if they used Quantity because some of the parameters have specific physical meaning. But I agree that's not really a reason to put them in, as they work fine as models even if the units aren't quite right or something.

So I can try to port over all of those that make sense to include (some of them like gaussians and various polynomials are superceded by exostomg astropy.modeling models, but many of them could be moved over.) So I'll try to do that soon (i.e. pre-1.0 feature freeze). Will do it in a separate PR or issue.

@eteq , is it feasible to move the models that need proper unit handling to astropy.analytic_functions first? And then wrap it with modeling in the future?

@pllim - yeah, that might make sense for some of them. Will have to go through the list of models a bit more to decide for sure

Merged in #3094

Thanks @skendrew and @pllim--it will be great to use these in the modeling framework once I've gotten Quantity working with models ^^;

No prob!