Update figure comparing Debye and Einstein models. Added experimental data for silver.

src/2_debye_model.md

@@ -101,13 +101,15 @@ $$
 where $x=\frac{\hbar\omega}{k_{\rm B}T}$ and $\Theta_{\rm D}\equiv\frac{\hbar\omega_{\rm D}}{k_{\rm B}}$, the _Debye temperature_.
 
 ```python
+pyplot.rcParams['axes.titlepad'] = 20 
 
 T = np.array([1.35,2.,3.,4.,5.,6.,7.,8.,10.,12.,14.,16.,20.,28.56,36.16,47.09,55.88,65.19,74.56,83.91,103.14,124.2,144.38,166.78,190.17,205.3])
 c = np.array([0.,0.,0.,0.,0.,0.,0.0719648,0.1075288,0.2100368,0.364008,0.573208,0.866088,1.648496,4.242576,7.07096,10.8784,13.47248,15.60632,17.27992,18.6188,20.33424,21.63128,22.46808,23.05384,23.47224,23.68144])
+c *= 3/24.945 #24.954 is 3Nk_B
 
 def c_einstein(T, T_E):
     x = T_E / T
-    return 24.945 * x**2 * np.exp(x) / (np.exp(x) - 1)**2
+    return 3 * x**2 * np.exp(x) / (np.exp(x) - 1)**2
 
 def integrand(y):
     return y**4 * np.exp(y) / (np.exp(y) - 1)**2
@@ -115,24 +117,30 @@ def integrand(y):
 @np.vectorize
 def c_debye(T, T_D):
     x = T / T_D
-    return 24.945 * 3 * x**3 * quad(integrand, 0, 1/x)[0]
+    return 9 * x**3 * quad(integrand, 0, 1/x)[0]
 
-fig, ax = pyplot.subplots()
+temp = np.linspace(1, 215, 100)
+
+fit = curve_fit(c_einstein, T, c, 500)
+T_E = fit[0][0]
+#delta_T_E = np.sqrt(fit[1][0, 0])
+#print(f"T_E = {T_E:.5} ± {delta_T_E:.3} K")
 
+fit = curve_fit(c_debye, T, c, 500)
+T_D = fit[0][0]
+#delta_T_D = np.sqrt(fit[1][0, 0])
+#print(f"T_D = {T_D:.5} ± {delta_T_D:.3} K")
+
+fig, ax = pyplot.subplots()
 ax.scatter(T, c)
-#ax.set_title('Heat capacity of silver compared to the Debye and Einstein models')
-ax.plot(T, c_einstein(T, 151), label='Einstein model')
-ax.plot(T, c_debye(T, 215), label='Debye model')
-ax.set_ylim(bottom=0, top=26)
+ax.set_title('Heat capacity of silver compared to the Debye and Einstein models')
+ax.plot(temp, c_debye(temp, T_D), label=f'Debye model, $T_D={T_D:.5}$')
+ax.plot(temp, c_einstein(temp, T_E), label=f'Einstein model, $T_E={T_E:.5}$')
+ax.set_ylim(bottom=0, top=3)
 ax.set_xlim(0, 215)
-ax.set_xlabel('T(K)')
-ax.set_ylabel(r'C(J/mol-K)')
-ax.set_xticks([0, 100, 200])
-ax.set_yticks([24.945])
-ax.set_yticklabels(['$3R$'])
-ax.legend(loc='lower right')
-pyplot.hlines([24.945], 0, T[-1], linestyles='dashed')
-draw_classic_axes(ax, xlabeloffset=0.3)
+ax.set_xlabel('$T(K)$')
+ax.set_ylabel(r'$C/k_B$')
+ax.legend(loc='lower right');
 ```
 
 ## Exercises