Bug in CCF?

[mulaura]

Hello,

I think the dopplrShift or the crosscorrRV functionality of pyAsl is bugged by a factor of 2. I used both to generate a Dopplershifted copy and without implementing any Noise the CCF maximizes at a RV of almost always the correct RV multiplied by the factor 2.

import pandas as pd
import pylab as plt
from PyAstronomy import pyasl

# Load the spectrum
headers = ['X', 'Y']
df = pd.read_csv('Proxima_Phoenix_HR_3000_5.0_1.042_2_5.5.dat', sep= " ", names=headers)
X, Y  =  np.asarray(df["X"]* 1e-4), np.asarray(df["Y"] * (df["X"]* 1e-4) * 1e8 ) # in cgs units

plt.figure(figsize=(12,6))
for v in range(-50, +50, 20):
    Y_Doppler, X_Doppler = pyasl.dopplerShift(X,Y,v)
    rv, cc = pyasl.crosscorrRV(X_Doppler, Y_Doppler, X, Y, -100., 100., 10./10., skipedge=1200)
    plt.plot(rv, cc, label="v = " +str(v)+' km/s')

plt.legend(loc='best')
plt.show()

Resulting in:

Does anyone run upon the same issues?

Kind regards, Laura

[sczesla]

Hi Laura, Your issue results from a misunderstanding about the result of dopplerShift, which returns the shifted spectrum on the input wavelength axis and the shifted wavelength axis (on which interpolation was carried out). To get the shifted spectrum, you need to use the original input wvl axis with the first return value like in the example below (and you can basically ignore the second return value of dopplerShift).

Your use essentially applies the shift twice. Yet, I think your assumption about the return value is actually quite natural and might have been avoided with a different return value. Unfortunately, it is too late for an API change here.

I hope this resolves your problem, Stefan

import pandas as pd
import pylab as plt
from PyAstronomy import pyasl
import numpy as np

X = np.arange(4000, 4100, 0.02)
Y = np.ones_like(X) - 0.5 * np.exp(-(X-4050)**2/(2*0.5**2))

plt.figure(figsize=(12,6))
for v in range(-50, +50, 20):
    Y_Doppler, _ = pyasl.dopplerShift(X,Y,v)
    rv, cc = pyasl.crosscorrRV(X, Y_Doppler, X, Y, -200., 200., 10./10., skipedge=200)
    plt.plot(rv, cc, label="v = " +str(v)+' km/s')
    print(f"v = {v} km/s, CCF max at {rv[np.argmax(cc)]} km/s")

plt.legend(loc='best')
plt.show()

[mulaura]

from PyAstronomy import pyasl import numpy as np

X = np.arange(4000, 4100, 0.02) Y = np.ones_like(X) - 0.5 * np.exp(-(X-4050)*2/(20.5**2))

plt.figure(figsize=(12,6)) for v in range(-50, +50, 20): YDoppler, = pyasl.dopplerShift(X,Y,v) rv, cc = pyasl.crosscorrRV(X, Y_Doppler, X, Y, -200., 200., 10./10., skipedge=200) plt.plot(rv, cc, label="v = " +str(v)+' km/s') print(f"v = {v} km/s, CCF max at {rv[np.argmax(cc)]} km/s")

plt.legend(loc='best') plt.show()

Thanks!

Hi Laura, Your issue results from a misunderstanding about the result of dopplerShift, which returns the shifted spectrum on the input wavelength axis and the shifted wavelength axis (on which interpolation was carried out). To get the shifted spectrum, you need to use the original input wvl axis with the first return value like in the example below (and you can basically ignore the second return value of dopplerShift).

Your use essentially applies the shift twice. Yet, I think your assumption about the return value is actually quite natural and might have been avoided with a different return value. Unfortunately, it is too late for an API change here.

I hope this resolves your problem, Stefan

import pandas as pd
import pylab as plt
from PyAstronomy import pyasl
import numpy as np

X = np.arange(4000, 4100, 0.02)
Y = np.ones_like(X) - 0.5 * np.exp(-(X-4050)**2/(2*0.5**2))

plt.figure(figsize=(12,6))
for v in range(-50, +50, 20):
    Y_Doppler, _ = pyasl.dopplerShift(X,Y,v)
    rv, cc = pyasl.crosscorrRV(X, Y_Doppler, X, Y, -200., 200., 10./10., skipedge=200)
    plt.plot(rv, cc, label="v = " +str(v)+' km/s')
    print(f"v = {v} km/s, CCF max at {rv[np.argmax(cc)]} km/s")

plt.legend(loc='best')
plt.show()

Thanks! This works:)