function st = simulate_buffer(max_time)
   st = initial_state();
   %Sts = [st];
   %while !stop(Sts(end), max_time)
   while !stop(st, max_time)
	   st = update(st);
	   %Sts = [Sts, st];
   end
end

function st = initial_state()
   st.tm_stamp = 0;
   st.srv_time = +inf;              % server is idle
   st.srv_last = 0;                 % last served time
   st.enq_time = +inf;              % buffer is free
   st.next_arr = next_arrival();
   st.n_accept = 0;
   st.n_reject = 0;  
   st.usd_time = 0; 
   st.wtq_time = 0; 
   st.queue_sz = 0;
   st.queue_av = 0;
   printf("starting at %6.2f; next arrival at %6.2f\n", 
         0, st.next_arr);
end

function st1 = update(st0)
   if       st0.srv_time == inf  && ...             % server is empty    
            st0.enq_time < inf            % buffer is not empty
      st1 = dequeue_event(st0);    % a  dequeue event  
   elseif   st0.srv_time < st0.next_arr   % a service ends before arrival  
      st1 = serviced_event(st0);   % a serviced event  
   elseif   st0.enq_time < inf && ...                  % queue is full
            st0.srv_time < inf && ...                   % buffer is full
            st0.srv_time > st0.next_arr   % the service ends after arrival
      st1 = reject_event(st0);     % a reject event
   elseif   st0.enq_time == inf           % queue is not full
       st1 = enqueue_event(st0);   % an enqueue event 
   else
      printf("error - unclear transition in state:\n")
      st1 = st0
      st1.time_stamp = + inf;
   end;
end

function st1 = reject_event(st0)
   st1 = st0;
   t0 = st0.tm_stamp;
   t1 = st0.next_arr;
   st1.tm_stamp = t1;
   st1.next_arr = t1 + next_arrival();    % compute next arrival time
   st1.n_reject = st0.n_reject + 1;  	   % increase rejected requests
   st1.queue_av = st0.queue_av + (t1-t0)*st0.queue_sz;
   printf("rejected at %6.2f: busy till %6.2f; enq at %6.2f; arrive \
at %6.2f\n", t1, st1.srv_time, st1.enq_time, st1.next_arr);
end 

function st1 = serviced_event(st0)
   st1 = st0;
   t0 = st0.tm_stamp;
   t1 = st0.srv_time;
   st1.tm_stamp = t1;
   st1.srv_time = +inf;     	   % server becomes idle
   st1.srv_last = t1;             % update last served 
   st1.queue_av = st0.queue_av + (t1-t0)*st0.queue_sz;
   printf("serviced at %6.2f: busy till %6.2f; enq at %6.2f; arrive \
at %6.2f\n", t1, st1.srv_time, st1.enq_time, st1.next_arr);
end
       
function st1 = enqueue_event(st0)
   st1 = st0;
   t0 = st0.tm_stamp;
   t1 = st0.next_arr;
   st1.tm_stamp = t1;
   st1.enq_time = t1;                      % update enqueued time
	st1.next_arr = t1 + next_arrival();    % compute next arrival time
   st1.n_accept = st0.n_accept + 1;  	  % increase rejected requests
   st1.queue_sz = st0.queue_sz + 1;
   st1.queue_av = st0.queue_av + (t1-t0)*st0.queue_sz;
   printf("enqueued at %6.2f: busy till %6.2f; enq at %6.2f; arrive \
at %6.2f\n", t1, st1.srv_time, st1.enq_time, st1.next_arr);
end
       
function st1 = dequeue_event(st0)
   st1 = st0;
   t0 = st0.tm_stamp;
   t1 = max(st0.enq_time, st0.srv_last);
   st1.tm_stamp = t1;
   st1.enq_time = +inf;                 % free the buffer
   st1.srv_time = t1 + service_time();   % update service time   
   st1.usd_time = st1.usd_time + (st1.srv_time - t1);
   st1.wtq_time = st0.wtq_time + (t1-st0.enq_time);
   st1.queue_sz = st0.queue_sz - 1;
   st1.queue_av = st0.queue_av + (t1-t0)*st0.queue_sz;
   printf("dequeued at %6.2f: busy till %6.2f; enq at %6.2f; arrive \
at %6.2f\n", t1, st1.srv_time, st1.enq_time, st1.next_arr);
end

function t_next = next_arrival()
% an exponential distribution with lambda = 0.2
% i.e. mean time between arrivals of 1/lambda = 5	
   lambda = 0.2; r = rand();	
   t_next = -log(1-r)/lambda;
end

function t = erlang()
% função de distribuição
   M = 0.15;
   T = 25;
   lambda = 0.4;
   k = 2;
   aceite = false;
   while ! aceite
      t = T * rand();
      r = M * rand();
      aceite = (r < (lambda^2 * t * exp(-lambda * t)));
   end
end

function t = service_time()
   mode = 1;
   if mode == 1  % an erlang (k=2, lambda= 0.4) distribution 
      t = erlang();        % average time = 2 / 0.4  = 5
   else           % a uniform distribution between 2 and 8
       t = 2 + 6*rand(); % average time = (2+8)/2 = 5
   end
end

function b = stop(st, max_time)
   b = (st.tm_stamp >= max_time);
end


