%------------------------------------------------------------------------%
%SBSRFOFC(Single Buffer,Single regime,First order,Finite capacity)       %
%------------------------------------------------------------------------%
%Author: H. Emre KANKAYA
%INPUT FORM:    
%SBSRFOFC(c,r,lambda,nu,n,B)
%Ex:SBSRFOFC(1,0.7,1,1,2,3)
%PARAMETERS:
%C          :service rate
%R          :rate per user
%LAMBDA     :on rate per user
%NU         :off rate per user
%N          :number of users
%B          :max of buffer (buffer size)
%
%BACKGROUND PROCESS:
%------------------
%lambda->
% /\ /    \
%O  1  ... N
% \/ \    /
%    <-nu*i
%
%BUFFER:
%-------
%i)single buffer
%ii)first order
%iii)single regime
%iv)finite capacity
%   B        0
%   |--------|_
%ir-->        _ c->  
%   |--------|
%------------------------------------------------------------------------%
function SBSRFOFC=sys1(c,r,lambda,nu,n,B)
%------------------------------------------------------------------------%
ctime=cputime;
%---------------------------CONSTRUCTION---------------------------------%
%A) construct the Q matrix
Q=zeros(n+1,n+1);
Q(1,1)=-n*lambda;
Q(1,2)=n*lambda;
Q(n+1,n+1)=-n*nu;
Q(n+1,n)=n*nu;
if n>1
    for i=2:n
        Q(i,i)=-((n+1-i)*lambda+(i-1)*nu);
        Q(i,i-1)=(i-1)*nu;
        Q(i,i+1)=(n+1-i)*lambda;
    end
end
%B) Construct R matrix
R=zeros(n+1,n+1);
for i=1:n+1
    R(i,i)=-c+(i-1)*r;
end
[v,d]=eig(Q',R');
%C) Reorder the eigenvalues and eigenvectors:
[y,index]=sort(diag(d));
dnew=zeros(n+1,n+1);
vnew=zeros(n+1,n+1);
for i=1:n+1
    dnew(i,i)=d(index(i),index(i));
    vnew(i,:)=(v(:,index(i)))';
end
Q
%------------------------------SOLUTIONS---------------------------------%
%A) Pi(steady-state) distribution (distribution at buffer size B)
Qnew=[Q ones(n+1,1)];
colnew=[zeros(1,n+1), 1];
sol=colnew/Qnew;


noofL=floor(c/r)+1;%Number of states with negative rate
noofU=n+1-noofL;%Number of states with positive rate
%B) BOUNDARY CONDITIONS: 
vu=vnew(:,n+2-noofU:n+1);
vLex1=zeros(n+1);
for i=1:n+1
    vLex1(i,:)=vnew(i,:)*exp(dnew(i,i)*B);
end
vLex2=vLex1(:,1:noofL);
vec2=[ sol(1:noofL) zeros(1,noofU)];
%C)SOLUTION TO THE PROBLEM
sol2=vec2/[vLex2 vu ];
disp('bitti')
ctime-cputime
%-----------------------------PLOTS&OUTPUTS------------------------------%
%A) PLOTS:
b1=floor(B)+1;
noofsamples=10000;
deltay=b1/noofsamples;
Fg=zeros(noofsamples+1,n+1);
no=floor(B/deltay);
for i1=0:no-1
    for j1=1:n+1
        Fg(i1+1,:)=Fg(i1+1,:)+sol2(j1)*vnew(j1,:)*exp(dnew(j1,j1)*i1*deltay);
    end
end
for i=no:noofsamples
        Fg(i,:)=sol;
end

i=0:deltay:b1;
hold on
plot(i,Fg)
str='';
hold on;
for k=1:n+1%3
    str1=strcat('state',num2str(k-1));
    str=strvcat(str,str1);
end%//3//
legend(str,0);

%B) OUTPUTS:
disp('# OF STATES WITH POSITIVE RATE:')
disp(noofL)
disp('# OF STATES WITH NEGATIVE RATE:')
disp(noofU)
disp('COEFFICIENTS:')
disp(sol2)
disp('ACCUMULATIONS AT 0:')
disp(Fg(1,:))
disp('ACCUMULATIONS AT B:')
disp(Fg(no,:)-Fg(no-1,:))
ctime-cputime