search

astropy / specutils

An Astropy coordinated package for astronomical spectroscopy. Maintainers: @rosteen @keflavich @eteq
http://specutils.readthedocs.io/en/latest/
166 stars 126 forks source link

Is the default implementation in`specutils.utils.wcs_utils.air_to_vac` incorrect? #1162

Open andycasey opened 3 weeks ago

andycasey commented 3 weeks ago

Hi folks,

I make a lot of model/data comparisons where model spectra are in vacuum and data are in wavelengths.

Today I noticed that the default method in specutils (scheme='inversion', method='Griesen2006') for converting air wavelengths to vacuum seems to be either incorrect, or inaccurate.

I attach a figure showing a ESO/HARPS spectrum of alpha Cen A, where the ESO/HARPS spectrum is recorded in air wavelengths. In the top panel I compare the data (in air wavelengths) with some approximate spectral model (in vacuum wavelengths) in pink. The subsequent panels also show the model in pink, but the data are shifted to vacuum wavelengths using different arguments in specutils.utils.wcs_utils.air_to_vac. I have added a small y-offset so we can see the conversion using different methods.

specutils-air-to-vac-issue

You can see that in either scheme (inversion/iteration), the Greisen 2006 conversion is vastly different from all other methods. This seems to be inaccurate or incorrect, and the unfortunately the Greisen 2006 method is the default in specutils.utils.wcs_utils.air_to_vac. The docstring for air_to_vac says that 'Griesen et al. (2006) report that the error in naively inverting is less than 10^-9.'. The difference in the figure between Greisen 2006 and the other methods is about 0.07 Angstroms.. far higher than 10^-9!

I haven't cross-checked the equations in the Greisen 2006 paper to see whether the implementation in specutils is correct or not. I wanted to open this issue for discussion first, particularly since users will get the wrong result by using the default parameters, and debugging air/vacuum issues with noisy RV/BERV shifts can be fiddly!

andycasey commented 3 weeks ago

A note from the International Association of Geodesy rendus (page 120), which is where the formula that Greisen et al. 2006 uses, is that the 0.25 ppm accuracy is quoted for wavelengths between 650 nm and 850 nm.

I'd be a little surprised if it was so inaccurate outside of that range and it gets adopted as the standard, but I can't see any obvious error in the specutils implementation of the equation.

I tried computing the vacuum wavelengths myself using the calculation for group refractivity (instead of phase refractivity) thinking that the equations might be mixed up, and there is a difference of about 0.07 Angstroms in my example, but it goes in the wrong way (e.g., it makes the 'Greisen' wavelengths 0.14 Angstroms away from all other methods).

Unless I have made some silly mistake here, I think the default scheme or method for air_to_vac should be changed.

rosteen commented 3 weeks ago

Thanks for the detailed report - I'm hoping I'll have some time to look into this/think about it next week.

weaverba137 commented 3 weeks ago

I did some work on this a long time ago this notebook has some code and plots.

andycasey commented 3 weeks ago

Thanks for the write-up @weaverba137! I'm glad to know that someone else has encountered the same behaviour.

The equations in the Greisen 2006 paper actually go back to the 90's: Ciddor and Hill (1999), and Ciddor et al. (1996).

I found a few papers referring to errors in those papers (e.g., https://opg.optica.org/ao/fulltext.cfm?uri=ao-59-31-9771&id=441874), but all of those papers are published in Optica, which for some reason, my poor university does not have a subscription with.