/* deadtime.c - Written by R.A. Davis December2, 1996 */

/* Purpose: To calculate deadtime fractions for the case of */
/*          random arrival time and fixed processing time */
/*          Follows the calculation of G.P. Heath in */
/*          NIM A278 (1989) 431. */

/* Usage:   Three operating modes: */

/* deadtime 1 N x ....to calculate deadtime for a single N,x pair */

/* deadtime 2 N xlo xhi xstep ....to calculate for fixed N over x range */

/* deadtime 3 x Nlo Nhi Nstep ....to calculate for fixed x over N range */

/* The last two modes produce a file deadtime.out in the format */
/* N,x,deadtime for mode 2 and */
/* x,N,deadtime for mode 3    */
/* This file may be read into paw by the command: */
/* vec/read a,b,c deadtime.out */



#include 
#include 
#include 
#include 
#include 
#include 
#include 

main(int argc, char *argv[]) {

  double x,xlo,xhi,xstep;
  int n,nlo,nhi,nstep,Status,OutFileDes;
  int nbytes,outlen;
  extern double dead(int n, double x);
  char OutputFile[] = "./deadtime.out";
  char output[20];
  char *format1 = "%i %f %f \n";
  char *format2 = "%f %i %f \n";

  printf("Starting deadtime calculation program.\n");

  /* Read mode of operation */

  if (atoi(argv[1]) == 1 ) {  /* Calculate single deadtime */

    printf("Single deadtime calculation mode.\n");

    if( (argc < 4) || (argc > 4) ) { /* Check for correct # of arguments */
      usage(1);
      exit(1);
    }
 
    n = atoi(argv[2]);
    x = strtod(argv[3],NULL);
    printf("x => %f\n",x);
    printf("N => %i\n",n);

    printf("Dead time fraction is %f\n",dead(n,x));

  }

  else if (atoi(argv[1]) == 2 ) { /* Hold N and vary x */

    printf("Holding N constant and varying x.\n");

    if( (argc < 6) || (argc > 6) ) { /* Check for correct # of arguments */
      usage(2);
      exit(1);
    }
 
    n = atoi(argv[2]);
    xlo = strtod(argv[3],NULL);
    xhi = strtod(argv[4],NULL);
    xstep = strtod(argv[5],NULL);
    printf("N => %i\n",n);
    printf("xlo => %f\n",xlo);
    printf("xhi => %f\n",xhi);
    printf("xstep => %f\n",xstep);

    /* Open output file */

    Status = creat(OutputFile,0666);
    if (Status < 0) {
      printf("Problem creating output file - exiting.\n");
      exit(3);
    }    

    OutFileDes = open(OutputFile,O_WRONLY);
    if (OutFileDes < 0) {
      printf("Problem opening output file - exiting.\n");
      exit(4);
    }    

    x = xlo;

    while ( x <= xhi ) {

      sprintf(output,format1,n,x,dead(n,x));
      printf("N, x, deadtime fraction ==> %s\n",output);
      outlen = strlen(output);
      nbytes = write(OutFileDes,output,outlen);
      if ( nbytes < 0 ) {
        printf("Problem writing to output file. Exiting.\n");
        exit(5);
      }
      x = x + xstep;

    }

   /* Close output file */

    close(OutFileDes);
  }

  else if (atoi(argv[1]) == 3 ) { /* Hold x and vary N */

    printf("Holding x constant and varying N.\n");

    if( (argc < 6) || (argc > 6) ) { /* Check for correct # of arguments */
      usage(3);
      exit(1);
    }
 
    x = strtod(argv[2],NULL);
    nlo = atoi(argv[3]);
    nhi = atoi(argv[4]);
    nstep = atoi(argv[5]);
    printf("x => %f\n",x);
    printf("nlo => %i\n",nlo);
    printf("nhi => %i\n",nhi);
    printf("nstep => %i\n",nstep);

    /* Open output file */

    Status = creat(OutputFile,0666);
    if (Status < 0) {
      printf("Problem creating output file - exiting.\n");
      exit(3);
    }    

    OutFileDes = open(OutputFile,O_WRONLY);
    if (OutFileDes < 0) {
      printf("Problem opening output file - exiting.\n");
      exit(4);
    }    

    n = nlo;

    while ( n <= nhi ) {

      sprintf(output,format2,x,n,dead(n,x));
      printf("x, N, deadtime fraction ==> %s\n",output);
      outlen = strlen(output);
      nbytes = write(OutFileDes,output,outlen);
      if ( nbytes < 0 ) {
        printf("Problem writing to output file. Exiting.\n");
        exit(5);
      }
      n = n + nstep;

    }

   /* Close output file */

    close(OutFileDes);
  }

  else {
    usage(0);
  }

  return(0);

}

double dead(int n,double x) {

  double S,dead;
  int j;
  extern double jfac(int j);

  /* Put in n=0 term by hand */

  S = 0;
  if ( n == 0) {
    S = 1;
  }
  j = 0;

  while (j <= n-1) {
 
    S = S + exp((double)n*x) * ((pow((double)(n - j), (double)j)) * pow((-x * exp(-x)),(double)j)) / jfac(j);

    j++;

      }

  dead = 1.0 - S / (1.0 + x * S);

  return(dead);
}

double jfac(j) {

  int i;
  double retval;

  i = 1;
  retval = 1.0;
  while (i <= j) {
    retval = retval * (double)i;
    i++;
  }
  return(retval);
}

int usage(int flag) {

  if (flag == 1) {
    printf("Single deadtime usage: deadtime 1 N x\n");
  }  

  else if (flag == 2) {
    printf("Fix N vary x usage: deadtime 2 N xlo xhi xstep\n");
  }

  else if (flag == 3) {
    printf("Fix x vary N usage: deadtime 3 x Nlo Nhi Nstep\n");
  }

  else {
    printf("May be used in one of three modes:\n");
    printf("  deadtime 1   Calculate deadtime for single N,x pair.\n");
    printf("  deadtime 2   Calculate deadtime over range of x for fixed N.\n");
    printf("  deadtime 3   Calculate deadtime over range of N for fixed x.\n");
  }

  return(0);

}