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

src/2_debye_model.md

@@ -101,34 +101,42 @@ $$
 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
-def integrand(y):
-    return y**4 * np.exp(y) / (np.exp(y) - 1)**2
+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=1):
+def c_einstein(T, T_E):
     x = T_E / T
     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
+
 @np.vectorize
-def c_debye(T, T_D=1):
+def c_debye(T, T_D):
     x = T / T_D
     return 9 * x**3 * quad(integrand, 0, 1/x)[0]
 
-T = np.linspace(0.01, 1.5, 500)
-fig, ax = pyplot.subplots()
+temp = np.linspace(1, 215, 100)
+
+fit = curve_fit(c_einstein, T, c, 500)
+T_E = fit[0][0]
 
-ax.plot(T, c_einstein(T), label="Einstein model")
-ax.plot(T, c_debye(T), label="Debye model")
-
-ax.set_ylim(bottom=0, top=3.4)
-ax.set_xlabel('$T$')
-ax.set_ylabel(r'$\omega$')
-ax.set_xticks([1])
-ax.set_xticklabels([r'$\Theta_D$'])
-ax.set_yticks([3])
-ax.set_yticklabels(['$3Nk_B$'])
-ax.legend(loc='lower right')
-pyplot.hlines([3], 0, 1.5, linestyles='dashed')
-draw_classic_axes(ax, xlabeloffset=0.3)
+fit = curve_fit(c_debye, T, c, 500)
+T_D = fit[0][0]
+
+fig, ax = pyplot.subplots()
+ax.scatter(T, c)
+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}K$')
+ax.plot(temp, c_einstein(temp, T_E), label=f'Einstein model, $T_E={T_E:.5}K$')
+ax.set_ylim(bottom=0, top=3)
+ax.set_xlim(0, 215)
+ax.set_xlabel('$T(K)$')
+ax.set_ylabel(r'$C/k_B$')
+ax.legend(loc='lower right');
 ```
 
 ## Exercises