📚 todo: Specify if vacuum or air in crop function

[minouHub]

💭 Description

In the crop function, if an astropy.unit is fed as a waverange, there is an ambiguity if the wavelength is in air or vacuum. For instance:

from radis import load_spec
from radis.test.utils import getTestFile
import astropy.units as u

# Work with a saved Spectrum object
s2 = load_spec(getTestFile("N2C_specair_380nm.spec"), binary=False)
s2.crop(0.376 * u.um, 0.381 * u.um)

w2 = s2.get_wavelength()
print(w2.min(), w2.max())

In this code, we do not know if the crop function will cut the spectrum in 'nm' or 'nm_air'. I suggest forcing the unit in the crop function and remove this astropy.units feature for this specific function. The use of crop would, in this case, be:

s2.crop(0.376, 0.381, "nm")

@erwanp I know you push for using astropy.units but I think it overcomplicates in this case.

[erwanp]

We'll have the same problem ("are nm in vacuum or in air by default ?") in every context, for instance in calc_spectrum input.

• Astropy also discussed the problem: Air <-> vacuum equivalency? (ANSWER: Use specutils) astropy/astropy#2768 and Functions to convert from air to vacuum wavelengths astropy/astropy#8427 and refers to Specutils
• Specutils has air_to_vac and vac_to_air conversions but as I understand their code, does not expect anything in particular. Up to the user.
• muler users expect wavelength in air Feature Request: Vacuum--> Air OttoStruve/muler#110
• In our own use cases (UV/Visible spectra) we expect wavelengths in air.

Based on this I would Make Air the default, unless stated otherwise in the Spectrum.

• i.e. when given "nm", even from a astropy unit, we read it as "nm_air"
• which means before sending the astropy unit input to the astropy converter, we adjust nm in air to nm in vacuum. There is a bit of machinery involved, but I think it would yield the right behavior for the user.
• we also give the user the option to use "nm_air" as a unit (in crop() for instance) ; to prevent any confusion

Eventually this would be the test :

# The guiding principle is that .crop()  and .plot() use the same waverange , so that it's very clear for the user

from radis import load_spec
from radis.test.utils import getTestFile
import astropy.units as u

#%% Work with a Spectrum object saved in nm_air
s2 = load_spec(getTestFile("N2C_specair_380nm.spec"), binary=False)

# Default crop (astropy units)
# ... we expect it happens in air
s_def = s2.crop(0.376 * u.um, 0.381 * u.um, inplace=False)
w_def = s_def.get_wavelength(medium="air")
# print(w_def.min(), w_def.max())
assert w_def.min() == 376
assert w_def.max() == 381

# Crop in air
s_air = s2.crop(376, 381, "nm_air", inplace=False)
w_air = s_air.get_wavelength(medium="air")
# print(w_def.min(), w_def.max())
assert w_air.min() == 376
assert w_air.max() == 381


# Crop in vacuum
s_vac = s2.crop(376, 381, "nm_vac", inplace=False)
w_vac = s_vac.get_wavelength(medium="vacuum")
# print(w_vac.min(), w_vac.max())
assert w_vac.min() == 376
assert w_vac.max() == 381



#%% Work with a Spectrum object saved in nm_vacuum

# WE DONT HAVE ANY AT THE MOMENT, HAVE WE ? 
# else: adapt the visible spectrum above ; ex s3 = s2.resample(s2.get_wavelength(wunit='nm_vac'), unit='nm_vac') ??



#%% Work with Spectrum object saved in cm-1 

s = load_spec(getTestFile("CO_Tgas1500K_mole_fraction0.01.spec"), binary=False)

# Default crop (astropy units)
# ... we expect it happens in air
s_def = s.crop(4.67 * u.um, 4.69 * u.um, inplace=False)
w_def = s_def.get_wavelength(medium="air")
# print(w_def.min(), w_def.max())
assert w_def.min() == 4670
assert w_def.max() == 4690

# Crop in air
s_air = s.crop(4670, 4690, "nm_air", inplace=False)
w_air = s_air.get_wavelength(medium="air")
# print(w_def.min(), w_def.max())
assert w_air.min() == 4670
assert w_air.max() == 4690


# Crop in vacuum
s_vac = s.crop(4670, 4690, "nm_vac", inplace=False)
w_vac = s_vac.get_wavelength(medium="vacuum")
# print(w_vac.min(), w_vac.max())
assert w_vac.min() == 4670
assert w_vac.max() == 4690

---

EDIT : Decided with @minouHub

• if Astropy.units are used with wunit!='default' ; FAIL. (unless - maybe - astropy.units are also in "nm" )