从昨天开始我一直在寻找答案,但没有运气。所以我有一个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
但是程序总是返回相同的通量值,即使波长不同。请帮忙。。。在
我对你的代码做了一些修改:
使用numpy ot创建}
wave_flux
作为ndarray
,使用np.hstack()
、np.repeat()
和{使用花式索引获取与搜索匹配的值
产生的代码是:
它将返回
^{pr2}$wave_flux
的子组,其波值在0
列中,通量值在1
列中:我将完全跳过二维表的创建,只使用^{} :
对于您发布的文件,由于wave_tg数组的硬编码长度,您发布的代码不起作用。因此,我建议您使用
^{pr2}$你似乎也在找一些波长不高的文件。您可能需要检查您是否正确地计算了相应的波长,或者检查是否查找了正确的波长。在
相关问题 更多 >
编程相关推荐