I have this code which calculates the energy produced in fusion processes (PP-chain and CNO cycle), as well as the relative energy, which i want to plot against a temperature interval. But when i run the program i get the error: 'The truth value of an array with more than one element is ambiguous'.
This is a problem with the if-statements in the code, i figured, as T is an array not a scalar. How do I solve this, i.e. implement the if conditions such that the code works and i get a ,correct, plot.
def pp_chain(rho,T):
N_A = 6.022e23 #Avogadro's constant inverse
mu = 1.6605e-27 #Units of [kg]
Joule = 1.60217e-13 #Converting Mev to Joule (SI)
#Mass fraction of each nuclei
X= 0.7
Y = 0.29
Y3He = 1e-10
Z7Li = 1e-7
Z7Be = 1e-7
Z14N = 1e-11
#Q-values of the different reactions, i.e. energy released
Q_pp = 1.177*Joule
Q_pd = 5.494*Joule
Q_33 = 12.86*Joule
Q_34 = 1.586*Joule
Q_e7 = 0.049*Joule
Q_17s = 17.346*Joule
Q_17 = 0.137*Joule
Q_8 = (8.367+2.995)*Joule
Q_CNO = (1.944 + 1.513+7.551 + 7.297+1.757 + 4.966)*Joule
#number density of each nuclei (in each reaction)
n_p = ((rho*X)/(mu))
nHe3 = (rho*Y3He)/(3*mu)
nHe4 = ((rho*Y)/(4*mu))
ne7 = (rho*Z7Be/(7*mu))
n17s = ((rho*Z7Li)/(7*mu))
nCNO = ((rho*Z14N)/(14*mu))
n_e = n_p+ 2*nHe4 + 2*nHe3 + 4*ne7 + 3*n17s + 7*nCNO
T9 = T/(1e9)
#The reaction rates (per unit mass)for each reaction
#in reaction lambda_e7.
lambda_pp = ((4.01e-15*T9**(-2/3)*np.exp(-3.380*T9**(-1/3))*(1+0.123*T9**(1/3)+1.09*T9**(2/3)+0.938*T9))/N_A)*0.000001
lambda_33 = ((6.04e10*T9**(-2/3)*np.exp(-12.276*T9**(-1/3))*(1+0.034*T9**(1/3)-0.522*T9**(2/3)-0.124*T9+0.353*T9**(4/3)+0.213*T9**(5/3)))/N_A)*0.000001
Ts = T9/(1+4.95e-2*T9)
lambda_34 = ((5.61e6*Ts**(5/6)*T9**(-3/2)*np.exp(-12.826*Ts**(-1/3)))/N_A)*0.000001
if T<1e6:
return lambda_e7 =((1.57e-7/n_e)/N_A)*0.000001
else:
return lambda_e7 = ((1.34e-10*T9**(-1/2)*(1-0.537*T9**(1/3)+3.86*T9**(2/3)+0.0027*T9**(-1)*np.exp(2.515e-3*T9**(-1))))/N_A)*0.000001
Tss = T9/(1+0.759*T9)
lambda_71Li = ((1.096e9*T9**(-2/3)*np.exp(-8.472*T9**(-1/3))-4.830e8*Tss*(5/6)*T9**(-3/2)*np.exp(-8.472*Tss**(-1/3))+1.06e10*T9**(-3/2)*np.exp(-30.442*T9**(-1)))/N_A)*0.000001
lambda_71Be = ((3.11e5*T9**(-2/3)*np.exp(-10.262*T9**(-1/3))+2.53e3*T9**(-3/2)*np.exp(-7.306*T9**(-1)))/N_A)*0.000001
lambda_CNO = ((4.90e7*T9**(-2/3)*np.exp(-15.228*T9**(-1/3)-0.092*T9**2)*(1+0.027*T9**(1/3)-0.778*T9**(2/3)-0.149*T9+0.261*T9**(4/3)+0.127*T9**(5/3)) + 2.37e3*T9**(-3/2)*np.exp(-3.011*T9**(-1))+2.19e4*np.exp(-12.53*T9**(-1)))/N_A)*0.000001
##The reaction rates for each reaction
r_pp = (n_p*n_p)/(2*rho)*lambda_pp
r_33 = (nHe3*nHe3)/(2*rho)*lambda_33
r_34 = ((nHe4*nHe3)/rho)*lambda_34
r_e7 = ((ne7*n_e)/rho)*lambda_e7
r_71Li = ((n17s*n_p)/rho)*lambda_71Li
r_71Be = ((ne7*n_p)/rho)*lambda_71Be
r_CNO = ((nCNO*n_p)/rho)*lambda_CNO
if r_pp<(r_34+2*r_33):
R1 = r_pp/(2*r_33+r_34)
r_33 = r_33*R1
r_34 = r_34*R1
if (r_34<(r_e7 + r_71Be)):
R2 = r_34/(r_e7 + r_71Be)
r_e7 = r_e7*R2
r_71Be = r_71Be*R2
if (r_e7 <(r_71Li)):
R3 = r_e7/(r_71Li)
r_71Li = R3*r_71Li
#Energy generation per unit mass from PP1,PP2, PP3 and CNO, found by multiplying energy released
#with reaction rates of the reactions that occur in each branch
PP0 = r_pp*(Q_pp + Q_pd)
PP1 = r_33*Q33+ r_33*2*(Qpp+Qpd)
PP2 = r_34*Q_34 + r_e7*Q_e7 + r_71Li*Q_17s + r_34*(Qpp+Qpd)
PP3 = r_34*Q_34 + r_71Be*Q_17 + r_34*(Qpp + Qpd)
CNO = r_CNO*Q_CNO
epsilon = PP1+PP2+PP3+CNO
#Relative energy
E_rel1= PP1/epsilon
E_rel2= PP2/epsilon
E_rel3 = PP3/epsilon
E_rel4 = CNO/epsilon
return E_rel1, E_rel2,E_rel3,E_rel4
#Defining the different temperatures used in function
rho = 1.62e5
T = np.arange(1e4,1e9,100)
#Calling the function for energy generated in eanch branch and printing the result
energy_prod = pp_chain(rho, T)
rel_energy1 = energy_prod[0]
rel_energy2 = energy_prod[1]
rel_energy3 = energy_prod[2]
rel_energy4 = energy_prod[3]
plt.plot(T, rel_energy1,label = 'PP1')
plt.plot(T, rel_energy2,label = 'PP2')
plt.plot(T, rel_energy3,label = 'PP3')
plt.plot(T, rel_energy4,label = 'CNO')
plt.title('Energy prooduced in each PP chain/CNO cycle')
plt.ylabel('Relative Energy [J]')
plt.xlabel('T [K]')
plt.xscale('log')
#plt.yscale('log')
plt.legend()
plt.show()
I did attempt to use np.all/np.any which it said i should use in the error. Doing this i did get a plot, except two of the relative energies (CNO and PP3?) did not behave as expected, thus i concluded that i shouldn't use that. After this i tried with np.where, though i was unsure how to implement it with the if statements for r_pp and r_34 as the conditions return two values if condition is true.
