title: Specific heat fit using Debye+Einstein model in matplotlib
body: <p>I am trying to find a fit to a specific heat data using gamma<em>T+m</em>Debye_model+(1-m)*Einstein model as given below.</p>
  
  <p>Cel+ph(T ) = γ T + [αCDebye(T ) + (1 − α)CEinstein(T )]</p>
  
  <p>where the Debye and Einstein models are given by eq. 3 and 4 in the attachment.</p>
  
  <p>I have tried the following code in jupyter notebook following some examples on the web but i have no idea how can i combine these functions together to carry out the fit.</p>
  
  <p>The data is linked <a href="https://www.dropbox.com/s/u0r2m3zwl8w77at/HC_ScPtBi.dat?dl=0" rel="nofollow noreferrer">https://www.dropbox.com/s/u0r2m3zwl8w77at/HC_ScPtBi.dat?dl=0</a>
  Column 1 is Temperature and Column 3 is Y data of interest.
  Model is in <a href="https://www.dropbox.com/s/9452fq7eydajr5o/Debye.pdf?dl=0" rel="nofollow noreferrer">https://www.dropbox.com/s/9452fq7eydajr5o/Debye.pdf?dl=0</a></p>
  
  <p>Code is in <a href="https://www.dropbox.com/s/hk9b1t0agvt36zn/Untitled2.ipynb?dl=0" rel="nofollow noreferrer">https://www.dropbox.com/s/hk9b1t0agvt36zn/Untitled2.ipynb?dl=0</a></p>
  
  <pre><code>from matplotlib import pyplot 
  import numpy as np
  from scipy import integrate
  from scipy.optimize import curve_fit
  from scipy.integrate import quad
  
  data=np.genfromtxt('HC_ScPtBi.dat', skip_header=1)
  R=8.314
  n=3
  M=1
  T=data[10:290,0]
  c=data[10:290,2]
  def plot_data():
      pyplot.scatter(T, c)
      pyplot.xlabel('$T [K]$')
      pyplot.ylabel('$C$')
  
  plot_data()
  def c_einstein(T, T_E):
      x = T_E / T
      return 3 *n*R*x**2 * np.exp(x) / (np.exp(x) - 1)**2
  
  popt0, pcov0 = curve_fit(c_einstein, T, c, 250)
  T_E = popt0[0]
  delta_T_E = np.sqrt(pcov0[0, 0])
  print(f"T_E = {T_E:.5} ± {delta_T_E:.3} K")
  print(popt0)
  
  plot_data()
  #temps = np.linspace(10, T[-1], 100)
  pyplot.plot(T, c_einstein(T, *popt0));
  
  def integrand(y):
      return y**4 * np.exp(y) / (np.exp(y) - 1)**2
  
  @np.vectorize
  def c_debye(T, T_D):
      x = T / T_D
      return 9 *n*R*x**3 * quad(integrand, 0, 1/x)[0]
  
  popt1, pcov1 = curve_fit(c_debye, T, c, 150)
  T_D = popt1[0]
  delta_T_D = np.sqrt(pcov1[0, 0])
  print(f"T_D = {T_D:.5} ± {delta_T_D:.3} K")
  print(popt1)
  plot_data()
  pyplot.plot(T, c_einstein(T, *popt0), label='Einstein')
  pyplot.plot(T, c_debye(T,  *popt1), label='Debye')
  pyplot.legend();
  </code></pre>