N_T_Ar_DW = int((T_max_Ar_DW - T_min_Ar_DW) // T_step_Ar_DW + 1)
one_on_L_Ar_DW = 1.0/np.arange(L_min_Ar_DW, L_max_Ar_DW + 1)
coeff_fit_Ar_DW = np.zeros((N_T_Ar_DW,2)) # last two indices: slope and intercept
coeff_err_Ar_DW = np.zeros((N_T_Ar_DW,2)) # last two indices: slope and intercept
K_eff = np.zeros(N_T_Ar_DW)
k_BOLTZ = 1.38064852e-23 # J/K
fctr_DW = 3.0*k_BOLTZ/(alat_Ar_DW*1e-10)/np.pi/1e9*1e20
iLtmp = 1
mrkr = ['o', '^', 'd', 'v', 's']
colore = cm.get_cmap('inferno')
fig, ax1 = plt.subplots()
Temprange = (np.arange(T_min_Ar_DW, T_max_Ar_DW + T_step_Ar_DW, T_step_Ar_DW)*0.1 + 0.2)*119.8
for iT, T in enumerate(Temprange):
# fit and find covariance of the fit parameters:
coeff_fit_Ar_DW[iT], covar_Ar_DW = np.polyfit(one_on_L_Ar_DW[iLtmp:], 0.5*A_self_Ar_DW[iLtmp:,iT,0] / T / fctr_DW, 1, cov=True)
# compute error on fit parameters:
coeff_err_Ar_DW[iT] = np.sqrt(np.diag(covar_Ar_DW))
print('K_eff between',
-1.0/(coeff_fit_Ar_DW[iT,0]+coeff_err_Ar_DW[iT,0]),
' and ',
-1.0/(coeff_fit_Ar_DW[iT,0]-coeff_err_Ar_DW[iT,0]),
'GPa', 'best = ', -1.0/(coeff_fit_Ar_DW[iT,0]), 'GPa')
K_eff[iT] = -1.0/coeff_fit_Ar_DW[iT,0]
#plt.text(0, 0.03+iT*0.007, '$K^\mathrm{eff} = $'+str(np.round(-1.0/coeff_fit_Ar_DW[iT,0],1))+'GPa, @$T=$'+str(T)+'K',
# fontsize=11)
str1 = str("$T =$"+str(np.round(T,2))+"K")
str2 = str('fit $K^\mathrm{eff}$ = '+str(np.round(-1.0/coeff_fit_Ar_DW[iT,0],1))+' GPa')
ax1.errorbar(one_on_L_Ar_DW,
0.5*A_self_Ar_DW[:,iT,0]/T/fctr_DW,
0.5*A_self_Ar_DW[:,iT,1]/T/fctr_DW,
ls='none',
linewidth=1,marker=mrkr[iT], markersize=3.5, color=colore(30+iT*50),
label=str1)
ax1.plot(np.arange(14,51)*0.01,
(np.arange(14,51)*0.01*coeff_fit_Ar_DW[iT,0] + coeff_fit_Ar_DW[iT,1]),
ls='-', linewidth=1,
color=colore(50+iT*50), label=None)
ax1.plot(np.arange(0,11)*0.1,
(np.arange(0,11)*0.1*coeff_fit_Ar_DW[iT,0] + coeff_fit_Ar_DW[iT,1]),
ls=':', linewidth=1,
color=colore(50+iT*50), label=None)
# These are in unitless percentages of the figure size. (0,0 is bottom left)
left, bottom, width, height = [0.62, 0.61, 0.25, 0.25]
ax2 = fig.add_axes([left, bottom, width, height])
print(coeff_fit_Ar_DW[:,0])
print(coeff_err_Ar_DW[:,0])
ax2.errorbar(Temprange,
np.abs(coeff_fit_Ar_DW[:,0]),
coeff_err_Ar_DW[:,0],
ls='-', linewidth=1,
marker='D', color='green', markersize=3)
ax1.set_xlabel('$1/\ell$')
ax1.set_xticks([1,1/2,1/3,1/4,1/5,1/6,1/7,0])
ax1.set_xticklabels(['1','$\\frac{1}{2}$','$\\frac{1}{3}$','$\\frac{1}{4}$','$\\frac{1}{5}$','','... ','0'])
#ax1.set_ylabel('$B / \\frac{3kT}{\pi a}$ [GPa$^{-1}$]')
ax1.set_ylabel('$B / \\frac{8\pi kT}{a}$ [GPa$^{-1}$]')
ax1.legend(frameon=True, loc='lower left', fontsize=11, markerscale=1.0, handletextpad=0.1)
ax1.tick_params(which='both', direction='in')
ax2_font = {'size':'9.5'}
ax2.set_xlabel('$T$ [K]', labelpad=2, **ax2_font)
#ax2.set_xticks([400, 800, 1200])
ax2.set_ylabel('$|m|$ [GPa$^{-1}$]', labelpad=2, **ax2_font)
ax2.tick_params(which='both', direction='in', labelsize=9.5)
#plt.savefig('B_Ar_DW.pdf', format='pdf', bbox_inches = 'tight')
plt.savefig('B_Ar_DW_def_8pi2.pdf', format='pdf', bbox_inches = 'tight')
