Написание пользовательской функции для оценки уравнения Саха с помощью циклов for для ожидаемого результата

Question

Я пытаюсь создать функцию, которая оценивает функцию Саха для определенных значений температуры и электронного давления. Вопрос немного углублен, поэтому я предоставлю как можно больше подробностей о прошлом коде, использованном перед этим разделом.

Код предыдущих секций

Оценка функции разделения (часть 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

Вопрос, с которым мне нужна помощь ( и мой код пока):

Оценка уравнения Саха (часть 3):

Инструкции для этого раздела следующие:

Самый простой способ получить это соотношение - установить ?_? = 1 (т.е. нейтральный атом) на некоторое значение (например, единицу), последовательно оценить популяции следующей стадии ионизации из уравнения Саха в a для l oop, и в конце разделите их на сумму всех ? в одной шкале. Вы найдете функцию numpy np.sum полезной для получения суммы по всем этапам. Мы хотим, чтобы температура T составляла 5000 К, а электронное давление Pe составляло 100,0 Н / м ^ 2.

К вашему сведению: I - ступень ионизации, Z_1 - функция разделения из части 1, Z_I - функция разделения для ступени I + 1, Pe - электронное давление, chiI - энергии ионизации (для Кальция в моем коде), T - температура, а функция, для которой «фракция» установлена ​​равной, - это уравнение Саха.

начать что-то вроде:

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

вычислить долю населения саха N_I / N

вход: энергии ионизации, температура, электронное давление , ионный каскад

1. Вычислить функции разделения

2. L oop для каждой ступени ионизации, для которой у вас есть энергия, вычисляя фракцию с помощью уравнение Саха. Обратите внимание, что для первой ступени должно быть установлено значение 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'

Если мой код абсолютно неверен, пожалуйста, скажите мне, поскольку я потратил так много времени, пытаясь выяснить это.

Редактировать

Этот вопрос до сих пор не полностью отвечен, пожалуйста комментируйте!

Answer 1

Если я хорошо понял вашу проблему, мой подход состоит в том, чтобы вычислить «дроби» и «суммы дробей» в одном l oop для различных энергий и нормализовать только тогда, когда мы находимся за пределами l oop.

Кроме того, будьте осторожны с объемом вашего кода. Я поместил некоторые переменные, которые вы объявили, вне функции внутри нее, потому что нет причин поддерживать их вне области действия функции.

Будьте осторожны и не используйте одну и ту же переменную дважды. Ваша функция принимает аргумент I, но затем имеет переменную I в a для l oop.

Как сказано в чате, вы хотите написать строки и комментарии, чтобы вы знали, куда вы идете, даже прежде чем касаться какого-либо кода. Вот база для завершения:

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))