/* Write a density of states */
/* (c) MJR 2024, 2025 */

/* syntax of command line option:
 *
 * --dos_{fmt}{={smear}}{:npoints}{afgmr}
 *
 * accepted values of fmt: gnu dos raw
 */

/* Plotting routines expect bin_width=(emax-emin)/(nbins-1)
 *                          x=emin+i*bin_width; i=0,nbins-1
 *
 * So for bin-based routines, bin i extends from
 * emin+(i-0.5)*bin_width to emin+(i+0.5)*bin_width
 */

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

#include "c2xsf.h"

#define PLT_MIRROR 1
#define PLT_ROTATE 2
#define PLT_PH 4

static void gauss_dos(struct kpts *k, struct es *e,
		      double e_min, double e_max, double smear,
		      double *bins, int nbins, int gaps);

static void dos_plot_gnu(FILE *outfile, double e_min, double e_max,
			 double *bins, int nbins,
			 struct es *e, char *title, int plot_flags);

static void dos_plot_eps(FILE *outfile, double e_min, double e_max,
			 double *bins, int nbins,
			 struct es *e, char *title, int plot_flags);
static void dos_plot_raw(FILE *outfile, double e_min, double e_max,
			 double *bins, int nbins,
			 struct es *e, char *title, int plot_flags);

int blochl(struct unit_cell *c, struct contents *m,
	    struct kpts *kp, struct es *elect, struct symmetry *rs,
	    double e_min, double e_max, double *y, int n);


void dos_plot(FILE* outfile, struct unit_cell *c, struct contents *m,
	      struct kpts *k, struct symmetry *s, struct es *e,
	      char *fmt, double *range){
  int i,ns,nbins,plot_flags,gaps;
  int add,nspins_saved,fix,tet;
  double x,e_min,e_max,*bins,bin_width,smear;
  char *cptr,*fmt2;

  if ((!e->eval)||(k->n*e->nbands*e->nspins==0)) {
    fprintf(stderr,"No eigenvalues to write %d %d!\n",e->nbands,e->nspins);
    return;
  }

  plot_flags=0;
  add=0;
  gaps=0;
  tet=1;
  fix=0;
  smear=-1;
  nbins=-1;
  if (!strncmp(fmt,"--dos_",6))
    cptr=fmt+6;
  else if (!strncmp(fmt,"--phdos_",8)){
    plot_flags=PLT_PH;
    cptr=fmt+8;
  }
  else
    error_exit("Unexpected format string in dos_plot");

  fmt2=cptr;
  if ((strncmp(cptr,"gnu",3))&&(strncmp(cptr,"eps",3))&&
      (strncmp(cptr,"raw",3)))
    error_exit("Unexpected output format in dos_plot");
  cptr+=3;
  if (*cptr=='='){
    cptr++;
    i=0;
    sscanf(cptr,"%lf%n",&smear,&i);
    cptr+=i;
  }
  if (*cptr==':'){
    cptr++;
    i=0;
    sscanf(cptr,"%d%n",&nbins,&i);
    cptr+=i;
  }
  while (*cptr){
    if (*cptr=='m') plot_flags|=PLT_MIRROR;
    if (*cptr=='r') plot_flags|=PLT_ROTATE;
    if (*cptr=='g') gaps=1;
    if (*cptr=='a') add=1;
    if (*cptr=='f') fix=1;
    if (*cptr=='G') tet=0;
    cptr++;
  }

  if (k->n<8) {
    fprintf(stderr,"Too few k-points for Bloechl method\n");
    tet=0;
  }
  
  if (tet==1) fix=1;
  if (smear<0) smear=0.2;
  if (nbins<=0) nbins=200;
  if (debug){
    if (tet)
      fprintf(stderr,
	      "Calculating DOS using Bloechl's method and %d points\n",nbins);
    else
      fprintf(stderr,
	      "Calculating DOS using %d bins and a smearing width of %feV\n",
	      nbins,smear);
    fprintf(stderr,"%d bands, %d spin%s and %d k-points\n",e->nbands,
	    e->nspins,(e->nspins==1)?"":"s",k->n);
    if ((tet==0)&&(gaps)) fprintf(stderr,"Gaps preserved\n");
  }

  if (range){
    e_min=*range;
    e_max=*(range+1);
  }
  else{
    e_min=1e99;
    e_max=-e_min;

    for(i=0;i<k->n*e->nbands*e->nspins;i++){
      e_min=min(e_min,e->eval[i]);
      e_max=max(e_max,e->eval[i]);
    }
    x=(e_max-e_min)/(2*nbins);
    if (!fix){
      e_min=e_min-x-3*smear;
      e_max=e_max+x+3*smear;
    }
    else{
      e_min=e_min-4*x;
      e_max=e_max+4*x;
    }
  }
  
  bins=malloc(e->nspins*nbins*sizeof(double));
  if (!bins) error_exit("malloc error for bins");

  if (tet){
    if (blochl(c,m,k,e,s,e_min,e_max,bins,nbins)){
      tet=0;
      fprintf(stderr,"Warning: Bloechl failed, falling back to Gaussians\n");
    }
  }
  if (!tet)
    gauss_dos(k,e,e_min,e_max,smear,bins,nbins,gaps);
  
  /* Are bands doubly occupied? */
  if ((e->nspins==1)&&(e->nspinors==1))
    for(i=0;i<nbins;i++) bins[i]*=2;

  /* Unit conversion */
  if (flags&AU){
    e_max/=H_eV;
    e_min/=H_eV;
    for(i=0;i<e->nspins*nbins;i++) bins[i]*=H_eV;
  }

  bin_width=(e_max-e_min)/(nbins-1);
  for(ns=0;ns<e->nspins;ns++){
    x=0;
    for(i=0;i<nbins;i++) x+=bins[ns*nbins+i];
    fprintf(stderr,"DOS sum,      spin=%d: %f\n",ns,x*bin_width);
    if ((nbins>=4)&&((nbins&1)==0)){
      x=bins[ns*nbins]+bins[ns*nbins+nbins-1];
      for(i=1;i<nbins-1;i+=2)
	x+=4*bins[ns*nbins+i]+2*bins[ns*nbins+i+1];
      x/=3;
      fprintf(stderr,"DOS integral, spin=%d: %f\n",ns,x*bin_width);
    }
  }    

  nspins_saved=e->nspins;
  if (add){
    if (debug) fprintf(stderr,"Combining spin densities\n");
    for(i=0;i<nbins;i++)
      for(ns=1;ns<e->nspins;ns++)
	bins[i]+=bins[ns*nbins+i];
    e->nspins=1;
  }

  
  if (!strncmp(fmt2,"gnu",3))
    dos_plot_gnu(outfile,e_min,e_max,bins,nbins,e,m->title,plot_flags);
  else if (!strncmp(fmt2,"eps",3))
    dos_plot_eps(outfile,e_min,e_max,bins,nbins,e,m->title,plot_flags);
  else if (!strncmp(fmt2,"raw",3))
    dos_plot_raw(outfile,e_min,e_max,bins,nbins,e,m->title,plot_flags);
  else
    error_exit("Unexpected format - internal error");

  e->nspins=nspins_saved;
  free(bins);
}

static void dos_plot_gnu(FILE *outfile, double e_min, double e_max,
			 double *bins, int nbins,
			 struct es *e, char *title, int plot_flags){
  double bin_width,scale;
  int ns,i;
  char *e_unit;

  if (flags&AU){
    e_unit="Ha";
    scale=1/H_eV;
  }
  else{
    e_unit="eV";
    scale=1;
  }

  if (nbins<2) error_exit("Must have at least two points for a plot");
  
  bin_width=(e_max-e_min)/(nbins-1);

  if (title) fprintf(outfile,"set title \"%s\"\n",title);
  if (plot_flags&PLT_PH){
    fprintf(outfile,"set xlabel '{/Symbol w} (cm^{-1})'\n");
    fprintf(outfile,"set ylabel 'g({/Symbol w})'\n");
  }
  else{
    fprintf(outfile,"set xlabel 'E(%s)'\n",e_unit);
    fprintf(outfile,"set ylabel 'g(E)'\n");
  }
  if (e->e_fermi)
    fprintf(outfile,
	    "set arrow from %f, graph 0 to %f, graph 1 nohead dt '.'\n"
	    "set link x2\n"
	    "set format x2 ''\n"
	    "set x2tics add right offset character -.4,%.1f ('E_F' %f)\n",
	    *e->e_fermi*scale,*e->e_fermi*scale,
	    ((e->nspins>1)&&(plot_flags&&PLT_MIRROR))?0:-1.7,
	    *e->e_fermi*scale);
  if ((e->nspins>1)&&(plot_flags&&PLT_MIRROR)){
    fprintf(outfile,
	    "set arrow from graph 0, first 0 to graph 1, first 0 nohead\n");
    fprintf(outfile,"set link x2\nset x2tics axis format ''\n");
  }

  fprintf(outfile,"plot ");
  for(ns=0;ns<e->nspins;ns++){
    fprintf(outfile,"'-' notitle with line");
    if (ns+1<e->nspins) fprintf(outfile,",");
  }
  fprintf(outfile,"\n");

  for(ns=0;ns<e->nspins;ns++){
    if ((ns&1)&&(plot_flags&&PLT_MIRROR))
      for(i=0;i<nbins;i++)
	fprintf(outfile,"%f %f\n",e_min+i*bin_width,-bins[ns*nbins+i]);
    else
      for(i=0;i<nbins;i++)
	fprintf(outfile,"%f %f\n",e_min+i*bin_width,bins[ns*nbins+i]);
    fprintf(outfile,"end\n");
  }
}

double tic(double range){
  double sc,tmp,out;

  out=0.1*range;
  sc=pow(10,floor(log10(out)));
  tmp=ceil(out/sc); /* single digit */
  if ((tmp==3)||(tmp==4)) tmp=5;
  if (tmp>=6) tmp=10;
  return tmp*sc;
}

void yaxis_eps(FILE *outfile, double ymin, double ymax,
		      double ytic, int pt_h, char *title, int flag){
  double tmp,y;
  int i;

  fprintf(outfile,"\n%% y-axis\n");
  for(i=0;i<=0.5+(ymax-ymin)/ytic;i++){
    tmp=pt_h*(i*ytic)/(ymax-ymin);
    y=ymin+i*ytic;
    if (flag) y=fabs(y);
    fprintf(outfile,"%d %.1f moveto (%g) rshow\n",-10,tmp-5,y);
    fprintf(outfile,"0 %.1f moveto -5 0 rlineto stroke\n",tmp);
  }
  fprintf(outfile,"-35 %.1f moveto\n",0.5*pt_h);
  if (title)
    fprintf(outfile,"gsave 90 rotate (%s) ctrshow grestore\n",title);
}

static void xaxis_eps(FILE *outfile,double xmin, double xmax,
		      double xtic, int pt_w, char *title, int flag){
  double tmp,x;
  int i;

  fprintf(outfile,"\n%% x-axis\n");
  for(i=0;i<=0.5+(xmax-xmin)/xtic;i++){
    tmp=pt_w*(i*xtic)/(xmax-xmin);
    x=xmin+i*xtic;
    if (flag) x=fabs(x);
    fprintf(outfile,"%.1f -17 moveto (%g) ctrshow\n",tmp,x);
    fprintf(outfile,"%.1f 0 moveto 0 -5 rlineto stroke\n",tmp);
  }
  fprintf(outfile,"%f -32 moveto\n",0.5*pt_w);
  if (title)
    fprintf(outfile,"(%s) ctrshow\n",title);
}

static void dos_plot_eps(FILE *outfile, double e_min, double e_max,
			 double *bins, int nbins,
			 struct es *e, char *title, int plot_flags){
  double xmin,xmax,ymin,ymax,ytic,xtic,tmp,tmp2,bin_width,scale;
  int i,pt_h,pt_w,ns,itmp,first_bin,last_bin;

  if (flags&AU)
    scale=1/H_eV;
  else
    scale=1;

  if (nbins<2) error_exit("Must have at least two points for a plot");

  bin_width=(e_max-e_min)/(nbins-1);
  
  pt_h=300;
  pt_w=500;
  
  if ((e->nspins>1)&&(plot_flags&PLT_MIRROR)){
    if (!(plot_flags&PLT_ROTATE)) pt_h=500;
    ymin=1e99;
    ymax=-1e99;
    for (ns=0;ns<e->nspins;ns++){
      if (ns&1)
	for(i=0;i<nbins;i++)
	  ymin=min(ymin,-bins[ns*nbins+i]);
      else
	for(i=0;i<nbins;i++)
	  ymax=max(ymax,bins[ns*nbins+i]);
    }
  }
  else{
    if (plot_flags&PLT_ROTATE) pt_w=300;
    ymin=0;
    ymax=-1e99;
    for (i=0;i<e->nspins*nbins;i++)
      ymax=max(ymax,bins[i]);
  }

  ytic=tic(ymax-ymin);
  ymin=ytic*floor(ymin/ytic);
  ymax=ytic*ceil(ymax/ytic);
  xtic=tic(e_max-e_min);
  xmin=xtic*floor(e_min/xtic);
  xmax=xtic*ceil(e_max/xtic);
  
  fprintf(outfile,"%%!PS-Adobe-2.0 EPSF-2.0\n");
  fprintf(outfile,"%%%%BoundingBox: 0 0 %d %d\n",pt_w+100,pt_h+100);
  fprintf(outfile,"%%%%LanguageLevel: 2\n"
	  "%%%%EndComments\n\n"
	  "/ctrshow { dup stringwidth pop -0.5 mul 0 rmoveto show } def\n"
	  "/rshow { dup stringwidth pop neg 0 rmoveto show } def\n\n");


  fprintf(outfile,"50 50 translate\n");
  if (title)
    fprintf(outfile,"/Helvetica-Bold 14 selectfont\n"
	    " %.1f %d moveto (%s) ctrshow\n",0.5*pt_w,pt_h+20,title);
  fprintf(outfile,"/Helvetica 12 selectfont\n");
  fprintf(outfile,"0 0 moveto 0 %d lineto %d %d lineto %d 0 lineto "
	  "closepath stroke\n",pt_h,pt_w,pt_h,pt_w);

  if (plot_flags&PLT_ROTATE){
    if (plot_flags&PLT_PH){
      xaxis_eps(outfile,ymin,ymax,ytic,pt_w,NULL,1);
      fprintf(outfile,"(g(w)) stringwidth pop -0.5 mul 0 rmoveto\n"
	      "(g\\() show /Symbol 12 selectfont (w) show"
	      "/Helvetica 12 selectfont (\\)) show\n");
      yaxis_eps(outfile,xmin,xmax,xtic,pt_h,NULL,0);
      fprintf(outfile,"gsave 90 rotate\n"
	      " (w(cm-1)) stringwidth pop -0.5 mul 0 rmoveto\n"
	      " /Symbol 12 selectfont (w) show /Helvetica 12 selectfont\n"
	      " (\\(cm) show 0 -5 rmoveto (-1) show 0 5 rmoveto (\\)) show\n"
	      "grestore\n");
    }
    else{
      xaxis_eps(outfile,ymin,ymax,ytic,pt_w,"g(E)",1);
      if (flags&AU)
	yaxis_eps(outfile,xmin,xmax,xtic,pt_h,"E(Ha)",0);
      else
	yaxis_eps(outfile,xmin,xmax,xtic,pt_h,"E(eV)",0);
    }
  }
  else{
    if (plot_flags&PLT_PH){
      xaxis_eps(outfile,xmin,xmax,xtic,pt_w,NULL,0);
      fprintf(outfile,"(w(cm-1)) stringwidth pop -0.5 mul 0 rmoveto\n"
	      " /Symbol 12 selectfont (w) show /Helvetica 12 selectfont\n"
	      " (\\(cm) show 0 -5 rmoveto (-1) show 0 5 rmoveto (\\)) show\n");
      yaxis_eps(outfile,ymin,ymax,ytic,pt_h,NULL,1);
      fprintf(outfile,"gsave 90 rotate\n"
	      "(g(w)) stringwidth pop -0.5 mul 0 rmoveto\n"
	      "(g\\() show /Symbol 12 selectfont (w) show"
	      "/Helvetica 12 selectfont (\\)) show\n"
	      "grestore\n");
    }
    else{
      if (flags&AU)
	xaxis_eps(outfile,xmin,xmax,xtic,pt_w,"E(Ha)",0);
      else
	xaxis_eps(outfile,xmin,xmax,xtic,pt_w,"E(eV)",0);
      yaxis_eps(outfile,ymin,ymax,ytic,pt_h,"g(E)",1);
    }
  }

  if (e->e_fermi){
    tmp=(*e->e_fermi*scale-xmin)/(xmax-xmin);
    if (plot_flags&PLT_ROTATE)
      fprintf(outfile,"%d %.1f",pt_w+5,tmp*pt_h-3);
    else
      fprintf(outfile,"%.1f %d",tmp*pt_w+8,pt_h-18);
    fprintf(outfile," moveto (E) show 0 -5 rmoveto (F) show\n");
  }

  if (plot_flags&PLT_ROTATE){
    itmp=pt_w;
    pt_w=pt_h;
    pt_h=itmp;
    fprintf(outfile,"\n%% co-ordinates rotated and inverted beyond this point\n"
	    "gsave -90 rotate %d 0 translate\n",-pt_w);
    fprintf(outfile,"-1 1 scale %d 0 translate\n",-pt_w);
  }

  if (e->e_fermi){
    fprintf(outfile,"%% E Fermi line\n");
    tmp=pt_w*(*e->e_fermi*scale-xmin)/(xmax-xmin);
    fprintf(outfile,"[2 5] 0 setdash\n"
		    " %.1f 0 moveto 0 %d rlineto stroke\n"
	    "[] 0 setdash\n",
	    tmp,pt_h);
  }

  if (ymin<0){
    fprintf(outfile,"%% central axis of mirror plot\n");
    tmp2=-pt_h*ymin/(ymax-ymin);
    fprintf(outfile,"0 %.1f moveto %d 0 rlineto stroke\n",tmp2,pt_w);
    for(i=0;i<(xmax-xmin)/xtic;i++){
      tmp=pt_w*(i*xtic)/(xmax-xmin);
      fprintf(outfile,"%.1f %.1f moveto 0 -5 rlineto stroke\n",tmp,tmp2);
    }
  }

  for(ns=0;ns<e->nspins;ns++){
    first_bin=0;
    while((bins[ns*nbins+first_bin]==0)&&(first_bin<nbins)) first_bin++;
    if (first_bin>0) first_bin--;
    last_bin=nbins-1;
    while((bins[ns*nbins+last_bin]==0)&&(last_bin>0)) last_bin--;
    last_bin++;
    fprintf(outfile,"\n%% spin %d\n",ns);
    if (ns==0)
      fprintf(outfile,"0.8 0 0.8 setrgbcolor\n");
    else if (ns==1)
      fprintf(outfile,"0 0.8 0 setrgbcolor\n");
    else if (ns==2)
      fprintf(outfile,"0.8 0 0 setrgbcolor\n");
    else
      fprintf(outfile,"0 0 0 setrgbcolor\n");
    fprintf(outfile,"0 0 moveto\n");
    if ((ns&1)&&(plot_flags&PLT_MIRROR)){
      fprintf(outfile,"%f %f moveto\n",
	      pt_w*(e_min+first_bin*bin_width-xmin)/(xmax-xmin),
	      pt_h*(-bins[ns*nbins]-ymin)/(ymax-ymin));
      for(i=first_bin;i<last_bin;i++)
	fprintf(outfile,"%f %f lineto\n",
		pt_w*(e_min+i*bin_width-xmin)/(xmax-xmin),
		pt_h*(-bins[ns*nbins+i]-ymin)/(ymax-ymin));
    }
    else{
      fprintf(outfile,"%f %f moveto\n",
	      pt_w*(e_min+first_bin*bin_width-xmin)/(xmax-xmin),
	      pt_h*(bins[ns*nbins]-ymin)/(ymax-ymin));
      for(i=first_bin+1;i<last_bin;i++)
	fprintf(outfile,"%f %f lineto\n",
		pt_w*(e_min+i*bin_width-xmin)/(xmax-xmin),
		pt_h*(bins[ns*nbins+i]-ymin)/(ymax-ymin));
    }
    fprintf(outfile,"stroke\n");
  }
  if (plot_flags&PLT_ROTATE){
    fprintf(outfile,"grestore\n%% co-ordinates restored\n");
  }
  
}

static void dos_plot_raw(FILE *outfile, double e_min, double e_max,
			 double *bins, int nbins,
			 struct es *e, char *title, int plot_flags){
  double bin_width,scale;
  int i,ns;
  char *e_unit;

  if (flags&AU){
    e_unit="Ha";
    scale=1/H_eV;
  }
  else{
    e_unit="eV";
    scale=1;
  }

  if (plot_flags&PLT_PH)
    e_unit="cm-1";

  if (nbins<2) error_exit("Must have at least two points for a plot");

  bin_width=(e_max-e_min)/(nbins-1);
  if (title)
    fprintf(outfile,"# %s\n",title);
  fprintf(outfile,"# DOS energy (%s), density%s\n# %d bins\n",e_unit,
	  (e->nspins>1)?" per spin":"",nbins);
  if (e->e_fermi)
    fprintf(outfile,"# Fermi energy %f %s\n",*e->e_fermi*scale,e_unit);
  for(i=0;i<nbins;i++){
    fprintf(outfile,"%12f",e_min+i*bin_width);
    for(ns=0;ns<e->nspins;ns++)
      fprintf(outfile," %10f",bins[ns*nbins+i]);
    fprintf(outfile,"\n");
  }
}

void gauss_dos(struct kpts *k, struct es *e,
	       double e_min, double e_max, double smear,
	       double *bins, int nbins, int gaps){
  double gsmear,bin_width,x,band_max,band_min,tmp,old,z,bbin_min,bbin_max,wt2;
  double total_weight;
  int i,ioff,ns,ik,nb,j,jw,nbin;

  /* Assume "smearing width" is meant to be FWHM
   *
   * For gaussian exp(-0.5(x/s)**2), s=1, FWHM is 2sqrt(2ln2)=2.35482
   *     gaussian exp(-(x/s)**2), s=1, FWHM is 2sqrt(ln2)=1.66511
   * So for FWHM of one, need s=1/1.66511
   * and smearing width depends linearly on s
   * Note erf integrates exp(-x**2)
   */

  gsmear=smear/1.66511;

  for(i=0;i<e->nspins*nbins;i++) bins[i]=0;
  
  if (nbins<2) error_exit("Must have at least two points for a plot");
  bin_width=(e_max-e_min)/(nbins-1);
  jw=max(1,3*smear/bin_width);

  if (debug)
    fprintf(stderr,"minimum bin energy %f, maximum %f,"
	    " bin width %f (eV)\n",e_min,e_max,bin_width);

  if (gaps==0){
    for(ns=0;ns<e->nspins;ns++){
      for(ik=0;ik<k->n;ik++){
	for(nb=0;nb<e->nbands;nb++){
	  ioff=ik*e->nspins*e->nbands+nb+ns*e->nbands;
	  nbin=floor((e->eval[ioff]-e_min)/bin_width+0.5);
	  x=k->kpts[ik].wt/bin_width;
	  if (smear==0){
	    if ((nbin>=0)&&(nbin<nbins))
	      bins[ns*nbins+nbin]+=x;
	  }
	  else{
	    old=1;
	    for(j=jw;j>-jw;j--){
	      z=(e_min+(j+nbin-0.5)*bin_width)-e->eval[ioff];
	      tmp=erf(z/gsmear);
	      if ((nbin+j<nbins)&&(nbin+j>=0))
		bins[ns*nbins+nbin+j]+=0.5*x*(old-tmp);
	      //	  fprintf(stderr,"%2d %f %f\n",j,z,0.5*(old-tmp));
	      old=tmp;
	    }
	    if ((nbin-jw<nbins)&&(nbin-jw>=0))
		bins[ns*nbins+nbin-jw]+=0.5*x*(old+1);
	  }
	}
      }
    }
  }
  else{
    for(ns=0;ns<e->nspins;ns++){
      for(nb=0;nb<e->nbands;nb++){
	band_max=-1e99;
	band_min=-band_max;
	for(ik=0;ik<k->n;ik++){
	  x=e->eval[ik*e->nspins*e->nbands+nb+ns*e->nbands];
	  band_max=max(band_max,x);
	  band_min=min(band_min,x);
	}
	bbin_min=floor((band_min-e_min)/bin_width+0.5);
        bbin_max=floor((band_max-e_min)/bin_width+0.5);
	for(ik=0;ik<k->n;ik++){
	  ioff=ik*e->nspins*e->nbands+nb+ns*e->nbands;
	  nbin=floor((e->eval[ioff]-e_min)/bin_width+0.5);
	  x=k->kpts[ik].wt/bin_width;
	  /* Check what fraction is in occupied band */
	  wt2=0;
	  old=1;
	  for(j=jw;j>-jw;j--){
	    z=(e_min+(j+nbin-0.5)*bin_width)-e->eval[ioff];
	    tmp=erf(z/gsmear);
	    if ((nbin+j>=bbin_min)&&(nbin+j<=bbin_max))
	      wt2+=0.5*(old-tmp);
	    old=tmp;
	  }
	  if ((nbin-jw>=bbin_min)&&(nbin-jw<=bbin_max))
	    wt2+=0.5*(old+1);
	  wt2=1/wt2;
	  /* Add to histogram, weighting if some bins will be rejected */
	  old=1;
	  for(j=jw;j>-jw;j--){
	    z=(e_min+(j+nbin-0.5)*bin_width)-e->eval[ioff];
	    tmp=erf(z/gsmear);
	    if ((nbin+j>=0)&&(nbin+j<nbins)&&
		(nbin+j>=bbin_min)&&(nbin+j<=bbin_max))
	      bins[ns*nbins+nbin+j]+=0.5*x*(old-tmp)*wt2;
	    old=tmp;
	  }
	  if ((nbin-jw>=0)&&(nbin-jw<nbins)){
	    if ((nbin-jw>=bbin_min)&&(nbin-jw<=bbin_max))
	      bins[ns*nbins+nbin-jw]+=0.5*x*(old+1)*wt2;
	  }
	}
      }
    }
  } /* if gaps */ 

  total_weight=0;
  for(ik=0;ik<k->n;ik++)
    total_weight+=k->kpts[ik].wt;

  if (total_weight!=1.0)
    for(i=0;i<e->nspins*nbins;i++)
      bins[i]/=total_weight;

}

/* Start of Bloechl routines */

/* Sort four doubles into ascending order by network sort */
static void sort4(double x[4]){
  double y[4];

  y[0]=min(x[0],x[2]);
  y[2]=max(x[0],x[2]);

  y[1]=min(x[1],x[3]);
  y[3]=max(x[1],x[3]);

  x[0]=min(y[0],y[1]);
  x[1]=max(y[0],y[1]);

  x[2]=min(y[2],y[3]);
  x[3]=max(y[2],y[3]);

  y[1]=min(x[1],x[2]);
  y[2]=max(x[1],x[2]);

  x[1]=y[1];
  x[2]=y[2];
}

/* Bloechl's tetrahedral interpolation DoS function */
static double bl(double e,double t[4]){
  double e21,e31,e41,e32,e42,e43;
  if (e<=t[0])
    return 0;
  else if (e<=t[1]){
    e21=t[1]-t[0];
    e31=t[2]-t[0];
    e41=t[3]-t[0];
    return 3*(e-t[0])*(e-t[0])/(e21*e31*e41);
  }
  else if (e<=t[2]){
    e21=t[1]-t[0];
    e31=t[2]-t[0];
    e41=t[3]-t[0];
    e32=t[2]-t[1];
    e42=t[3]-t[1];
    return (3*e21+6*(e-t[1])-3*(e31+e42)*(e-t[1])*(e-t[1])/(e32*e42))/(e31*e41);
  }
  else if (e<t[3]){
    e41=t[3]-t[0];
    e42=t[3]-t[1];
    e43=t[3]-t[2];
    return 3*(e-t[3])*(e-t[3])/(e41*e42*e43);
  }
  return 0;
}


int blochl(struct unit_cell *c, struct contents *m,
	    struct kpts *kp, struct es *elect, struct symmetry *rs,
	    double e_min, double e_max, double *y, int n){
  int i,ib,ns,n1,n2,grid[3],ngpts,*mapping,kcube[8],ktet[4],tt;
  int ikx,iky,ikz,ikx1,iky1,ikz1,tet;
  double k_off[3],e[4],shortest,x;
  struct atom a;

  /* The possible sets of tetrhedron vertices for a cube with vertices ordered
   * (0,0,0) (1,0,0) (0,1,0) (1,1,0)
   * (0,0,1) (1,0,1) (0,1,1) (1,1,1)
   */
  int c2t[4][6][4]={{
    {0,1,3,7},
    {0,1,5,7},
    {0,2,3,7},
    {0,2,6,7},
    {0,4,5,7},
    {0,4,6,7}
    },{
    {1,0,2,6},
    {1,0,4,6},
    {1,3,2,6},
    {1,3,7,6},
    {1,5,4,6},
    {1,5,7,6}
    },{               /* Form next two by +/-4 from first two */
    {4,5,7,3},
    {4,5,1,3},
    {4,6,7,3},
    {4,6,2,3},
    {4,0,1,3},
    {4,0,2,3}
    },{
    {5,4,6,2},
    {5,4,0,2},
    {5,7,6,2},
    {5,7,3,2},
    {5,1,0,2},
    {5,1,3,2}
    }};

  for(i=0;i<elect->nspins*n;i++) y[i]=0;

  grid_detect(kp,grid,k_off);

  ngpts=grid[0]*grid[1]*grid[2];

  mapping=malloc(ngpts*sizeof(int));
  if (!mapping) error_exit("Malloc error for mapping");
  if (kgrid_expand(c,m,kp,rs,grid,k_off,mapping)){
    if (debug) fprintf(stderr,"kgrid_expand failed\n");
    free(mapping);
    return 1;
  };

  /* Find shortest diagonal */
  init_atoms(&a,1);
  a.frac[0]=a.frac[1]=a.frac[2]=1;
  addabs(&a,1,c->recip);
  shortest=vmod2(a.abs);
  tt=0;
  a.frac[1]=a.frac[2]=1; a.frac[0]=-1;
  addabs(&a,1,c->recip);
  x=vmod2(a.abs);
  if (x<shortest){
    shortest=x;
    tt=1;
  }
  a.frac[0]=a.frac[1]=1; a.frac[2]=-1;
  addabs(&a,1,c->recip);
  x=vmod2(a.abs);
  if (x<shortest){
    shortest=x;
    tt=2;
  }
  a.frac[0]=a.frac[2]=1; a.frac[1]=-1; 
  addabs(&a,1,c->recip);
  x=vmod2(a.abs);
  if (x<shortest){
    shortest=x;
    tt=3;
  }
  if (debug){
    fprintf(stderr,"Shortest diagonal for Bloechl is ");
    if (tt==0)
      fprintf(stderr," 1  1  1\n");
    else if (tt==1)
      fprintf(stderr,"-1  1  1\n");
    else if (tt==2)
      fprintf(stderr," 1  1 -1\n");
    else if (tt==3)
      fprintf(stderr," 1 -1  1\n");
  }
  
  for(ikx=0;ikx<grid[0];ikx++){
    ikx1=(ikx+1)%grid[0];
    for(iky=0;iky<grid[1];iky++){
      iky1=(iky+1)%grid[1];
      for(ikz=0;ikz<grid[2];ikz++){
	ikz1=(ikz+1)%grid[2];
	kcube[0]=ikz +grid[2]*iky +grid[1]*grid[2]*ikx;
	kcube[1]=ikz +grid[2]*iky +grid[1]*grid[2]*ikx1;
	kcube[2]=ikz +grid[2]*iky1+grid[1]*grid[2]*ikx;
	kcube[3]=ikz +grid[2]*iky1+grid[1]*grid[2]*ikx1;
	kcube[4]=ikz1+grid[2]*iky+ grid[1]*grid[2]*ikx;
	kcube[5]=ikz1+grid[2]*iky+ grid[1]*grid[2]*ikx1;
	kcube[6]=ikz1+grid[2]*iky1+grid[1]*grid[2]*ikx;
	kcube[7]=ikz1+grid[2]*iky1+grid[1]*grid[2]*ikx1;
	for(tet=0;tet<6;tet++){
	  for(i=0;i<4;i++)
	    ktet[i]=mapping[kcube[c2t[tt][tet][i]]];
	  for(ns=0;ns<elect->nspins;ns++){
	    for(ib=0;ib<elect->nbands;ib++){
	      for(i=0;i<4;i++)
		e[i]=elect->eval[ns*elect->nbands+ib+
				 ktet[i]*elect->nbands*elect->nspins];
	      sort4(e);
	      n1=n*(e[0]-e_min)/(e_max-e_min);
	      n2=n*(e[3]-e_min)/(e_max-e_min)+1;
	      n1=max(0,n1);
	      n2=min(n-1,n2);
	      for(i=n1;i<n2;i++)
		y[ns*n+i]+=bl(e_min+i*(e_max-e_min)/(n-1),e);
	    }
	  }
	}
      }
    }
  }

  free(mapping);
  
  for(i=0;i<elect->nspins*n;i++) y[i]/=6*ngpts;
  return 0;
}