Use astropy.units?

[karllark]

We could use astropy.units to specify the units of the spectra. Then we could ensure they are in the units we expect. It may even be transparent through the use of the astropy.models as long as we setup things correctly. Worth thinking about.

[jdtsmith]

Units are a bit of a pain. Wavelength is straightforward. But right now you have two choices for flux input in original PAHFIT: MJy/sr (in which case all integrated strengths/uncertanties/etc. are in W/m^2/sr), or "some other f_nu unit", in which case integrated values are [Unit]-Hz. Would we want to support non-f_nu units? Any reason to? I don't know if astropy.units can handle composite unit systems like f_lambda vs. f_nu.

[karllark]

Yes, astropy.units can handle a wide range of units. Basically, as long as the user specifies the wavelength and spectral units, we can use astropy.units to convert our preferred units. And convert back to the user specified units for output. And we could also have default units if they are not specified (probably with a warning stating what defaults are being assumed).

[jdtsmith]

Interesting. I suspect we'd prefer everything in f_nu units internally, so we can integrate f_nu dnu wherever it's needed. Even with user-specified units, we'd still have to make choices for output units, since many outputs are (weighted) integrals over flux density. E.g. if that give erg/s/cm^2/angstrom, we still need to provide some composite flux unit output (e.g. erg/s/cm^2).

[jdtsmith]

One other niggling issue with the LM optimization PAHFIT currently does is it doesn't like very small numbers. So depending on your units choice, the output can differ. I often recommend people to normalize the spectrum to do the fit in numbers near unity, then rescale. Hoping we can avoid that, or perhaps it won't be necessary. But a good TEST would be to give PAHFIT the same exact spectrum in a variety of units, and see if it produces the same output (after unit conversion).

[jdtsmith]

We discussed this issue on the call today. One interesting question is what JWST pipeline output spectra will look like, in terms of their units attachment. Another issue is that specifying input spectral units (e.g. MJy/sr) is not sufficient to specify an output integrated power unit (e.g. W/m^2/sr). We could either make that choice for the user, or have a default they could (rarely) modify.

[jdtsmith]

Should we re-open this?

[karllark]

Probably should reopen. Have no memory of why I closed it.

[karllark]

For JWST, units of cubes will be MJy/sr. Units of extracted spectra will be in Jy. Both products will be outputs from the pipeline. Of course, offline work will likely improve both products - especially the extracted spectra.

[jdtsmith]

I think we should abstract out all of our units handling — input -> modeling space -> output — outside of the main model itself. So the main model doesn't even think about units at all, other than to understand them to be "some f_nu variant".

Such a smart pre/post-processor, ala pahfit/units.py, could do a few things;

• use astropy.units to check that the passed flux unit conforms with "f_nu" or "f_nu/sr" flavor.
• Comes up with a scaled internal modeling unit that places the flux scale at numbers of order unity (from experience with LM fitters, this is far preferable to scales in which numbers are 1e-20 or similar).
• Takes the resulting feature strengths, etc., and converts them back to "output" units which are based on the inputs.

[jdtsmith]

Reviving this old issue; @drvdputt discovered new_flux_unit. This works:

flux=(np.random.random(100)*5+5)*1e-17 * u.erg/u.s/u.cm**2/u.angstrom
sp=Spectrum1D(spectral_axis=np.linspace(5,25,100)*u.um, flux=flux)
print(np.mean(sp.flux))
sp = sp.new_flux_unit(u.mJy)
print(np.mean(sp.flux))

7.162028594446536e-17 erg / (Angstrom cm2 s)
62.047523994297244 mJy

So I say for flux, we accept either flux or surface brightness (in any units), then turn them into our "working units" of mJy or MJy/sr.

We can tell the difference between flux and SB like this:

uparts = sp.flux.unit.decompose()
isSB = (u.rad, -2) in zip(uparts.bases, uparts.powers)

For spectral_axis, anything (GHz, nm, etc.) that can be converted into micron is fine (our working unit).

For equivalent power units, for flux density, use W/m^2, and for SB, W/m^2/sr.

All outputs are in our working wavelength, flux-density/surface-brightness, power units (and we will tag the table columns accordingly).

One issue is that while MJy/sr and mJy are nice lowish numbers in most places, not too different than one, this is not so for power. So we could opt for power units of 1e-10 W/m^2/sr and 1e-22 W/m^2 (the latter equivalent to the former uniformly over a solid angle 0.21" on a side: about a small JWST IFU pixel). A bit more complicated, but much easier to read in tables, as these will be numbers like 500 or 20.

These custom units can easily be created as (e.g.):

sbp_unit = CompositeUnit(1e-10,(u.W, u.m, u.sr), (1, -2, -1))

[jdtsmith]

Proposed PAHFIT working + output units summary:

Type	Working Unit
Wavelength/FWHM/spectral_axis	µm
Temperature	K
Flux density	mJy
Flux power	1e-22 W/m2
— OR —
Surface Brightness/Intensity	MJy/sr
Surface Brightness Power	1e-10 W/m2/sr

[jdtsmith]

I could be talked in to µJy for flux density though :).

[jdtsmith]

I have prepared a units.py file to encapsulate this, in the features PR #241. @drvdputt have a look. I'd like to get dev merged soon, so we can add stuff like this.