Changing the line for nu=1:nu_max to for nu=0:(nu_max-1) makes julia agree with python.
On Wednesday, July 16, 2014 1:55:40 PM UTC-7, Andrius Popovas wrote: 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 >
