import math


def read_wien2k_bands(filename,a,b,c):
    """
    Read the Wien2K band structure data and return
    - List of 1D positions along the path
    - A dictionary with the data. The keys consist of tuples of
      the point number and the band number
    - The minimum number of bands at any point
    - The number of k-points
    """
    numpoints = 0
    minbands  = 1000000
    points    = []
    data      = {}
    line      = ""
    with open(filename,"r") as fp:
        #
        # Skip over initial stuff, looking for the Gamma point
        #
        while line.find("G") == -1:
            line = fp.readline()
            if not line:
                raise Exception("ERROR: Not a valid band structure file: {}".format(filename))
                
        #
        # Read bands
        #
        while True:
            #
            # Read k-point line
            #
            numpoints += 1
            nbands = int(line[74:80])
            minbands = min(minbands,nbands)
            sx = float(line[ 0:19])
            sy = float(line[19:38])
            sz = float(line[38:57])
            points.append((sx,sy,sz))
            #
            # Read band energies
            #
            for ii in range(0,nbands):
                line = fp.readline()
                items = line.split()
                iband = int(items[0])
                energy = float(items[1])
                data[(numpoints,iband)] = energy
            line = fp.readline()
            if not line:
                break
    prev = (0.0,0.0,0.0)
    xcoords = []
    r = 0.0
    for current in points:
        # Current is a point in fractional coordinates (xf,yf,zf)
        # this needs to be converted to Cartesian coordinates (xc,yc,zc)
        (xf,yf,zf) = current
        xc         = xf*a[0]+yf*b[0]+zf*c[0]
        yc         = xf*a[1]+yf*b[1]+zf*c[1]
        zc         = xf*a[2]+yf*b[2]+zf*c[2]
        (xp,yp,zp) = prev
        current_c  = (xc,yc,zc)
        d=(xc-xp)*(xc-xp)+(yc-yp)*(yc-yp)+(zc-zp)*(zc-zp)
        d=math.sqrt(d)
        r+=d
        xcoords.append(r)
        prev = current_c
    return (xcoords,data,minbands,numpoints)