#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
#include <stdio.h>
#include <gsl/gsl_math.h>
#define min( a, b ) ( ((a) < (b)) ? (a) : (b) )

double AOrB(double qval,const double aorb){return cos(aorb*qval);}

void DEntries(gsl_matrix *DMatrix,int pnum){
 int i;
 double Cx,Cy,Sx,Sy;
 const double sqrt3on2=sqrt(3.0)/2.0,sqrt3=sqrt(3.0);
 
 for (i=0;i<pnum;i++){
   Cx=cos(gsl_matrix_get(DMatrix,i,0)/2.0);Cy=cos(gsl_matrix_get(DMatrix,i,1)*sqrt3on2);
   Sx=sin(gsl_matrix_get(DMatrix,i,0)/2.0);Sy=sin(gsl_matrix_get(DMatrix,i,1)*sqrt3on2);
   
   gsl_matrix_set(DMatrix,i,2,4.0*(2.0-Cx*(Cx+Cy)));gsl_matrix_set(DMatrix,i,3,3.0*(1.0-Cx*Cy));
   gsl_matrix_set(DMatrix,i,4,2.0*(1.0+Cx*(Cy-2.0*Cx) ));gsl_matrix_set(DMatrix,i,5,sqrt3*Sx*Sy); 
 }
}

void SortMatrix(gsl_matrix *DMatrix, int pnum){
  int dummy=1,i;
  double a,b,c,d,e,f;
  
  while(dummy==1){
    dummy=0;
    for(i=0;i<pnum-1;i++){
      if(gsl_matrix_get(DMatrix,i+1,0)>gsl_matrix_get(DMatrix,i,0)){
        a=gsl_matrix_get(DMatrix,i+1,0);b=gsl_matrix_get(DMatrix,i+1,1);
        c=gsl_matrix_get(DMatrix,i+1,2);d=gsl_matrix_get(DMatrix,i+1,3);
        e=gsl_matrix_get(DMatrix,i+1,4);f=gsl_matrix_get(DMatrix,i+1,5);
        
        
        gsl_matrix_set( DMatrix,i+1,0,gsl_matrix_get(DMatrix,i,0) );gsl_matrix_set(DMatrix,i+1,1,gsl_matrix_get(DMatrix,i,1) );
        gsl_matrix_set( DMatrix,i+1,2,gsl_matrix_get(DMatrix,i,2) );gsl_matrix_set(DMatrix,i+1,3,gsl_matrix_get(DMatrix,i,3) );
        gsl_matrix_set( DMatrix,i+1,4,gsl_matrix_get(DMatrix,i,4) );gsl_matrix_set(DMatrix,i+1,5,gsl_matrix_get(DMatrix,i,5) );
        
        gsl_matrix_set( DMatrix,i,0, a);gsl_matrix_set( DMatrix,i,1,b );
        gsl_matrix_set( DMatrix,i,2, c);gsl_matrix_set( DMatrix,i,3,d );
        gsl_matrix_set( DMatrix,i,4, e);gsl_matrix_set( DMatrix,i,5,f );
        
        dummy=1;
      }
    }
  }
  
  dummy=1;
  while(dummy==1){
    dummy=0;
    for(i=0;i<pnum-1;i++){
      if(gsl_matrix_get(DMatrix,i+1,1)>gsl_matrix_get(DMatrix,i,1) && gsl_matrix_get(DMatrix,i+1,0)==gsl_matrix_get(DMatrix,i,0)){
        a=gsl_matrix_get(DMatrix,i+1,0);b=gsl_matrix_get(DMatrix,i+1,1);
        c=gsl_matrix_get(DMatrix,i+1,2);d=gsl_matrix_get(DMatrix,i+1,3);
        e=gsl_matrix_get(DMatrix,i+1,4);f=gsl_matrix_get(DMatrix,i+1,5);
        
        
        gsl_matrix_set( DMatrix,i+1,0,gsl_matrix_get(DMatrix,i,0) );gsl_matrix_set(DMatrix,i+1,1,gsl_matrix_get(DMatrix,i,1) );
        gsl_matrix_set( DMatrix,i+1,2,gsl_matrix_get(DMatrix,i,2) );gsl_matrix_set(DMatrix,i+1,3,gsl_matrix_get(DMatrix,i,3) );
        gsl_matrix_set( DMatrix,i+1,4,gsl_matrix_get(DMatrix,i,4) );gsl_matrix_set(DMatrix,i+1,5,gsl_matrix_get(DMatrix,i,5) );
        
        gsl_matrix_set( DMatrix,i,0, a);gsl_matrix_set( DMatrix,i,1,b );
        gsl_matrix_set( DMatrix,i,2, c);gsl_matrix_set( DMatrix,i,3,d );
        gsl_matrix_set( DMatrix,i,4, e);gsl_matrix_set( DMatrix,i,5,f );
        dummy=1;
      }
    }
  }
}

void OutputMatrix(gsl_matrix *DMatrix,int pnum,char fname[1024]){
  int i;
  FILE *fp=fopen(fname,"w");fprintf(fp,"#qx(1) qy(2) Dg00(3) Dh00(4) D11(5) D12(6)\n");
  for (i=0;i<pnum;i++){fprintf(fp,"%f %f %f %f %f %f\n",gsl_matrix_get(DMatrix,i,0),gsl_matrix_get(DMatrix,i,1),gsl_matrix_get(DMatrix,i,2),gsl_matrix_get(DMatrix,i,3),gsl_matrix_get(DMatrix,i,4),gsl_matrix_get(DMatrix,i,5));}
  fclose(fp);
}

void ConstructDmatrix(int pnum,char fname[1024]){
 int sqrtpnum=sqrt(pnum),d1=0,d2=0,t;
 const double oneonsqrt3=1.0/sqrt(3.0),latticeC=4.0*M_PI/sqrt(3.0);
 const double qx=2.0*M_PI/sqrtpnum,qy=latticeC/sqrtpnum,Emargin=0.001*latticeC/sqrtpnum;
 double y,x,a,b;
 gsl_matrix *DMatrix=gsl_matrix_alloc(pnum,6);

 
 //The matrix is initialized with an absurd value for qx and qy so that it is obvious if one of its is not treated correctly
 for(x=0;x<pnum;x++){gsl_matrix_set(DMatrix,(int)x,0,-1000.0);gsl_matrix_set(DMatrix,(int)x,1,-1000.0);}


//Draw interior of BZ:  /////
 for(y=qy;y<latticeC/2.0-Emargin*qy;y+=qy){
   for(x=qx;x<(latticeC-y)*oneonsqrt3-Emargin*qx;x+=qx){
     gsl_matrix_set(DMatrix,4*d2  ,0, x);gsl_matrix_set(DMatrix,4*d2  ,1, y);
     gsl_matrix_set(DMatrix,4*d2+1,0,-x);gsl_matrix_set(DMatrix,4*d2+1,1, y);
     gsl_matrix_set(DMatrix,4*d2+2,0, x);gsl_matrix_set(DMatrix,4*d2+2,1,-y);
     gsl_matrix_set(DMatrix,4*d2+3,0,-x);gsl_matrix_set(DMatrix,4*d2+3,1,-y);
     
//     printf("%f %f\n",x,y);
//     printf("%f %f\n",-x,y);
//     printf("%f %f\n",x,-y);
//     printf("%f %f\n",-x,-y);
    d2++;
  }  
 }
 d1=4*d2;
 
 //Here we have the horizontal part /////
 d2=0;y=0.0;
 for(x=qx;x<latticeC*oneonsqrt3-Emargin*qx;x+=qx){
   gsl_matrix_set(DMatrix,d1+2*d2  ,0, x);gsl_matrix_set(DMatrix,d1+2*d2  ,1, y);
   gsl_matrix_set(DMatrix,d1+2*d2+1,0,-x);gsl_matrix_set(DMatrix,d1+2*d2+1,1, y);
   
//    printf("%f %f\n",x,y);
//    printf("%f %f\n",-x,y);
   d2++;
 }
 d1+=2*d2;
 
 //...and finally the vertical part. /////
 d2=0;x=0.0;
 for(y=-latticeC/2.0+qy;y<latticeC/2.0-Emargin*qy;y+=qy){
   gsl_matrix_set(DMatrix,d1+d2,0, x);gsl_matrix_set(DMatrix,d1+d2,1, y);
//    printf("%f %f\n",x,y);
   d2++;
 }
 d1+=d2;

//Now for the boundaries... /////
 //The top boundary, /////
 d2=0;
 b=latticeC/2.0;x=1.0;y=qx;
 gsl_matrix_set(DMatrix,d1,0,0.0);gsl_matrix_set(DMatrix,d1,1, b);d1++;
//  printf("%f %f\n",a,b);
  
 while(y<2.0*M_PI/3.0+Emargin*qx){
   gsl_matrix_set(DMatrix,d1+d2  ,0, qx*x);gsl_matrix_set(DMatrix,d1+d2  ,1, b);
   gsl_matrix_set(DMatrix,d1+d2+1,0,-qx*x);gsl_matrix_set(DMatrix,d1+d2+1,1, b);
//    printf("%f %f\n",a+qx*x,b);
//    printf("%f %f\n",a-qx*x,b);
   d2+=2;x++;y+=qx;
 }
 d1+=d2;
 
//the diagonal boundaries /////

  d2=0;x=0.0;
 while(x<(oneonsqrt3*latticeC/2.0)+Emargin*qx){x+=qx;}
 
 for(x=x;x<latticeC*oneonsqrt3-Emargin*qx;x+=qx){
   t=(int)((latticeC-x/oneonsqrt3)/qy+Emargin);
   if( t!=0 && fabs(t*qy-(latticeC-x/oneonsqrt3))<Emargin){
     gsl_matrix_set(DMatrix,d1+2*d2  ,0,x);gsl_matrix_set(DMatrix,d1+2*d2  ,1, t*qy);
     gsl_matrix_set(DMatrix,d1+2*d2+1,0,x);gsl_matrix_set(DMatrix,d1+2*d2+1,1,-t*qy);
//      printf("%f %f\n%f %f\n",x,t*qy,x,-t*qy);
     d2++;
   } 
 }
 d1+=2*d2;

 DEntries(DMatrix,pnum);
 SortMatrix(DMatrix,pnum);
 OutputMatrix(DMatrix,pnum,fname);
}

void main(int argc,char **argv){
  int pnum;
  char fname[1024];
  if(argc<2){printf("Give me a particle Number!\n");exit(1);}
  sscanf(argv[1],"%d",&pnum);
  if( pow((int)sqrt(pnum),2)!=pnum){printf("Are you sure you gave me a square number?\n");exit(1);}
  sprintf(fname,"Dvals_pnum_%d.dat",pnum);
  ConstructDmatrix(pnum,fname);
}