#include <cmath>
#include <cassert>
#include <iostream>
#include <vector>

using namespace std;

inline std::pair<double,double> eps(double x, double y, double tp)
{ // returns pair of energies which correspond to two band of graphene
  double a = 2*(cos(x)+cos(y)+cos(x-y));
  double v = sqrt(3+a);
  return std::make_pair(v-tp*a,-v-tp*a);
}

// This class takes care of binning real numbers in the interval[a:b]
// into integer index between [0:Nbin-1]
class Bins{
  int Nbin;
  double a, b;
public:
  Bins(int Nbin_, double a_, double b_) : Nbin(Nbin_), a(a_), b(b_){};
  int to_int(double x) const
  { // gives an integer between [0:Nbin-1] for energies between [a:b]
    //    assert(x>=a && x<b); // checks that energy is inside the interval
    if (!(x>=a && x<b)) cerr<<"Point is outside the interval ["<<a<<","<<b<<"],  "<<x<<endl;
    return static_cast<int>((x-a)/(b-a)*(Nbin)); 
  }
  double to_real(int i) const
  { // gives real number that corresponds to bin i
    return a + (b-a)*(i+0.5)/Nbin;
  }
  double min() const {return a;}
  double max() const {return b;}
};

void grapheneDOS(vector<double>& Dos, int Nq, double tp, const Bins& bin)
{ // computes DOS by calling function eps which returns energy of the band
  int Nbin = Dos.size();
  for (int i=0; i<Nq; i++){
    double x = -M_PI + 2*M_PI*i/(Nq+0.0);   // kx in [-pi,pi]
    for (int j=0; j<Nq; j++){
      double y = -M_PI + 2*M_PI*j/(Nq+0.0); // ky in [-pi,pi]
      pair<double,double> e = eps(x,y,tp);  // energy of the two bands
      // updates DOS
      Dos[bin.to_int(e.first)]++;
      Dos[bin.to_int(e.second)]++;
    }
  }
  
  double norm = Nbin/(bin.max()-bin.min())/(2*Nq*Nq+0.0);
  for (int i=0; i<Dos.size(); i++) Dos[i] *=norm;
}

int main(int argc, char *argv[], char *env[])
{
  int Nq=2000;          // numer of k-points in each dimension
  double tp=0;//-1/12.; // next nearest neighbor hopping
  int Nbin = 200;       // number of beans
  
  int i=0;
  while (++i<argc){
    std::string str(argv[i]);
    if (str=="-Nq" && i<argc-1) Nq = atoi(argv[++i]);
    if (str=="-tp" && i<argc-1) tp = atof(argv[++i]);
    if (str=="-Nbin" && i<argc-1) Nbin = atoi(argv[++i]);
    if (str=="-h" || str=="--help"){
      std::clog<<"********************* DOS for graphene ***************\n";
      std::clog<<"**                                                  **\n";
      std::clog<<"**      Copyright Kristjan Haule, Feb 2007          **\n";
      std::clog<<"******************************************************\n";
      std::clog<<"\n";
      std::clog<<"honey [Options] \n" ;
      std::clog<<"Options:   -Nq         Number of momentum points in each dimension ("<<Nq<<")\n";
      std::clog<<"           -tp         Next neares neighbor hopping ("<<tp<<")\n";
      std::clog<<"           -Nbin       Number of beans for calculation of DOS ("<<Nbin<<")\n";
      std::clog<<"*****************************************************\n";
      return 0;
    }
  }


  // Creates a class which can convert real number into integer index
  // in array, or, can give real number from index in array
  Bins bin(Nbin, -3-6*tp, 3-6*tp+1e-7); // the limits of energy to graphene
  
  vector<double> Dos(Nbin);
  for (int i=0; i<Dos.size(); i++) Dos[i]=0; // vectors from STD unfortunately does not support simple assignment

  // computes DOS
  grapheneDOS(Dos, Nq, tp, bin);

  // printing
  for (int i=0; i<Dos.size(); i++)
    cout<<bin.to_real(i)<<" "<<Dos[i]<<endl;
  
  return 0;
}