%Модель M/M/1/n с обслуживанием только парами
%Для графического интерфейса пришлось использовать num2str, что 
%сильно ухудшило скорость вычислений.

1;
global n;  %максимальное число требований в системе (больше пяти)
global strla; %интенсивность обслуживания. функция от t
global strmu; %интенсивность поступления. функция от t
global filename; %имя сохраняемых графиков
global begt; %начальное время для построения графика на периоде
global oP; % список состояний которые строятся на графике

function r=f1(x,t)
global n;
global strla;
global strmu;
%Здесь можно la и mu записать как функции.
%Если медленно считает, здесь можно записать функции явно
% и будет нормально
la=str2num(strrep (strla,"t",num2str(t,16))); 
mu=str2num(strrep (strmu,"t",num2str(t,16))); 
%mu=4-cos(2*pi*t);
%la=2+sin(2*pi*t);

r(1)=-la*x(1)+mu*x(3);
r(2)=la*x(1)-la*x(2)+mu*x(4);
for k=3:(n-2)
r(k)=la*x(k-1)-(la+mu)*x(k)+mu*x(k+2);
endfor;
r(n-1)=la*x(n-2)-(la+mu)*x(n-1);
r(n)=la*x(n-1)-mu*x(n);
endfunction

function r=initf(nn,muu,laa,ind,begtt,outP)
global n;
global strmu;
global strla;
global filename; 
global begt;
global oP;
n=str2num(get(nn,"string"))+1;
strmu=get(muu,"string"); 
strla=get(laa,"string"); 
filename=get(ind,"string"); 
begt=str2num(get(begtt,"string")); 
oP=get(outP,"string");

t=linspace(0,begt,1000);
tt=linspace(begt,begt+1,1000);
x0=zeros(n,1);
x0(1)=1;
y1=lsode("f1",x0,t);
y11=lsode("f1",y1(end,:),tt);
x0=zeros(n,1);
x0(end)=1;
y2=lsode("f1",x0,t);
y22=lsode("f1",y2(end,:),tt);


tmp='';
for (k=str2num(oP))
tmp=[tmp;';p',num2str(k),'|X(0)=0;'];
endfor

tmp2='';
for (k=str2num(oP))
tmp2=[tmp2;';p',num2str(k),'|X(0)=',int2str(n-1),';'];
endfor

hf = figure ();
hold on;
 plot(t,y1(:,str2num(oP)+1),tmp,t,y2(:,str2num(oP)+1),tmp2)
  xlabel ("t");
ylabel ("probability");
grid minor;
 axis("tight");
hold off;
print ([filename,"p.pdf"]);


hf = figure ();
hold on;
 plot(tt,y11(:,str2num(oP)+1),tmp,tt,y22(:,str2num(oP)+1),tmp2)
  xlabel ("t");
ylabel ("probability");
grid minor;
 axis("tight");
hold off;
print ([filename,"pS.pdf"]);


hf = figure ();
hold on;
plot(t,y1*(0:(n-1))',"k;X(0)=0;",t,y2*(0:(n-1))',["--k;X(0)=",int2str(n-1),";"])
xlabel ("t");
ylabel ("mean");
grid minor;
 axis("tight");
hold off;
print ([filename,"mean.pdf"]);


hf = figure ();
hold on;
plot(tt,y11*(0:(n-1))',"k;X(0)=0;",tt,y22*(0:(n-1))',["--k;X(0)=",int2str(n-1),";"])
xlabel ("t");
ylabel ("mean");
grid minor;
 axis("tight");
hold off;
print ([filename,"meanS.pdf"]);

endfunction

f = figure();
hold on;
b1 = uicontrol (f, "style", "pushbutton","string", "Solve", "position",[10 10 150 30],"callback","initf(nn,muu,laa,ind,begtt,outP)");
t1 = uicontrol (f, "style", "text", "string", "t", "position",[10 40 60 40]);
begtt = uicontrol (f, "style", "edit", "string", "10", "position",[70 40 300 40]);
t1 = uicontrol (f, "style", "text", "string", "n", "position",[10 120 60 40]);
nn = uicontrol (f, "style", "edit", "string", "10", "position",[70 120 300 40]);
t1 = uicontrol (f, "style", "text", "string", "mu", "position",[10 160 60 40]);
muu = uicontrol (f, "style", "edit", "string", "4-cos(2*pi*t)", "position",[70 160 300 40]);
t1 = uicontrol (f, "style", "text", "string", "lambda", "position",[10 200 60 40]);
laa = uicontrol (f, "style", "edit", "string", "2+sin(2*pi*t)", "position",[70 200 300 40]);
t1 = uicontrol (f, "style", "text", "string", "i", "position",[10 240 60 40]);
outP = uicontrol (f, "style", "edit", "string", "[0,1,2]", "position",[70 240 300 40]);
t1 = uicontrol (f, "style", "text", "string", "filename", "position",[10 280 60 40]);
ind = uicontrol (f, "style", "edit", "string", "BDP", "position",[70 280 300 40]);
hold off;
