// C++/CSIM Model of M/M/1 queue using Markov Chain approach
//                track the state probability for P[0] to P[4]  using PASTA

#include "cpp.h"			// class definitions

#define NARS 50000
#define IAR_TM 2.0
#define SRV_TM 1.0


double NT;		// records the cumulative area under the curve N(t)
double t, next_t;       // represents current time (t) and next event time (next_t)

double arrival_time;    // holding time for arrival event
double departure_time;  // holding time for departure event
double N;		// records the average # of customer, NT/t for large t
long   state, next_state;   // state of the Markov Chain
int    next_event;      // arrival = 1, departure = 2

double P[6];		// P[5] = sum of Prob[5] to Prob[\infty]
double sum;
double total_arrival;   // track total number of arrivals

extern "C" void sim(int, char **);

void sim(int argc, char *argv[])
{
	set_model_name("M/M/1 MC Approach");
	create("sim");

	// initialization
        NT = 0.0;
        t = 0.0;
	next_t = 0.0;
        state = 0;
        P[0]=P[1]=P[2]=P[3]=P[4]=P[5]=0.0;
        total_arrival = 0.0;


	// generate events
       t = clock;
       for (int i=0; i<NARS; i++){
	 if (state == 0){ // current state = 0
             hold(expntl(IAR_TM));
             next_t = clock;
	     NT += ((double) state) * (next_t - t);

             // update state prob
             P[0] += 1.0;
             total_arrival += 1.0;

             t = next_t;      // update current time
             state += 1;    // update state for arrival event

         }
         else { // current state != 0
             arrival_time = expntl(IAR_TM);
             departure_time = expntl(SRV_TM);
             if (arrival_time < departure_time) 
                     next_event = 1;
             else if (arrival_time > departure_time) 
                     next_event = 2;
             else {// randomly pick arrival or departure
                 if ((uniform(0.0,1.0)) < 0.5)
                      next_event = 1;
                 else next_event = 2;
             }
             
            if (next_event == 1)
               hold(arrival_time);
            else
               hold(departure_time);

             next_t = clock;
             NT += ((double) state) * (next_t - t);

             // update state prob
             if ( (next_event == 1) && (state <=4) ){ 
               P[state] += 1.0;
               total_arrival += 1.0;
             }
             else if (next_event == 1){
               P[5] += 1.0;
               total_arrival += 1.0;
             }

             t = next_t;      // update current time
             if (next_event == 1)
                state += 1;    // update state for arrival event
             else
                state -= 1;
         }
   
       } // end for loop

       // average number of customer
       N = NT/clock;


       printf("Average number of customer = %lf\n", N);
       printf("Average Response Time = %lf\n", N/(1.0/IAR_TM)); // Little's Result

       sum = 0.0;
       for (int j=0; j<=5;j++){
        printf("Prob[%2d] = %lf\n", j, P[j]/total_arrival);
        sum += P[j]/total_arrival;
       }
       printf("Sum all probabilities = %lf\n", sum);

}