function [max_delay, arr_seed, dep_seed] = lunch_truck(arr_seed,dep_seed)
%Programmer: Raghu Pasupathy
%Last Modified: Fall 2007

global TNOW TMAX PEAK_START WARMUP EVENT_LIST;
    %disp(sprintf('Generating output data...'));
    initialize(); %initialize event list
    ARRIVAL_TIMES =[]; DEPARTURE_TIMES = [];    
    TNOW = EVENT_LIST(1,1);     %advance simulation clock
    nextevent_type = EVENT_LIST(1,2);   
    while ( TNOW <= TMAX )
        switch nextevent_type
            case 1 %arrival event
                ARRIVAL_TIMES = [ ARRIVAL_TIMES; TNOW]; %update arrival time array
                [arr_seed, dep_seed] = arrival(arr_seed, dep_seed);
            case 2 %departure event
                DEPARTURE_TIMES = [ DEPARTURE_TIMES; TNOW]; %update departure time array
                [dep_seed] = service_end(dep_seed);       
        end
        TNOW = EVENT_LIST(1,1);             %advance simulation clock
        nextevent_type = EVENT_LIST(1,2);   %set the next event type
    end
    M = [ARRIVAL_TIMES(1:size(DEPARTURE_TIMES,1)), DEPARTURE_TIMES, DEPARTURE_TIMES-ARRIVAL_TIMES(1:size(DEPARTURE_TIMES,1))];
    i = 1;
    while ( M(i,1) < WARMUP )
        i = i + 1;
    end
    dlmwrite('dog.dat',M(i:size(M,1),1:3));
    %disp(sprintf('Done.'));
    
    % accumulate required data
    max_delay = 0; sum = 0; sum2 = 0;
    for i = 1:size(DEPARTURE_TIMES,1)
        sum = sum + DEPARTURE_TIMES(i) - ARRIVAL_TIMES(i);
        sum2 = sum2 + (DEPARTURE_TIMES(i) - ARRIVAL_TIMES(i))^2;
        if ( ARRIVAL_TIMES(i) > PEAK_START )
            if ( DEPARTURE_TIMES(i) - ARRIVAL_TIMES(i) > max_delay )
                max_delay = DEPARTURE_TIMES(i) - ARRIVAL_TIMES(i);
            end
        end
    end
    avg_delay = sum / size(DEPARTURE_TIMES,1);
    se = sqrt( ( sum2 / size(DEPARTURE_TIMES,1) - avg_delay^2 ) / size(DEPARTURE_TIMES,1) );
            
%--------------------------------------------------------------
function initialize()
global lambda mu TNOW TMAX PEAK_START WARMUP EVENT_LIST NO_SYSTEM; 
    lambda = 1;           %set arrival rate
    mu = 0.8;               %set mean service time
    TNOW = 0;               %initialize simulation time
    TMAX = 300;             %set maximum simulation time
    PEAK_START = 180;
    WARMUP = 0;
    EVENT_LIST = [ 0, 1 ];  %initialize event list
    NO_SYSTEM = 0;
%--------------------------------------------------------------
function [arr_seed, dep_seed] = arrival(arr_seed, dep_seed)
global lambda mu TNOW NO_SYSTEM;
    delete();           % delete first row of event list
    [ arr_seed, u ] = u16807d( arr_seed );
    trnsfer1 = TNOW + ( -1 / lambda ) * log( 1 - u );   %generate next arrival time
    trnsfer2 = 1;       % event type is arrival
    trnsfer = [trnsfer1, trnsfer2];   
    insert( trnsfer ); %append the next arrival record to the event list
    NO_SYSTEM = NO_SYSTEM + 1;
    if ( NO_SYSTEM == 1 )           %if server is empty, schedule a departure
        [ dep_seed, u ] = u16807d( dep_seed );
        trnsfer1 = TNOW + ( -mu ) * log( 1 - u );   % generate the departure time
        trnsfer2 = 2;           %event type is departure
        trnsfer = [trnsfer1, trnsfer2];
        insert( trnsfer );      %append the next arrival record to the event list
    end
%---------------------------------------------------------------
function [dep_seed] = service_end(dep_seed)
global mu TNOW NO_SYSTEM;
    delete();           % delete first row of event list
    NO_SYSTEM = NO_SYSTEM - 1;
    if ( NO_SYSTEM >= 1 )           %if server is not empty, schedule a departure
        [ dep_seed, u ] = u16807d( dep_seed );
        trnsfer1 = TNOW + ( -mu ) * log( 1 - u );   % generate the departure time
        trnsfer2 = 2;           %event type is departure
        trnsfer = [trnsfer1, trnsfer2];
        insert( trnsfer );  %update event list
    end
%---------------------------------------------------------------
function delete()
global EVENT_LIST;
% delete first record in the event list
    if ( size(EVENT_LIST,1) >= 2 )
        EVENT_LIST = EVENT_LIST(2:size(EVENT_LIST,1), 1:size(EVENT_LIST,2)); % delete first row
    else
        EVENT_LIST = [];    
    end
%-----------------------------------------------------------------
function insert(trnsfer)
global EVENT_LIST;
    EVENT_LIST = [ EVENT_LIST; trnsfer ];   % append the next record
    [Y I] = sort(EVENT_LIST(1:size(EVENT_LIST,1),1));  %sort the first column (event time) in the event list
    Z = EVENT_LIST(1:size(EVENT_LIST,1),2);  %get the second column (event type) in the event list
    Z = Z(I);               %order the second column correctly
    EVENT_LIST = [Y,Z];    % form the sorted event list
%----------------------------------------------------------------------
function [iseed,u16807d]=u16807d( iseed )
%..........................................................................
%     bruce schmeiser     january 1989.                                   .
%     a linear congruential pseudorandom number generator                 .
%       using constant 16807 and modulus (2**31)-1.                       .
%              iseed = iseed*16807 (mod 2^31 -1)                          .
%     for implementations that don't require double precision, see        .
%       s.k. park and k.w. miller, "random numbers generators: good       .
%         ones are hard to find," cacm, 31, 10 (October 1988), 1192-1201. .
%     in correct implementations, starting with a seed value of 1 will    .
%       result in a seed value of 1043618065 on call number 10001.        .
%..........................................................................
%     input:  iseed.   integer.                                           . 
%                        chosen from [1,2147483646] on the first call.    .
%                        thereafter, the value returned from the last call.
%     output: iseed.   integer.                                           .
%                        to be used in the next call.                     .
%     output: u16807d.  real.                                             .
%                        a pseudorandom number in (0,1).                  .
%..........................................................................
%iseed=0;     
u16807d=0;
while (u16807d<=0 | u16807d>=1)
    iseed = mod (iseed * 16807,2147483647);
    u16807d = iseed / 2147483648;
end
%--------------------------------------------------------------------------