/* Plot a histogram using the arrays x and y to represent data       */
/* values and frequencies respectively. If centre is false, x values     */
/* denote the lower edge of the bin, and if centre is true, they        */
/* they denote the centre of the bin                                */

#include "plplot.h"
#include <math.h>

void plbin(nbin,x,y,centre)
int nbin;
int centre;
float x[], y[];
{

      int i;
      float xmin, xmax, vpwxmi, vpwxma, vpwymi, vpwyma;

      int level;
      glev(&level);
      if (level<3)fatal("Please set up window before calling PLBIN.");

/*    Check x[i] are in ascending order */

      for(i=0;i<nbin-1;i++)
        if (x[i] >= x[i+1])
           fatal("Elements of X(*) must be increasing in PLBIN.");

      gvpw(&vpwxmi,&vpwxma,&vpwymi,&vpwyma);
      if (!centre) {
        for(i=0;i<nbin;i++){
          pljoin(x[i],vpwymi,x[i],y[i]);
          pljoin(x[i],y[i],x[i+1],y[i]);
          pljoin(x[i+1],y[i],x[i+1],vpwymi);
        }
        if (x[nbin-1] < vpwxma) {
          pljoin(x[nbin-1],vpwymi,x[nbin-1],y[nbin-1]);
          pljoin(x[nbin-1],y[nbin-1],vpwxma,y[nbin-1]);
          pljoin(vpwxma,y[nbin-1],vpwxma,vpwymi);
        }
      }
      else {
        if (nbin < 2) return;
        xmin = vpwxmi;
        xmax = max(0.5*(x[1]+x[2]),vpwxmi);
        if (xmin < xmax) {
          pljoin(xmin,vpwymi,xmin,y[1]);
          pljoin(xmin,y[1],xmax,y[1]);
          pljoin(xmax,y[1],xmax,vpwymi);
        }
        for (i=1;i<nbin-1;i++) {
          xmin = xmax;
          xmax = min(0.5*(x[i]+x[i+1]),vpwxma);
          pljoin(xmin,vpwymi,xmin,y[i]);
          pljoin(xmin,y[i],xmax,y[i]);
          pljoin(xmax,y[i],xmax,vpwymi);
        }
        xmin = xmax;
        xmax = vpwxma;
        if (xmin < xmax) {
          pljoin(xmin,vpwymi,xmin,y[nbin]);
          pljoin(xmin,y[nbin],xmax,y[nbin]);
          pljoin(xmax,y[nbin],xmax,vpwymi);
        }
      }
}