/**************************************************************************
*               This script reads in a CPMD 3.4.1 output of               *
*               a bandstructure calculation  and prints the               *
*               eigenvalues in a format readable by graphic               *
*                                 software                                *
*                                                                         *
*                   cc -o bandstructure bandstructure.c                   *
*                                                                         *
*                          Veronika Brazdova 2002                         *
*                          vb@chemie.hu-berlin.de                         *
*                                                                         *
***************************************************************************/

#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <string.h>

#define BUFFER 150

/// Return codes for function main:
//  0 successful run
//  1 input file not found
//  2 input file not in a correct format

int main(int argc, char *argv[])
{
  FILE   *fp;
  double xshift, yshift, kpoint, value1, value2;
  int    i, states;
  int    skip, found;
  char   line[BUFFER];

  // argv[0] = name of the executable
  // argv[1] = input file
  // argv[2] = shift in the x axis (default 0)
  // argv[3] = how many low lying states should be left out (default 0)
  // argv[4] = shift in the Fermi level (default 0)

  /// Handle the command line arguments
  if (argc < 2)
  {
    fprintf(stderr,\
      "bandstructure cpmd_output x_axis_shift Fermi_level_shift no_skipped_states\n"); 
    return 0;
  }

  if (!(fp = fopen(argv[1],"r")))
  {
    fprintf(stderr,"bandstructure: file %s can not be open\n",argv[1]);
    return 1;
  }
  
  if (argc < 3)
    xshift = 0;
  else
    xshift = atof(argv[2]);

  if (argc < 4)
    yshift = 0;
  else
    yshift = atof(argv[3]);

  if (argc < 5)
    skip = 0;
  else
  {
    skip = atoi(argv[4]);
    if (skip%2 == 1)
      skip = skip + 1;
  } 

  found = 0;
  while (fgets(line,BUFFER,fp))
  {
    if (strstr(line," NUMBER OF STATES: ") != NULL)
    {
      sscanf(line," NUMBER OF STATES:    %i  ",&states);
      found = 1;
      break;
    }
  }
  if (!found)
  {
    fprintf(stderr,"Number of states not found.\n");
    return 2;
  }
 
  found = 0;
  if (states%2 == 1)
    states = states/2 +1;
  else 
    states = states/2;

  while (fgets(line,BUFFER,fp))
    if (strstr(line,"FINAL RESULTS"))
       break;
  while (fgets(line,BUFFER,fp)) 
  {
    if (strstr(line,"EIGENVALUES(EV)") != NULL)
    {
      found = 1;
      while (fgets(line,BUFFER,fp))
      {
        if(strstr(line," CHEMICAL") != NULL)
          break;
        sscanf(line," K POINT: %lf %*lf %*lf %*lf %*lf",&kpoint); 
        for (i = 0; i < skip/2; i++)
        {
          fgets(line,BUFFER,fp);
        }
        for (i = 0; i < (states - skip/2); i++)
        {
          fgets(line,BUFFER,fp);
          sscanf(line," %*i %lf %*lf %*i %lf %*lf ",&value1,&value2);
          printf("%6.3lf %15.9lf\n%6.3lf %15.9lf\n",\
            kpoint+xshift,value1+yshift,kpoint+xshift,value2+yshift); 
        }
      }
    }
  }

  if (!found)
  {
    fprintf(stderr,"File %s does not contain the eigenvalues.\n",argv[1]);
    return 2;
  }
  return 0;

} /// main