#%%
import numpy as np
import sys,os
import math
import matplotlib.pyplot as plt

os.chdir('/Users/tadano/Work/calc/BGG/BTE')

# lattice constant in nm
a = 1.095396826356
# Tau values to read
fildata = 'cell+4/tau50K_SCP.dat'
data = np.loadtxt(fildata, usecols=[0,2,3,4,5])
qind = [int(t) - 1 for t in data[:,0]]

aa = np.array([[a, 0.0, 0.0], [0.0, a, 0.0], [0.0, 0.0, a]])
bb = np.linalg.inv(aa)
factor = 2.0 * np.pi
for i in range(3):
    for j in range(3):
        bb[i,j] = factor * bb[i, j]

q = np.loadtxt('qmesh999.dat')
nq = len(q)
qnorm = np.zeros(nq)
wavelen = np.zeros(nq)

for iq in range(nq):
    qtmp = bb[:,0] * q[iq,0] + bb[:,1] * q[iq,1] + bb[:,2] * q[iq,2]
    qnorm[iq] = math.sqrt(np.dot(qtmp, qtmp))
    if qnorm[iq] > 0.0:
        wavelen[iq] = factor / qnorm[iq]




fig = plt.figure()
ax = plt.subplot(111)
plt.plot(wavelen[qind], data[:,4], linestyle='None', marker='o', ms=3)

x = np.linspace(1, 20, num=10000)
plt.plot(x, x, linestyle='-', lw = 1.0)
ax.set_xscale('log')
ax.set_yscale('log')
plt.show()