// count_pulses.c - utility function for use with the Ablation_Station_XYZ
//                  LabView program that controls the diamond milling process.
//                  It receives a time-sampled sequence from the energy meter
//                  and counts the number of pulses in the sequence. It 
//                  also computes the average pulse height and returns
//                  the value in one of its arguments.
//
// author: pratt at phys.uconn.edu, richard.t.jones at uconn.edu
// version: february 9, 2015

// uncomment this line to create an executable test module
// for these functions, reads data from a file called pulses.dat
#define TESTING 1

//#include <Windows.h>
//#include <extcode.h>
#include <stdlib.h>
#include <stdio.h>

#define PRESAMPLE_WIDTH 100
#define ENERGY_THRESHOLD 0.01
#define MINIMUM_PULSE_WIDTH 100
#define MAXIMUM_PULSE_WIDTH 150

int //__declspec(dllexport) 
count_pulses(double *pulse_data,
             double sample_size, 
             double *average_pulse_height)
{
    int pulse_count = 0;
    double pulse_height_sum = 0;
    double pulse_height_max = 0;
    int pulse_width = 0;
    int trigger_armed;

    double presample_buffer[PRESAMPLE_WIDTH];
    double presample_sum;
    int i, ipre;
    for (ipre = 0; ipre < PRESAMPLE_WIDTH; ++ipre) {
        presample_buffer[ipre] = 999;
    }
    presample_sum = 999 * PRESAMPLE_WIDTH;

    trigger_armed = 1;
    for (i = 0; i < sample_size; ++i) {
        double energy = pulse_data[i] - presample_sum / PRESAMPLE_WIDTH;
        if (trigger_armed) {
            if (energy > ENERGY_THRESHOLD) {
                trigger_armed = 0;
                pulse_height_max = energy;
                pulse_width = 1;
            }
            else {
                presample_sum += pulse_data[i];
                presample_sum -= presample_buffer[ipre % PRESAMPLE_WIDTH];
                presample_buffer[ipre % PRESAMPLE_WIDTH] = pulse_data[i];
                ++ipre;
            }
        }
        else if (pulse_width > MINIMUM_PULSE_WIDTH && energy < ENERGY_THRESHOLD
                 || pulse_width > MAXIMUM_PULSE_WIDTH)
        {
            pulse_height_sum += pulse_height_max;
            ++pulse_count;
            trigger_armed = 1;
        }
        else {
            if (energy > pulse_height_max) {
                pulse_height_max = energy;
            }
            ++pulse_width;
        }
    }

    if (average_pulse_height) {
        if (pulse_count > 0) {
            *average_pulse_height = pulse_height_sum / pulse_count;
        }
        else {
            *average_pulse_height = 0;
        }
    }
    return pulse_count;
}

#if TESTING

int main()
{
   FILE *fp = fopen("pulses.dat","r");
   if (fp == 0) {
      fprintf(stderr, "unable to open input file pulses.dat");
      exit(1);
   }
   else {
      int pulse_count;
      double average_pulse_height;
      double databuf[1000000];
      int ndata = 0;
      while (fscanf(fp, "%lf", &databuf[ndata]) != EOF) {
         ++ndata;
      }
      printf("read %d values from input file\n", ndata);
      pulse_count = count_pulses(databuf, ndata, &average_pulse_height);
      printf("%d pulses found, average pulse height %lf\n",
             pulse_count, average_pulse_height);
      fclose(fp);
   }
}

#endif