Large velocity values cause strange error in pyasl.dopplerShift

sczesla / PyAstronomy

A collection of astronomy-related routines in Python

152 stars 35 forks source link

Large velocity values cause strange error in pyasl.dopplerShift #41

Closed texadactyl closed 6 years ago

texadactyl commented 6 years ago

E.g.

VELOCITY = 1e6
nflux1, wlprime1 = pyasl.dopplerShift(wvl, flux, VELOCITY, edgeHandling="firstlast")

File "/usr/local/lib/python2.7/dist-packages/PyAstronomy/pyasl/asl/dopplerShift.py", line 96, in dopplerShift nflux[firsts] = firstval UnboundLocalError: local variable 'firstval' referenced before assignment

If there are velocity boundary value restrictions, they should be tested for. Agree?

texadactyl commented 6 years ago

I noticed that this issue is only for the edgeHandling="firstlast" case. Leaving off the edgeHandling= parameter or using edgeHandling="fillValue" with fillValue=6000. did not cause this exception with a velocity of 1e6 or even 3e8 although their results were unsatisfactory.

edgeHandling="fillValue", fillValue=6000.: Solid straight-line red on top and solid straight-line green on bottom (parallel). No appearance of blue.
No edgeHandling= specification: initial (blue) curve with a solid straight-line green on bottom. No appearance of red.

texadactyl commented 6 years ago

In line 64 of PyAstronomy/src/pyasl/asl/dopplerShift.py, there is a hard-coded speed of light value. I would suggest using the constant [SpeedOfLight].

I noticed that the 2 for-loops (beginning at line 80 and 88 respectively) fail to set a value for firstval and lastval respectively when there are no NAN values. Perhaps, initialize them to illegal index values and then check before using. If the check fails, raise with more specific information like at line 69.

sczesla commented 6 years ago

Hi,

You are quite right with what you say. The idea behind this particular implementation was to apply the Doppler shift and maintain the original wavelength axis. A very frequent use case at least for me. As the shifted spectrum necessarily does not cover the entire original range, the edgeHandling option was introduced.

Large shifts, such as those you try to apply, (can) entirely remove the overlap between input and output spectrum, which leads to nonsense results. Admittedly, this case is not currently recognized, which should be changed. May I ask you for your use case? Such large shifts may need some more scrutiny in terms of the transformation of the flux density also...

As far as the speed of light is concerned, yes, a constant would be somewhat nicer (preventing typos etc.). However, as this is a definition, which is not likely to change, I opted for hard coding it to have it self-contained.

texadactyl commented 6 years ago

I am suggesting that you simply put in protections for out-of-bounds conditions that I noted earlier. A big effort should not be required.

Maybe I am picky but I am a software engineer since the 1960s. My interest in Astronomy is related to my experiences as a student. I do not as yet have a use case; terms like "doppler shift" catch my eye so I investigated.

jason-neal commented 6 years ago

One other thing to note is that speeds approaching 3e8 the realitivistic formuation is required. https://en.wikipedia.org/wiki/Relativistic_Doppler_effect. (I am not suggesting that the dopplerShift function be turned relativistic as it is fine for most uses)

sczesla commented 6 years ago

@texadactyl : Thanks for pointing it out. I agree it is a good idea to (try to) prevent misuses by oneself and others. I implemented a check.

@jason-neal : Totally agree. I am now checking for 5% of c as a default. I also think we should keep it easy and provide a relativistic treatment elsewhere if needed or somebody is in the mood for coding it.