时间戳

首页 > 解决方案 > 编写用户定义的函数来评估带有 for 循环的 Saha 方程以获得预期输出

python-3.x - 编写用户定义的函数来评估带有 for 循环的 Saha 方程以获得预期输出

问题描述

我正在尝试创建一个函数来评估某些温度和电子压力值的 Saha 函数。这个问题有点深入，所以我将尽可能详细地介绍本节之前使用的过去代码。

前几节代码

评估分区函数（第 1 部分）：

k= 8.617333262145179e-05  
T=10000. 
g=1.0
Ca_ion_energies = np.array([6.1131554, 11.871719, 50.91316, 67.2732, 84.34])  #in eV
Ca_partition_values= []

def partfunc_E(chiI,T):
    for chiI in Ca_ion_energies:
        elem = 0
        for i in np.arange(chiI):
            elem = elem + (g*np.exp(-(i/(k*T))))
        Ca_partition_values.append(elem)
    return Ca_partition_values

print(partfunc_E(Ca_ion_energies,T))

输出：

[1.455902590894594, 1.45633321917395, 1.4563345239240013, 1.4563345239240013, 1.4563345239240013]

评估玻尔兹曼方程（第 2 部分）：

chiI = np.array([6.1131554, 11.871719, 50.91316, 67.2732, 84.34])  #in eV
k= 8.617333262145179e-05
T=10000.

def boltz_E(chiI,T,I,i):
    Z_1 = partfunc_E(chiI,T)
    ratio = np.exp(-i/(k*T)) / Z_1
    return ratio [I-1] 

print(Ca_ion_energies)
print("i Fraction in level i for I=1 (neutral)")
print("- -------------------------------------")
for n in range(0,10):
    print(n,boltz_E(chiI,10000,1,n))

输出：

[ 6.1131554 11.871719  50.91316   67.2732    84.34     ]
i Fraction in level i for I=1 (neutral)
- -------------------------------------
0 0.6868591389658425
1 0.21522358567610525
2 0.06743914320048579
3 0.021131689732463026
4 0.006621500359539954
5 0.002074811222693332
6 0.0006501308428703751
7 0.0002037149733085943
8 6.383298193775377e-05
9 2.0001718660577703e-05

我需要帮助的问题（以及到目前为止的代码）：

评估 Saha 方程（第 3 部分）：

本节的说明如下：

获得该比率的最简单方法是将 _=1（即中性原子）设置为某个值（例如统一），在 for 循环中从 Saha 方程连续评估下一个电离阶段种群，最后将它们除以由所有在同一尺度上的总和。您会发现 numpy np.sum 函数对于获取所有阶段的总数很有用。我们希望温度 T 为 5000K，电子压力 Pe 为 100.0 N/m^2。

仅供参考：I 是电离阶段，Z_1 是第 1 部分的配分函数，Z_I 是 I+1 阶段的配分函数，Pe 是电子压力，chiI 是电离能量（对于我的代码中的钙），T 是温度而设置“分数”等于的函数就是萨哈方程。

它应该像这样开始：

def saha_E(chiI,T,Pe,I):

计算萨哈人口分数 N_I/N

输入：电离能量、温度、电子压力、离子级

  1. 计算分区函数

  2. 循环遍历您有能量的每个电离阶段，通过 saha 方程计算分数。请注意，第一阶段应设置为 1。

  3. 将每个阶段除以总数

  4. 返回请求阶段的分数

我的代码尝试：

k= 8.617333262145179e-05  
T=10000.  
g=1.0
Ca_ion_energies = np.array([6.1131554, 11.871719, 50.91316, 67.2732, 84.34])  
N_I = 1
h = 6.626e-34
m = 9.11e-31
fractions = []
fraction_sum = []

def saha_E(chiI,T,Pe,I):
    Z_1 = partfunc_E(chiI,T)
    Z_I = partfunc_E(chiI+1,T)
    for I in Ca_ion_energies:
        fraction = (N_I*(Z_I/Z_1)*((2*k*T)/((h**3)*Pe))*((2*np.pi*m*k*T)**(3/2))*np.exp(-I/(k*T)))
        fractions.append(fraction)
        fraction_sum.append(np.sum(fractions))
        for i in fractions:
            i/fraction_sum
            return fraction

print("For ionisation energies (in eV) of:",chiI)
print()
print("I Fraction in stage I")
print("- -------------------")
for I in range(0,6):
    print(I,saha_E(chiI,5000,100.0,I))

我还被指示输出应该类似于：

For ionisation energies (in eV) of: [  6.11  11.87  50.91  67.27  84.34]

I  Fraction in stage I
-  -------------------
1 0.999998720736
2 1.27926351211e-06
3 7.29993420039e-52
4 1.3474665329e-113
5 1.54848994685e-192

首先，我认为我的代码不正确，但这是我能做的最好的，这就是为什么我需要一些帮助，而且，这段代码给了我以下错误：

TypeError: unsupported operand type(s) for /: 'list' and 'list'

如果我的代码完全错误，请告诉我，因为我已经花了很多时间试图弄清楚这一点。

编辑

这个问题还没有完全回答，请继续评论！

标签： python-3.xfor-loopuser-defined-functionsphysicsastronomy

解决方案

如果我很好地理解了您的问题，我的方法是在单个循环中计算各种能量的“分数”和“分数总和”，并且仅在我们在循环之外进行归一化。

另外，请注意代码的范围。我将您在函数外部声明的一些变量推入其中，因为没有理由让它们在函数范围之外保持活动状态。

还要小心不要两次使用相同的变量。您的函数接受一个I参数，但I在 for 循环中有一个变量。

正如聊天中所说，您想编写文档字符串和注释，以便在接触任何代码之前就知道您要去哪里。这是一个要完成的基础：

import numpy as np 

# Constants.
k = 8.617333262145179e-05  
g = 1.0
h = 6.626e-34
m = 9.11e-31
Ca_ion_energies = np.array([6.1131554, 11.871719, 50.91316, 67.2732, 84.34]) # in eV.


# Partition function.
def partfunc_E(chiI, T):
    """This function returns the partition of blablabla.

       args:
       ------
       :chiI: (array or list) the energy levels of a chosen ion.
       :T: (float) the temperature at which kT will be calculated."""

    Ca_partition_values = []
    for energy_level in chiI: # For each energy level.
        elem = 0
        for i in np.arange(energy_level): # From 0 to current energy level.
            elem += g*np.exp(-(i/(k*T)))
        Ca_partition_values.append(elem)
    return np.array(Ca_partition_values) # Conversion to numpy array to support operations later.

print(partfunc_E(Ca_ion_energies, T=10000))


# Boltzmann equation.
def boltz_E(chiI, T, I, i):
    Z_1 = partfunc_E(chiI, T)
    ratio = np.exp(-i/(k*T)) / Z_1
    return ratio[I-1] 

print(Ca_ion_energies)
print("i Fraction in level i for I=1 (neutral)")
print("- -------------------------------------")
for n in range(0,10):
    print(n, boltz_E(Ca_ion_energies, T=10000, I=1, i=n))


# Saha equation.
def saha_E(chiI, T, Pe, i):
    p = partfunc_E(chiI, T)
    Z_ratios = np.array([p[n]/p[0] for n in range(len(chiI))])
    fractions = []
    fractions_sum = []
    for n, I in enumerate(chiI):
        fraction = Z_ratios[n]*((2*k*T)/((h**3)*Pe))*((2*np.pi*m*k*T)**(3/2))*np.exp(-I/(k*T))
        fractions.append(fraction)
        fractions_sum.append(np.sum(fractions))

    # Let's normalize the array before returning it.
    fractions = np.divide(fractions, fractions_sum)
    return fractions[i]


print("For ionisation energies (in eV) of:", Ca_ion_energies)
print()
print("I Fraction in stage n")
print("- -------------------")
for n in range(0, 4):
    print(n, saha_E(Ca_ion_energies, T=5000, Pe=100.0, i=n))

推荐阅读