function [] = queue_simulator()

R_bar = []

for trial = 1:5


% Parameters
N = 10000;
lambda = 1;
s = .75;
mu = 1 / s;

% Calculate true value of R_bar
U = lambda * s;
R_true = s / (1-U);

% Generate exponential interarrival and service times using inverse CDF
r = rand(N,1);
interarrival_times = -log(r) ./ lambda;

r = rand(N,1);
service_times = -log(r) ./ mu;

% Use interarrival times to generate actual arrival time of each customer
arrival_times = cumsum(interarrival_times);

% Create vectors to track the time each customer enters service and departs
enter_service_times = zeros(N,1);
completion_times = zeros(N,1);

% Set the enter_service_time and completion_time for the first customer
enter_service_times(1) = arrival_times(1);
completion_times(1) = enter_service_times(1) + service_times(1);

% Process all of the other customers
for i = 2:N
    
    % A customer enters service either when it arrives or when the
    % previous customer completes, whichever is later
    enter_service_times(i) = max(arrival_times(i), completion_times(i-1));
    
    % Completion time 
    completion_times(i) = enter_service_times(i) + service_times(i);
end

% Waiting time is the time between arriving and entering service
waiting_times = enter_service_times - arrival_times;

% Residence time is the time between arriving and finally departing
residence_times = completion_times - arrival_times;

% Estimated average waiting and residence times
W_bar = sum(waiting_times) / length(waiting_times);
R_bar(end+1) = sum(residence_times) / length(residence_times);


end


R_bar