自昨天以来,我一直在寻找答案,但是没有运气。因此,我有一个1D光谱(.fits)文件,其中每个波长处的通量值。我已将它们转换为2D数组(x,y)=(波长,通量),并希望编写一个程序,该程序将在某些指定的波长(x)下返回通量(y)。我已经试过了:

#modules
import scipy
import numpy as np
import pyfits as pf

#Target Global Vaiables
hdulist_tg = pf.open('cutmask1-2.0001.fits')
hdr_tg = hdulist_tg[0].header
flux_tg = hdulist_tg[0].data
crval_tg = hdr_tg['CRVAL1']             #Starting wavelength
cdel_tg = hdr_tg['CDELT1']              #Wavelength axis width
wave_tg = crval_tg + np.arange(3183)*cdel_tg        #Create an x-axis
wavelist = [6207,6315,6369,6438,6490,6565,6588]

wave_flux=[]
diff = 10
for wave in wave_tg:
    for flux in flux_tg:
        wave_flux.append((wave,flux))

for item in wave_flux:
    wave = item[0]
    flux = item[1]

#Where I got my actual wavelength that exists in wave_tg
    diffmatch = np.abs(wave - wavelist[0])
    if diffmatch < diff:
        flux_wave = flux
        diff = diffmatch
        wavematch = wave

print wavelist[0],flux_wave,wavematch


但是即使波长不同,程序也总是返回相同的通量值。请帮忙...

最佳答案

我将完全跳过二维表的创建,而只使用interp

fluxvalues = np.interp(wavelist, wave_tg, flux_tg)


对于您发布的文件,由于wave_tg数组的硬编码长度,您发布的代码不起作用。因此,我建议您宁愿使用

wave_tg = crval_tg + np.arange(len(flux_tg))*cdel_tg


另外,由于某种原因,您发布的文件似乎并没有达到您要查找的波长。您可能需要检查您是否正确计算了相应的波长,或者检查了正确的波长。

10-07 18:59
查看更多