#! /usr/bin/env python
# ===================================================================
# file:   auto.py
# author: k.koepernik@ifw-dresden.de
# date:   19 Jun 2017
from __future__ import print_function
import sys
import numpy as np
import numpy.linalg as LA

import pyfplo.slabify as sla

print( '\npyfplo version=: {0}\nfrom: {1}\n'.format(sla.version,sla.__file__))




# ===================================================================
# 
# ===================================================================
def work():


    hamdata='../../+hamdata'
    s=sla.Slabify()
    s.object='3d'
    s.prepare(hamdata)

    
    homo=56

    nk=10
    tol=1e-4
    
    # setup primitive reciprocal cell BoxMesh 
    G=s.hamdataRCell()
    print( 'G=\n',G)
    x=G[:,0]
    y=G[:,1]
    z=G[:,2]
    dx=LA.norm(x)
    dy=LA.norm(y)
    dz=LA.norm(z)
    # with most isotropic mesh subdivision
    nx=max(2,nk)
    ny=max(2,int(dy/dx*nx))
    nz=max(2,int(dz/dx*nx))
    
    
    # shift the origin by a small amount
    k0=-(x/(nx-1)+y/(ny-1)+z/(nz-1))*0.5

    # define the BoxMesh
    box=sla.BoxMesh()
    box.setBox(xaxis=x,yaxis=y,zaxis=z,origin=k0)
    box.setMesh(nx=nx,xinterval=[-dx/2,dx/2],
                ny=ny,yinterval=[-dy/2,dy/2],
                nz=nz,zinterval=[-dz/2,dz/2])

    # and try to find Weyl points
    s.findWeylPoints(box,[homo],tol)

    try:
        from weylpoints import wps
    except:
        raise RuntimeError('file weylpoints.py does not exist')

    # the conventional cell
    A=s.hamdataCCell()
    # its lattice parameters
    a=LA.norm(A[:,0])
    b=LA.norm(A[:,1])
    c=LA.norm(A[:,2])
    boa=b/a
    coa=c/a
    # print in units of 2p/a_i
    for i,w in enumerate(wps):
        print( 'WP{0:4d}: {1:10.6f} {2:10.6f} {3:10.6f} {4:5.1f} {5:10.6f} {6:10d} {7}'\
            .format(i,w.k[0],w.k[1]*boa,w.k[2]*coa,w.chirality,w.energy,
                    w.homo,w.spin))
    
    return


# ===================================================================
# 
# ===================================================================

if __name__ == '__main__':
    work()
