//========================================================= file = h2d1.c =====
//=  A simple "straight C" H2/D/1 queue simulation                            =
//=============================================================================
//=  Notes: 1) This program is adapted from Figure 1.6 in Simulating          =
//=            Computer Systems, Techniqyes and Tools by M. H. MacDougall     =
//=            (1987).                                                        =
//=         2) The values of SIM_TIME, ARR_TIME, COV, and SERV_TIME need to   =
//=            be set.                                                        =
//=---------------------------------------------------------------------------=
//=  Build: gcc h2d1.c -lm, bcc32 h2d1.c, cl h2d1.c                           =
//=---------------------------------------------------------------------------=
//=  Execute: h2d1                                                            =
//=---------------------------------------------------------------------------=
//=  History: KJC (05/27/99) - Genesis                                        =
//=============================================================================

//----- Include files ---------------------------------------------------------
#include <stdio.h>             // Needed for printf()
#include <stdlib.h>            // Needed for exit() and rand()
#include <math.h>              // Needed for log()

//----- Constants -------------------------------------------------------------
#define SIM_TIME   1.0e6       // Simulation time
#define ARR_TIME   1.25        // Mean time between arrivals
#define COV        2.00        // Coefficient of variation for arrivals
#define SERV_TIME  1.00        // Mean service time

//----- Function prototypes ---------------------------------------------------
double hyper(double x, double cov);  // Generate hyperexponential RV

//===== Main program ==========================================================
void main(void)
{
  double   end_time = SIM_TIME; // Total time to simulate
  double   Ta = ARR_TIME;       // Mean time between arrivals
  double   cov = COV;           // CoV of arrivals
  double   Ts = SERV_TIME;      // Mean service time

  double   time = 0.0;         // Simulation time
  double   t1 = 0.0;           // Time for next event #1 (arrival)
  double   t2 = SIM_TIME;      // Time for next event #2 (departure)
  unsigned long int n = 0;     // Number of customers in the system

  unsigned long int c = 0;     // Number of service completions
  double   b = 0.0;            // Total busy time
  double   s = 0.0;            // Area of number of customers in system
  double   tn = time;          // Variable for "last event time"
  double   tb;                 // Variable for "last start of busy time"
  double   x;                  // Throughput
  double   u;                  // Utilization
  double   l;                  // Mean number in the system
  double   w;                  // Mean residence time

  // Main simulation loop
  while (time < end_time)
  {
    if (t1 < t2)                 // *** Event #1 (arrival)
    {
      time = t1;
      s = s + n * (time - tn);   // Update area under "s" curve
      n++;
      tn = time;                 // tn = "last event time" for next event
      t1 = time + hyper(Ta, cov);
      if (n == 1)
      {
        tb = time;               // Set "last start of busy time"
        t2 = time + Ts;
      }
    }
    else                         // *** Event #2 (departure)
    {
      time = t2;
      s = s + n * (time - tn);   // Update area under "s" curve
      n--;
      tn = time;                 // tn = "last event time" for next event
      c++;                       // Increment number of completions
      if (n > 0)
        t2 = time + Ts;
      else
      {
        t2 = SIM_TIME;
        b = b + time - tb;      // Update busy time sum if empty
      }
    }
  }

  x = c / time;   // Compute throughput rate
  u = b / time;   // Compute server utilization
  l = s / time;   // Compute mean number in system
  w = l / x;      // Compute mean residence or system time

  // Output results
  printf("=============================================================== \n");
  printf("=           *** Results from H2/D/1 Simulation ***            = \n");
  printf("=============================================================== \n");
  printf("=  Total simulated time         = %3.4f sec  \n", end_time);
  printf("=============================================================== \n");
  printf("=  INPUTS: \n");
  printf("=    Mean time between arrivals = %f sec      \n", Ta);
  printf("=    CoV for arrivals           = %f          \n", cov);
  printf("=    Mean service time          = %f sec      \n", Ts);
  printf("=============================================================== \n");
  printf("=  OUTPUTS: \n");
  printf("=    Number of completions      = %ld cust    \n", c);
  printf("=    Throughput rate            = %f cust/sec \n", x);
  printf("=    Server utilization         = %f %%       \n", 100.0 * u);
  printf("=    Mean number in system      = %f cust     \n", l);
  printf("=    Mean residence time        = %f sec      \n", w);
  printf("=============================================================== \n");
}

//===========================================================================
//=  Function to generate hyperexponentially distributed random variables   =
//=    - Input:  Mean value of distribution and coefficient of variation    =
//=    - Output: Returns with hyperexponentially distributed rv             =
//=-------------------------------------------------------------------------=
//=  Uses Morse's method taken from SMPL from Simulating Computer Systems   =
//=  Systems, Techniques and Tools by M. H. MacDougall (1987)               =
//===========================================================================
double hyper(double x, double cov)
{
  double p;                     // Probability value for Morse's method
  double z1, z2;                // Uniform random numbers from 0 to 1
  double hyp_value;             // Computed exponential value to be returned
  double temp;                  // Temporary double value

  // Pull a uniform random number (0 < z1 < 1)
  do
  {
    z1 = ((double) rand() / RAND_MAX);
  }
  while ((z1 == 0) || (z1 == 1));

  // Pull another uniform random number (0 < z2 < 1)
  do
  {
    z2 = ((double) rand() / RAND_MAX);
  }
  while ((z2 == 0) || (z2 == 1));

  // Compute hyperexponential random variable using Morse's method
  temp = cov*cov;
  p = 0.5 * (1.0 - sqrt((temp - 1.0) / (temp + 1.0)));
  if (z1 > p)
    temp = x / (1.0 - p);
  else
    temp = x / p;
  hyp_value = -0.5 * temp * log(z2);

  return(hyp_value);
}