sczesla / PyAstronomy

A collection of astronomy-related routines in Python
154 stars 35 forks source link

Interpolation to wrong wavelength vector in pyasl.dopplerShift #5

Closed DanielAndreasen closed 9 years ago

DanielAndreasen commented 9 years ago

Hi Stefan,

I have looking at your very simple code in pyasl.dopplerShift for more than an hour, and found an error. When you use the interpolation in line 72 (at my version), you use the old wavelength vector. Instead you should use the shifted wavelength vector. Without the right wavelength, it shift wrong with a factor of 2.

So, change

nflux = f(wvl)

to

nflux = f(wlprime)

Cheers, Daniel

sczesla commented 9 years ago

Hi Daniel,

thanks for your scrutiny. I believe the code is correct, when used as intended, but also unclear (or counter-intuitive).

Line 71:

f = sci.interp1d(wlprime, flux, bounds_error=False, fill_value=fv)

essentially gives the shifted spectrum and line 72:

nflux = f(wvl)

the shifted spectrum on the original wavelength axis.

If you start with wvl and flux as your starting values and obtain

nflux, wlprime = dopplerShift(...)

from the doppler shift, then the shifted spectrum is given by (wvl, nflux) and/or (wlprime, flux). If you use (wlprime, nflux), you get a double shift, indeed.

When writing this, I thought it was a good idea, because I mostly need the shifted spectrum at the input wavelengths, but sometime one would need the updated wavelength axis. However, realizing this "cross-over" between input and output, it may not have been such a good idea.

Cheers, Stefan

DanielAndreasen commented 9 years ago

I see, thanks for explaining. For me it was not obvious at least, but it makes more sense now. Personally, I prefer the "solution" I gave, but I guess that is a matter of taste.

Cheers, Daniel

sczesla commented 9 years ago

Well, I think your approach is more obvious. Unfortunately, it is hard to change now, because it would break existing code. I will try to make it more explicit, nonetheless.

Cheers, Stefan