/* Compute occupancy and blocking for simulated telephone calls */
/* Usage: cc -o calls calls.c -lm                               */
/*        ./calls                                               */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define NUM_CALLS 5
#define LAMBDA 1.0
#define MEAN_CALL_DURATION 3.0
#define MAX_CAPACITY 5
#define END_OF_TIME 99999.9

#define DEBUG 1
#define OCCUPANCY 1

/***********************************************************************/
/*                 RANDOM NUMBER GENERATION STUFF                      */
/***********************************************************************/

/* Parameters for random number generation. */
#define MAX_INT 2147483647       /* Maximum positive integer 2^31 - 1 */

/* Generate a random floating point value uniformly distributed in [0,1] */
float Uniform01()
  {
    float randnum;

    /* get a random positive integer from random() */
    randnum = (float) 1.0 * random();

    /* divide by max int to get something in 0..1  */
    randnum = (float) randnum / (1.0 * MAX_INT); 

    return( randnum );
  }

/* Generate a random floating point number from an exponential    */
/* distribution with mean mu.                                     */

float Exponential(mu)
    float mu;
  {
    float randnum, ans;

    randnum = Uniform01();
    ans = -(mu) * log(randnum);

    return( ans );
  }

/***********************************************************************/
/*                 MAIN PROGRAM                                        */
/***********************************************************************/
extern void qsort();

int mycompare(num1, num2)
     float *num1, *num2;
{
  if( *num1 < *num2 )
    return( -1);
  else if( *num1 > *num2 )
    return( 1);
  else return(0);
}

int main()
  {

    float starts[NUM_CALLS+1], ends[NUM_CALLS+1];
    float durations[NUM_CALLS];
    float mytime, lasttime, occtime[NUM_CALLS];
    int i, j, num, calls, lost, wasblocked, startindex, endindex;
    float mean, var;

    /* Seed the PRNG */
    srandom(time(NULL));
    srandom(1234567);

    /* Make huge arrays with all the timing information */
    for( i = 0; i < NUM_CALLS; i++ )
      {
	if( i == 0 )
	  starts[i] = Exponential(1.0/LAMBDA);
	else starts[i] = starts[i-1] + Exponential(1.0/LAMBDA);
	durations[i] = Exponential(MEAN_CALL_DURATION);
	ends[i] = starts[i] + durations[i];
#ifdef DEBUG
	fprintf(stderr, "Caller %d start %8.6f duration %8.6f end %8.6f\n",
	       i, starts[i], durations[i], ends[i]);
#endif
      }
    starts[NUM_CALLS] = END_OF_TIME;
    ends[NUM_CALLS] = END_OF_TIME;

    /* sort the ending times into timestamp order */
    qsort(ends, NUM_CALLS, sizeof(float), mycompare);
#ifdef DEBUG2
    for( i = 0; i < NUM_CALLS; i++ )
      fprintf(stderr, "Start %8.6f  end %8.6f\n", starts[i], ends[i]);
#endif

    /* Initialization */
    for( i = 0; i < NUM_CALLS; i++ )
	occtime[i] = 0.0;

    startindex = 0;
    endindex = 0;
    calls = 0;
    lost = 0;
    lasttime = 0.0;

    /* print a play-by-play report of call occupancy */
    printf("0.0 0\n");
    for( i = 0; i < NUM_CALLS; )
      {
	while( (startindex < NUM_CALLS) && (starts[startindex] < ends[endindex]) )
	  {
	    mytime = starts[startindex];
	    occtime[calls] += mytime - lasttime;
	    lasttime = mytime;
	    printf("%f %d\n", mytime, calls);
	    calls++;
	    if( calls > MAX_CAPACITY )
	      {
		lost++;
		calls--;
		fprintf(stderr, "Blocked call %d at time %8.6f\n", startindex, mytime);
		starts[startindex] *= -1;
	      }
	    printf("%f %d\n", mytime, calls);
	    startindex++;
	  }
	while( (ends[endindex] <= starts[startindex]) && (i < NUM_CALLS) )
	  {
	    mytime = ends[endindex];
	    occtime[calls] += mytime - lasttime;
	    lasttime = mytime;
	    wasblocked = 0;
	    for( j = 0; j < startindex; j++ )
	      {
		if( (starts[j] < 0) && (-starts[j]+durations[j]==mytime) )
		  {
		    wasblocked++;
		    fprintf(stderr, "Ignoring end time %8.6f for blocked call %d\n",
			   mytime, j);
		  }
	      }
	    printf("%f %d\n", mytime, calls);
	    calls--;
	    if( wasblocked )
	      calls++;
	    printf("%f %d\n", mytime, calls);
	    endindex++;
	    i++;
	  }
      }
    fprintf(stderr, "Call Simulation: CALLS %d  LAMBDA %3.1f  DUR %3.1f  C %d\n",
	    NUM_CALLS, LAMBDA, MEAN_CALL_DURATION, MAX_CAPACITY);

    fprintf(stderr, "Loss probability: %8.6f\n", 1.0*lost/NUM_CALLS);

#ifdef OCCUPANCY
    fprintf(stderr, "Occupancy distribution at time %8.6f\n", mytime);
    mean = 0.0;
    var = 0.0;
    for( i = 0; i <= MAX_CAPACITY; i++ )
      {
	fprintf(stderr, "%d %8.6f\n", i, occtime[i]/mytime);
	mean += i * occtime[i]/mytime;
	var += i * i * occtime[i]/mytime;
      }
    var -= mean * mean;
    fprintf(stderr, "Mean occupancy: %8.6f\n", mean);
    fprintf(stderr, "Variance of occupancy: %8.6f\n", var);
#endif
    return(0);
  }