Hello,
I have a code for calculating some physical quantities. I had a very rough
day trying to understand, why my results differ so much from what I should
get. Up until I have tried the same code in Python and I got the correct
result. Both codes are completely identical, but why is there a difference?
Julia code:
# THIS IS FOR HETERONUCLEAR PFs
Be = [1.8997,1.7153,1.973]
De = [0.0,0.0,0.0]
He = [0.0,0.0,0.0]
we = [2163.9,1812.33,2068.61]
alfa_e = [0.0175,0.01727,0.023]
beta_e = [0.0,0.0,0.0]
wexe = [13.114,12.61,20.2]
gamma_e=[0.0,0.0,0.0]
weye = [0.0,0.0,0.0]
T_e = [0,8215.01,24725.99]
w_e=[2,4,2]
mu = 6.860525081309955 #m1*m2/(m1+m2)
sigma = 1 #1 for hetero, 2 for homonuclear molecules
g_n = [1,2,2,2,2,2]
T = 2000 #temperature
HPLNCK=6.62554e-27
BOLTZK=1.38054e-16
CLIGHT=2.997925e10
const c2=(HPLNCK*CLIGHT)/BOLTZK
QQ=0
#over all electronic states in consideration
for el=1:length(T_e)
nu_max=trunc(we[el]/(2*wexe[el])) #Determine nu_max
nu_max2= trunc(Be[el]/alfa_e[el]-0.5) #check nu_max
if nu_max>nu_max2
nu_max=nu_max2
end
#over all vibrational states until nu_max
for nu=1:nu_max
aa=(w_e[el]*T)/(sigma*c2) #first term in summation, degeneracies
bb=-c2/T*(we[el]*(nu+0.5)-wexe[el]*(nu+0.5)^2+T_e[el]) #exponential part
cc=(Be[el]-alfa_e[el]*(nu+0.5))
QQ=QQ+aa*exp(bb)/cc
end
end
first=exp(c2/T*(we[1]/2-wexe[1]/4))
Q_AB=first*QQ
println(T," ",Q_AB)
Python code:
import numpy
Be = [1.8997,1.7153,1.973]
De = [0.0,0.0,0.0]
He = [0.0,0.0,0.0]
we = [2163.9,1812.33,2068.61]
alfa_e = [0.0175,0.01727,0.023]
beta_e = [0.0,0.0,0.0]
wexe = [13.114,12.61,20.2]
gamma_e=[0.0,0.0,0.0]
weye = [0.0,0.0,0.0]
T_e = [0,8215.01,24725.99]
w_e=[2,4,2]
mu = 6.860525081309955 #m1*m2/(m1+m2)
sigma = 1 #1 for hetero, 2 for homonuclear molecules
g_n = [1,2,2,2,2,2]
T = 2000 #temperature
#constants
HPLNCK=6.62554e-27
BOLTZK=1.38054e-16
CLIGHT=2.997925e10
c2=(HPLNCK*CLIGHT)/BOLTZK
QQ=0
#over all electronic states in consideration
for el in range(len(T_e)):
nu_max=int(we[el]/(2*wexe[el])) #Determine nu_max
nu_max2= int(Be[el]/alfa_e[el]-0.5) #check nu_max
if nu_max>nu_max2:
nu_max=nu_max2
#over all vibrational states until nu_max
for nu in range(nu_max):
aa=(w_e[el]*T)/(sigma*c2) #first term in summation, degeneracies
bb=-c2/T*(we[el]*(nu+0.5)-wexe[el]*(nu+0.5)**2+T_e[el])
#exponential part
cc=(Be[el]-alfa_e[el]*(nu+0.5))
QQ=QQ+aa*numpy.exp(bb)/cc
first=numpy.exp(c2/T*(we[0]/2-wexe[0]/4))
Q_AB=first*QQ
print "T= ",T,"Q= ",Q_AB
