1

我正在尝试在 python 中重现在 IDL 中正常工作的东西。（这是一个更大程序的一部分，所以不能只使用 IDL。）

我正在尝试做的是将尚未经过波长校准的光谱拟合到另一个已经具有的光谱，从而获得未校准光谱的波长色散解决方案。

IDL 代码：

function fit_wav, x, p, wav=w, flux=f ;our function for mpfitfun
    ;eval a polynomial function for wavelength and interpolate
    ;to the wavelength of the fiducial spectrum
    wav = poly(w, p)
    return, interpol(f, wav, x)
end

function normalize, arr ;scale an array into [0,1]
    return, (arr - min(arr)) / (max(arr) - min(arr))
end

readcol, '~/Desktop/fiducial.txt', wavelength, flux, format='F,F'
readcol, '~/Desktop/to_fit.txt', uncalibrated, format='F'

npix = n_elements(uncalibrated)
pix = [0d:npix-1]
p_init = [1.95, 0.5/npix, 0.]
coef = mpfitfun('fit_wav', wavelength, normalize(flux), 1., $
                p_init, functargs={wav:pix, flux:normalize(uncalibrated)})
print, coef

;plot the results
calib_wav = poly(pix, coef)
cgplot, calib_wav, normalize(uncalibrated), xra=minmax([calib_wav, wavelength])
cgplot, /overplot, wavelength, normalize(flux), color='red'
end

这个输出看起来像（基准是红色的，适合是黑色的）：

    1.9652186   0.00070232586  -1.5183943e-07

适合 IDL 中的 MPFITFUN

python代码基本相同：

from astropy.modeling import Fittable1DModel, Parameter, fitting
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

class WavModel(Fittable1DModel):
    '''Our model for fitting the wavelength dispersion solution
    based on an existing calibrated spectrum'''
    a = Parameter(name='a',default=0.)
    b = Parameter(name='b',default=1.)
    c = Parameter(name='c',default=0.)

    def __init__(self, a, b, c, spec=None, **kwargs):
        self.spec = spec
        self.pix = np.arange(spec.size, dtype=float)
        super(WavModel, self).__init__(a, b, c, **kwargs)

    def evaluate(self, x, a, b, c):
        wav = a + self.pix * b + self.pix * self.pix * c
        return interp1d(wav, self.spec, bounds_error=False, fill_value=0.)(x)

    fit_deriv = None

normalize = lambda arr: (arr - arr.min()) / (arr.max() - arr.min()) #scale arr into [0,1]

wavelength, flux = np.loadtxt('fiducial.txt').T
uncalibrated = np.loadtxt('to_fit.txt')

npix = uncalibrated.size
pix = np.arange(npix, dtype=float)
fitter = fitting.LevMarLSQFitter()
p_init = WavModel(1.95, 0.5/npix, 0., spec=normalize(uncalibrated))
coef = fitter(p_init, wavelength, normalize(flux)).parameters

print coef

#now we plot!
calib_wav = coef[0] + coef[1] * pix + coef[2] * pix * pix
plt.plot(calib_wav, normalize(uncalibrated), 'b')
plt.plot(wavelength, normalize(flux), 'r')
plt.savefig('python_test.png')

然而，THIS 的结果有些不同：

[  1.98065438e+00   3.70817592e-04   1.06795564e-06]

Python适合分散解决方案。

IDL 版本绝对是正确的版本，因为检测器的波长限制在 2.5 微米左右……

我在这里和这里粘贴了数据。关于如何让 python 看起来像 IDL 合适的任何建议？太感谢了！

4

0 回答 0