%-----------------------------------------------|
%--SINGLE REGIME - FIRST ORDER - FINITE BUFFER--|
%-----------------------------------------------|
%Author: H. Emre KANKAYA
ctime=cputime;
%-----------------PARAMETERS--------------------|
B=1;%Buffer Length                             %|
Q=[-3,1,1,1;1,-3,1,1;1,1,-2,0;1,1,1,-3];       %|
rates=[1,-1,11,2];                             %|
noofsamples=1000;%For Plotting                 %|
%-----------------------------------------------|


%-------------------PREPARE---------------------|
n=length(rates);%number of states              %|
                                               %|
%A)Reorder Q and R:                            %|
%------------------                            %|
[Y,I]=sort(rates);                             %|
Qnew=zeros(n);                                 %|
for i=1:n                                      %| 
    for j=1:n                                  %| 
        Qnew(i,j)=Q(I(i),I(j));                %|
    end                                        %|
end                                            %|
Q=Qnew;                                        %|
R=diag(Y);                                     %|
%-   -   -   -   -   -   -   -   -   -   -   - %|
%B)Find Eigenvalues and Eigenvectors:          %|
%------------------------------------          %|
[v,d]=eig(Q',R');                              %|
%# of 0 eigenvales=1                           %|
%# of - eigenvalues=noofL                      %|
%# of + eigenvalues=noofU-1                    %|                                                   
%-   -   -   -   -   -   -   -   -   -   -   - %|
%C)Reorder the eigenvalues and eigenvectors:   %|
%-------------------------------------------   %|
[y,index]=sort(diag(d));                       %|
dnew=zeros(n,n);%ordered eigenvalues           %|
vnew=zeros(n,n);%ordered eigenvectors          %|
for i=1:n                                      %|
    dnew(i,i)=d(index(i),index(i));            %|
    vnew(i,:)=(v(:,index(i)))';                %|
end       
dnew
vnew%|
%-----------------------------------------------|



%-------------------SOLUTION--------------------|
%A) Pi(steady-state) Distribution:             %|
%---------------------------------             %|
Qnew=[Q ones(n,1)];                            %|
colnew=[zeros(1,n), 1];                        %|
sol=colnew/Qnew;%steady state distribution     %|
%-   -   -   -   -   -   -   -   -   -   -   - %|
%B) Find # of + and - States:                  %|
%----------------------------                  %|
noofL=0;%# of states with - rate               %|
for i=1:n                                      %|
    if rates(i)<0                              %|
        noofL=noofL+1;                         %|
    end                                        %|
end                                            %|
noofU=n-noofL;%# of states with + rate         %|
%-   -   -   -   -   -   -   -   -   -   -   - %|
%C) Boundary Conditions:                       %|
%-----------------------                       %|
%Boundary Conditions at zero,                  %|
%for states with + rate                        %|
vu=vnew(1:n,noofL+1:noofL+noofU);              %|
%Boundary Conditions at B                      %|
%for states with - rate                        %|
vLex1=zeros(n);                                %|
for i=1:n                                      %|
    vLex1(i,:)=vnew(i,1:n)*exp(dnew(i,i)*B);   %|
end                                            %|
vLex2=vLex1(:,1:noofL);                        %|
%-   -   -   -   -   -   -   -   -   -   -   - %|
%D)Solution:                                   %|
%-----------                                   %|
totalmatrix=[vLex2 vu];                        %|
totalvec=[ sol(1:noofL) zeros(1,noofU)];       %|
sol2=totalvec/totalmatrix;                     %|
%-----------------------------------------------|


%-------------------PLOTS-----------------------|
no=noofsamples-2;                              %|
deltay=B/no;                                   %|
Fg=zeros(noofsamples+1,n);                     %|
for i=1:no+1                                   %|
    for j=1:n                                  %|
        Fg(i,:)=Fg(i,:)+sol2(j)*vnew(j,:)*exp(dnew(j,j)*(i-1)*deltay);
    end                                        %|
end                                            %|
for i=no+1:noofsamples+1                       %|
        Fg(i,:)=sol;                           %|
end                                            %|
i=0:deltay:noofsamples*deltay;                 %|
hold on;                                       %|
plot(i,Fg)                                     %|
str='';                                        %|
hold on;                                       %|
for k=1:n                                      %|
    str1=strcat('state',num2str(I(k)));        %|
    str=strvcat(str,str1);                     %|
end                                            %|
legend(str,0);                                 %|
%-----------------------------------------------|

%-------------------OUTPUTS---------------------|    
disp('--------------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+1,:)-Fg(no,:))                      %|
disp('SYSTEM IS SOLVED IN(sec.)')              %|
disp(-ctime+cputime)                            %|
disp('---------------------------------------')%|
%-----------------------------------------------|