In [1]:

from sys import path
path.append("/home/federico/Softwares/sportran/thermocepstrum/")
import numpy as np
import thermocepstrum as tc
import matplotlib.pyplot as plt

In [2]:

DT = 14.6 # femtoseconds

1. Load time series of the currents

In [3]:

elc = np.loadtxt('elclass.dat', skiprows=1)
elq = np.loadtxt('elquant.dat', skiprows=1)
timeseries = DT*np.arange(len(elc[:,0])) # femtoseconds

2. Create objects HeatCurrent

In [4]:

jjq = tc.HeatCurrent(elq, 'metal', DT, 1., 1., 0.3)
jjq.kappa_scale = 1.
jjc = tc.HeatCurrent(elc, 'metal', DT, 1., 1., 0.3)
jjc.kappa_scale = 1.

Using single component code.
Using single component code.

3. Perform cepstral analysis (caveat: units are not W/mK as stated, but $ang^2/fs$)

In [5]:

jjq.cepstral_analysis()
jjc.cepstral_analysis()

-----------------------------------------------------
  CEPSTRAL ANALYSIS
-----------------------------------------------------
  AIC_Kmin  = 13  (P* = 14, corr_factor = 1.000000)
  L_0*   =          -2.485457 +/-   0.059619
  S_0*   =           0.120455 +/-   0.007181
-----------------------------------------------------
  kappa* =           0.060227 +/-   0.003591  W/mK
-----------------------------------------------------

-----------------------------------------------------
  CEPSTRAL ANALYSIS
-----------------------------------------------------
  AIC_Kmin  = 12  (P* = 13, corr_factor = 1.000000)
  L_0*   =          -2.485301 +/-   0.057368
  S_0*   =           0.120473 +/-   0.006911
-----------------------------------------------------
  kappa* =           0.060237 +/-   0.003456  W/mK
-----------------------------------------------------

In [6]:

plt.rcParams.update({'font.size': 12})

In [7]:

plt.plot(jjq.dct.logpsdK, '.-', c='tab:blue', label='$C_n$')
plt.plot(np.arange(len(jjq.dct.logpsdK)), 0.*(np.arange(len(jjq.dct.logpsdK))),'--', color='tab:grey')
plt.fill_between(np.arange(len(jjq.dct.logpsdK)), jjq.dct.logpsdK - jjq.dct.logpsdK_THEORY_std, jjq.dct.logpsdK + jjq.dct.logpsdK_THEORY_std, color='tab:blue', alpha=0.4)
plt.axvline(x = jjq.dct.aic_Kmin + 1, ls='--', c='tab:red')
plt.xlim([0,40])
plt.ylim([-.2,+.2])
plt.xticks([0,10,20,30,40])
plt.yticks([-0.2,-0.1,0,0.1,0.2])
plt.ylabel('$C_n$')
plt.xlabel('$n$')
plt.tick_params('both', bottom=True, top=False, left=True, right=False, direction='in')
plt.savefig('Cn.pdf', bbox_inches = 'tight', pad_inches = 0)
print("Akaike Information Criterion optimal P*:  ", jjq.dct.aic_Kmin + 1)

Akaike Information Criterion optimal P*:   14

No description has been provided for this image

In [8]:

plt.plot(np.arange(jjq.Nfreqs) + 1, jjq.dct.logtau, '.-', c='tab:blue', label='$L_0(P^\ast)$')
plt.fill_between(np.arange(jjq.Nfreqs) + 1, jjq.dct.logtau - jjq.dct.logtau_THEORY_std, jjq.dct.logtau + jjq.dct.logtau_THEORY_std, color='tab:blue', alpha=0.4)
plt.axvline(x = jjq.dct.aic_Kmin + 1, ls='--', c='tab:red')
plt.xlim([0,40])
plt.ylim([-2.9,-1.7])
plt.xticks([0,10,20,30,40])
plt.yticks([-2.9, -2.3, -1.7])
plt.xlabel('$P^\\ast$')
plt.ylabel('$L_0^\\ast$')
plt.tick_params('both', bottom=True, top=False, left=True, right=False, direction='in')
plt.savefig('L0ast.pdf', bbox_inches = 'tight', pad_inches = 0)