% simulation program          two class: exhaustive/gated constant switch negative expont.message length %
% r pseudostations and r/2 stations in all
clear;
r=20;s=1*ones(1,r);s_load=0.6; heavy=0.9; ratio=0.5;% heavy:station asymm. ratio:two class asymm.  
for i=1:r/2,
     if i==1,
         l(2*i-1)=s_load/s(1)*heavy*ratio;    l(2*i)=s_load/s(1)*heavy*(1-ratio); 
     else
         l(2*i-1)=s_load/s(1)*(1-heavy)/(r/2-1)*ratio;    l(2*i)=s_load/s(1)*(1-heavy)/(r/2-1)*(1-ratio);
     end
     rule(2*i-1)=1;          rule(2*i)=0;  % 1=exhaustive,0:gated, service discipline
     u(2*i-1)=0;             u(2*i)=0.2*s(1);
end

load=zeros(1,r);s_u=0;x=0;s_load=0;max_l=0;% calculate the load
for i=1:r,
    load(i)=l(i)*s(i);s_u=s_u+u(i);s_load=s_load+load(i);x=x+load(i)+l(i)*u(i); %stability condition<1
    if max_l<l(i),max_l=l(i);end
end
c=s_u/(1-s_load);% mean polling period

t_sim=1e7;%simulation time

arr=zeros(r,l_wlist_p);%saving arrival instant of messages
%intermediate variables
sum_w=zeros(1,r);sum_l=zeros(1,r);sum_s=zeros(1,r);s_sum_w=zeros(1,r);xx=zeros(1,r);t2=zeros(1,r);
waits=zeros(1,r);%the number of waiting  messages
n_srvd=zeros(1,r);%served numbers
%simulation process
t=0;T=0;
for i=1:r,% the first message
   arr(i,1)=(-1.0)/l(i)*log(1-rand);
end

while t<t_sim, 
   T=T+1;% new polling period, begin to poll from station 1
   for i=1:r, %   begin serve a pseudostation     
      % % managing the messages comes in a determining times
      if arr(i,1)<t,
        waits(i)=1;%緩沖中的第一個數據產生時間是在此論旬點之前產生的
        tx=arr(i,1)+(-1.0)/l(i)*log(1-rand);
        arr(i,waits(i)+1)=tx;
        while tx<t,         %表明指針所指的下一個信元是在當前時刻之前到達的，應當入等待隊列;     
           waits(i)=waits(i)+1;% a new waiting message
           % if waits(i)>l_wlist_p  aa='fafads'
           %  end
           tx=arr(i,waits(i))+(-1.0)/l(i)*log(1-rand);
           arr(i,waits(i)+1)=tx; %arrival instant
        end
         sum_l(i)=sum_l(i)+waits(i);  % the sum of waiting length 
     end
     if waits(i)>0  % there is message waiting
         for j=1:waits(i), 
             n_srvd(i)=n_srvd(i)+1;%the number of served messages
             sum_w(i)=sum_w(i)+t-arr(i,j);%
             if sum_w(i)>100,s_sum_w(i)=s_sum_w(i)+sum_w(i);sum_w(i)=0;end %the sum of waiting times
             t=t-s(i)*log(1-rand);%service time
         end 
         arr(i,1)=arr(i,waits(i)+1);
         waits(i)=0;% complete of service of message come before the gate 
         if rule(i)==1, %  exhaustive service
             while arr(i,1)<t,
                 n_srvd(i)=n_srvd(i)+1;% the number of served messages
                 sum_w(i)=sum_w(i)+t-arr(i,1);%
                 if sum_w(i)>100,s_sum_w(i)=s_sum_w(i)+sum_w(i);sum_w(i)=0;end % the sum of waiting times
                 t=t-s(i)*log(1-rand);;%service time
                 arr(i,1)=arr(i,1)+(-1.0)/l(i)*log(1-rand);
             end
         end
     end %   if l_waitlist(i)>0  %
     t=t+u(i);%-1.0)*u(i)*log(1-rand);   %constant switch
 end %  serve  the next pseudostation 
end   % the simulation time complete

%statistical works
for i=1:r,
   w_avrg(i)=(s_sum_w(i)+sum_w(i))/n_srvd(i); %mean waiting times
   l_avrg(i)=sum_l(i)/T;% mean waiting lenght when the service begin
end
w_avrg(1:2:r) %the first pseudostation
w_avrg(2:2:r) %the second pseudostation