嗨,我试着用一个多项式或指数函数来拟合我的数据,但我都失败了。我使用的代码如下:
with open('argon.dat','r') as f:
argon=f.readlines()
eng1 = np.array([float(argon[argon.index(i)].split('\n')[0].split(' ')[0])*1000 for i in argon])
II01 = np.array([1-math.exp(-float(argon[argon.index(i)].split('\n')[0].split(' ')[1])*(1.784e-3*6.35)) for i in argon])
with open('copper.dat','r') as f:
copper=f.readlines()
eng2 = [float(copper[copper.index(i)].split('\n')[0].split(' ')[0])*1000 for i in copper]
II02 = [math.exp(-float(copper[copper.index(i)].split('\n')[0].split(' ')[1])*(8.128e-2*8.96)) for i in copper]
fig, ax1 = plt.subplots(figsize=(12,10))
ax2 = ax1.twinx()
ax1.set_yscale('log')
ax2.set_yscale('log')
arg = ax2.plot(eng1, II01, 'b--', label='Argon gas absorption at STP (6.35 cm)')
cop = ax1.plot(eng2, II02, 'r', label='Copper wall transp. (0.81 mm)')
plot = arg+cop
labs = [l.get_label() for l in plot]
ax1.legend(plot,labs,loc='lower right', fontsize=14)
ax1.set_ylim(1e-6,1)
ax2.set_ylim(1e-6,1)
ax1.set_xlim(0,160)
ax1.set_ylabel(r'$\displaystyle I/I_0$', fontsize=18)
ax2.set_ylabel(r'$\displaystyle 1-I/I_0$', fontsize=18)
ax1.set_xlabel('Photon Energy [keV]', fontsize=18)
plt.show()
这给了我我想做的是,不是像这样绘制数据,将它们拟合成指数曲线,然后将这些曲线相乘,以获得探测器效率(我试图用元素乘以元素,但我没有足够的数据点来获得平滑的曲线)我尝试使用polyfit,还试图定义一个指数函数来但我在这两种情况下都用了一条线
^{pr2}$以及
model = np.polyfit(eng1, II01 ,5)
y = np.poly1d(model)
#splineYs = np.exp(np.polyval(model,eng1)) # also tried this but didnt work
ax2.plot(eng1,y)
如有需要,数据取自http://www.nist.gov/pml/data/xraycoef/index.cfm 类似的工作也可以在图3中找到:http://scitation.aip.org/content/aapt/journal/ajp/83/8/10.1119/1.4923022
在@Oliver的回答之后,剩下的人都被逗乐了:
我用现有数据做乘法运算:
i = 0
eff1 = []
while i < len(arg):
eff1.append(arg[i]*cop[i])
i += 1
我最后得到的结果是(红色:铜,蓝色虚线:氩,蓝色:乘法)这是我想得到的结果,但是通过使用曲线的函数,这将是一条平滑的曲线,我想用它来结束(在@oliver的回答中有一个关于什么是错误或误解的评论)
因为
curvefit
给你一个常量(一条直线),因为你给它传递的是一个使用你定义的模型不相关的数据集!在让我先重新创建您的设置:
现在请注意,
^{pr2}$f1
基于文件argon.dat
中数据的第2列。它与第一列无关,尽管没有什么可以阻止您绘制第二列与第一列的修改版本,这就是您在绘制时所做的:备注:在您的模型中,有一个名为
b
的参数未使用。这总是一个坏主意传递给一个合适的算法。把它扔掉。在现在有个诀窍:您使用指数模型在第2列的基础上生成了}作为因变量。像这样:
f1
。所以您应该传递curve_fit
第二列作为自变量(在function's doc-string中被标记为xdata
),然后{这将非常有效。在
现在,当你想要绘制两个图的乘积的平滑版本时,你应该从自变量中的一个公共区间开始。对你来说,这就是光子能量。两个数据文件中的第二列取决于:有一个函数(一个用于氩,另一个用于铜),它将
μ/ρ
与光子能量联系起来。因此,如果你有很多能量的数据点,并且你成功地得到了这些函数,那么你将有很多数据点用于μ/ρ
。由于这些函数是未知的,我能做的最好的事情就是简单地插值。但是,数据是对数的,因此需要对数插值,而不是默认的线性插值。在现在,继续获取大量的光子能量数据点。在数据集中,能量点呈指数级增加,因此可以使用^{} 创建一组像样的新点:
两个数据集中的能量具有相同的最小值和最大值,这对我们是有利的。否则,您将不得不缩小这个日志空间的范围。在
接下来,我们(对数)插值关系
energy -> μ/ρ
:以下是对刚刚完成的工作的直观描述:
原始的数据集,用蓝点标记,已经被精细地插值。在
现在,最后的步骤变得容易了。因为已经找到了将
μ/ρ
映射到某个指数变量(我将其重命名为f1
和f2
的函数)的模型参数,它们可以用来生成现有数据的平滑版本,以及这两个函数的乘积:就这样。总结一下:
photon energy -> μ/ρ
μ/ρ
相关问题 更多 >
编程相关推荐