/* ===================================================================
     Sample program to calculate density and sum of states of H2O
  =================================================================== */

#include <stdio.h>
#define MAXVIV 20
#define MAXGRN 3000

int directCount(int nv, double *freq, double sizg, double emax,
 double *dens, double *sums, int *jfrq);

int main() {
  int nVib = 3, nGrain, roundFreq[MAXVIV], igr;
  double energyMax = 60000., sizeGrain = 100.,
         frequency[MAXVIV] = {3657., 1595., 3756.},
         rho[MAXGRN], W[MAXGRN];
  nGrain = directCount(nVib, frequency, sizeGrain, energyMax, rho, W,
                       roundFreq);
  if (nGrain < 0) {
    printf(" Error termination of directCount with error code = %d\n",
           nGrain);
    return 1;
  }
  printf("   E[cm-1],     rho[cm],           W\n");
  for (igr = 0; igr < nGrain; igr++) {
    printf("%10.2f,%12.4E,%12.4E\n", igr * sizeGrain, rho[igr], W[igr]);
  }
  return 0;
}

/* ---------------------------------------------------------------------
     Direct count (Beyer-Swinehart algorithm)
  ----- input -----
   nv     : number of vibrators
   freq[] : vibrational frequencies [cm-1]
   sizg   : grain size [cm-1]
   emax   : maximum energy [cm-1]
  ----- output ----
   dens[] : density of states (subscript 0 corresponds to E=0)
   sums[] : sum of states
   jfrq[] : frequencies in 'grain size' unit (rounded to integer)
   (function value) : > 0 ... number of grains
                      < 0 ... error status flag
                       (-1: error in freq, -2: error in emax)
--------------------------------------------------------------------- */

int directCount(int nv, double *freq, double sizg, double emax,
 double *dens, double *sums, int *jfrq) {
  int ng, iv, ig, jfr, jg;
  /* --- Set frequencies and max energy in 'grain size' unit */
  if (nv <= 0) { return -1; }
  for (iv = 0; iv < nv; iv++) {
    jfrq[iv] = (int)(freq[iv] / sizg + 0.5);
    if (jfrq[iv] == 0) { return -1; }
  }
  ng = (int)(emax / sizg + 1.5);
  if (ng <= 1) { return -2; }
  /* --- Reset arrays */
  for (ig = 0; ig < ng; ig++) {
    dens[ig] = 0;
    sums[ig] = 1;
  }
  dens[0] = 1;
  /* --- Beyer-Swinehart direct count */
  for (iv = 0; iv < nv; iv++) {
    jfr = jfrq[iv];
    for (ig = jfr; ig < ng; ig++) {
      jg = ig - jfr;
      dens[ig] += dens[jg];
      sums[ig] += sums[jg];
    }
  }
  /* --- Convert unit of density of states from grain-1 to (cm-1)-1 */
  for (ig = 0; ig < ng; ig++) {
    dens[ig] /= sizg;
  }
  return ng;
}