title: Suggestions on matlab/python conversion
body: <p>I'm trying to translate the matlab code below to a python code. The code calculates numerical the para-state of a deuterium molecule and then plots the result. When I try to translate it to python, it seems that I get stuck in a nested for-loop which calculates a sum. I have been searching on the internet the past days yet without success.</p>
  
  <p>Because it's a physics code I will mention some aspect from the code. So first we calculate the partition function (Z). After that there is a calculation of the energy which is a partial derivative of ln(Z) to beta. From this we can calculate the specific heat (approximately) as the derivative of energy to temperature. </p>
  
  <p>So the matlab code looks like this:</p>
  
  <pre><code>epsilon = 0.0038*1.60217662*10^-19;
  k = 1.38*10^-23;
  T = 1:.1:2000;
  beta = 1./(k*T);
  
  %partitionfunction
  clear Z Zodd;
  for i = 1:length(T)  
      clear p;    
      for s = 1:2:31;
          a = 2*s+1;
          b = s^2+s;
          p(s) = 3*a*exp(-b*epsilon*beta(i));
      end
      Zodd(i) = sum(p);
  end
  
  %energy
  ln_Zodd = log(Zodd); 
  
  for i = 1 : (length(T)-1)
      Epara(i) = -(ln_Zodd(i+1)-ln_Zodd(i))/(beta(i+1)-beta(i));
  end
  
  %heat capacity
  for i = 1 : (length(T)-2)
      Cpara(i) = (Epara(i+1)-Epara(i))/(T(i+1)-T(i));
  end   
  
  %plot
  x = k*T/epsilon;
  plot(x(1:6000),Cpara(1:6000)/k, 'r');
  axis([0 7 0 1.5]);
  ylabel('C_v/k');
  xlabel('kT/eps');
  </code></pre>
  
  <p>The corresponding python code:</p>
  
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  import math
  
  epsilon=0.0038*1.60217662*10**-19
  k = 1.38*10**-23
  T = np.arange(1,2000,0.1)
  beta = 1/(k*T)
  
  #partitionfunction
  for i in np.arange(1,len(T)): 
      for s in np.arange(1,31,2):
      p[s] = 3*(2*s+1)*math.exp(-(s**2+s)*epsilon*beta(i))
      Zodd[i] = sum(p)
  
  #energy
  ln_Zodd = math.log(Zodd)
  
  for i in np.arange(1,(len(T) - 1)):
      Epara[i]=- (ln_Zodd(i + 1) - ln_Zodd(i)) / (beta(i + 1) - beta(i))
  
  #heat capacity
  for i in np.arange(1,(len(T) - 2)):
      Cpara[i]=(Epara(i + 1) - Epara(i)) / (T(i + 1) - T(i))
  
  #plot
  x = k*T/epsilon
  plt.plot(x(np.arange(1,6000)),Cpara(np.arange(1,6000)) / k,'r')
  plt.axis([0, 7, 0, 1.5]) 
  plt.ylabel('C_v/k')
  plt.xlabel('kT/eps')   
  plt.show()
  </code></pre>
  
  <p>This should be the easiest way to calculate (approximate) this problem because the analytic expression is way more involved. I'm new to python so any suggestions or corrections are appreciated.</p>