%MVA for C queues,, closed, created from mvaproc, Jan 14 2010;
%param is class 1 demands queues 1 to K, class 1 think, class 2...,
%...class C
%  *****THIS CODE IS FOR TWO QUEUES ONLY*******
%processor-sharing discipline allows different service demands by class
%selectpar corresponds to param, with 1 or selected for sensitivity, 0 not
%column vectors
%K queues, C classes. 
%Npop =[N1,...NC] are integer populations of classes
%Code must be written to correspond to C

function [h]=modelCalc(param, controlled)

global K
global C  %corresponds to the provided special code for the model
Npop(1) = controlled(1);  %unpack the controlled parameters for use here
Npop(2) = controlled(2);
%
%T(i,c,NN1,NN2) is calculated delay at station i for class c, with pop
%NN1-1,NN2-1
%n(i,c,NN1,NN2) is mean customers
%U( ) is utilization
%f(c,NN1,NN2) is throughput
%h returns class 1 delays by station, numbers by station, utilizations by
%station, throughput, then class 2...
%if (C ~= 2) print('bad value of C', C); exit; end
%check C corresponds to assumption of code
%unpack the parameters, C groups
D(:,1) = param(1:K);
Z(1) = param(K+1);
D(:,2) = param(K+2:2*K+1);
Z(2) = param(2*K+2);
%set the array sizes large enough
T = zeros(K,C,Npop(1)+1, Npop(2)+1);
n = zeros(K,C,Npop(1)+1, Npop(2)+1);
U = zeros(K,C,Npop(1)+1, Npop(2)+1);
f = zeros(C,Npop(1)+1, Npop(2)+1);
%initialize the boundaries        
T(:,1,1,1) = 0;
n(:,1,1,1) = 0;
U(:,1,1,1) = 0;
T(:,2,1,1) = 0;
n(:,2,1,1) = 0;
U(:,2,1,1) = 0;
f(1,1,1) = 0;
f(2,1,1) = 0;


%fill boundary NN2=1, set class 2 to zero
for NN1 = 2:Npop(1)+1
    NN2=1;
        NN(1) = NN1 - 1;  %actual class 1 population
        NN(2) = NN2-1;  %never use   
        for i = 1:K %queues
            nsum(i,1) = 0; 
            
            for d = 1:C %competition over all classes
                nsum(i,1) = nsum(i,1) + n(i,d,NN1-1,NN2);
                
            end
        end
        %or, nsum(:,1) = sum(n(:,1,NN1-1,NN2)) etc
        for c = 1:1
            sumT(c) = Z(c);
            %processor-sharing discipline
            for i =1:K
                T(i,c,NN1,NN2) = nsum(i,c)*D(i,c) + D(i,c);
                sumT(c) = sumT(c) + T(i,c,NN1,NN2);
            end
            f(c,NN1,NN2)= NN(c)/sumT(c);
            for i = 1:K
                n(i,c,NN1,NN2) = f(c,NN1,NN2)*T(i,c,NN1,NN2);
                U(i,c,NN1,NN2) = f(c,NN1,NN2) * D(i,c);
            end
        end
    T(:,2,NN1,1) = 0;
f(2,NN1,1) = 0;
n(:,2,NN1,1) = 0;
U(:,2,NN1,1) = 0;
end



%fill boundary for NN1=1, set class 1 = 0
NN1 = 1;
    for NN2 = 2:Npop(2)+1
        NN(1) = NN1 - 1;  %actual class 1 population, never use
        NN(2) = NN2-1;     
        for i = 1:K            
            nsum(i,2) = 0;
            for d = 1:C                
                nsum(i,2) = nsum(i,2) + n(i,d,NN1,NN2-1);
            end
        end
        %or, nsum(:,1) = sum(n(:,1,NN1-1,NN2)) etc
        for c = 2:2
            sumT(c) = Z(c);
            %processor-sharing discipline
            for i =1:K
                T(i,c,NN1,NN2) = nsum(i,c)*D(i,c) + D(i,c);
                sumT(c) = sumT(c) + T(i,c,NN1,NN2);
            end
            f(c,NN1,NN2)= NN(c)/sumT(c);
            for i = 1:K
                n(i,c,NN1,NN2) = f(c,NN1,NN2)*T(i,c,NN1,NN2);
                U(i,c,NN1,NN2) = f(c,NN1,NN2) * D(i,c);
            end
        end
        T(:,1,1,NN2) = 0;
        f(1,1,NN2) = 0;
        n(:,1,1,NN2) = 0;
        U(:,1,1,NN2) = 0;
    end




%fill interior both NN1 and NN1 > 1
for NN1 = 2:Npop(1)+1
    for NN2 = 2:Npop(2)+1
        NN(1) = NN1 - 1;  %actual class 1 population
        NN(2) = NN2-1;     
        for i = 1:K
            nsum(i,1) = 0; 
            nsum(i,2) = 0;
            for d = 1:C
                nsum(i,1) = nsum(i,1) + n(i,d,NN1-1,NN2);
                nsum(i,2) = nsum(i,2) + n(i,d,NN1,NN2-1);
            end
        end
        %or, nsum(:,1) = sum(n(:,1,NN1-1,NN2)) etc
        for c = 1:C
            sumT(c) = Z(c);
            %processor-sharing discipline
            for i =1:K
                T(i,c,NN1,NN2) = nsum(i,c)*D(i,c) + D(i,c);
                sumT(c) = sumT(c) + T(i,c,NN1,NN2);
            end
            f(c,NN1,NN2)= NN(c)/sumT(c);
            for i = 1:K
                n(i,c,NN1,NN2) = f(c,NN1,NN2)*T(i,c,NN1,NN2);
                U(i,c,NN1,NN2) = f(c,NN1,NN2) * D(i,c);
            end
        end
    end
end

for c = 1:C
    for j=1:K        
        h((c-1)*(3*K+1)+j)= T(j,c,Npop(1)+1,Npop(2)+1);
        h((c-1)*(3*K+1)+K+j)= n(j,c,Npop(1)+1,Npop(2)+1);
        h((c-1)*(3*K+1)+2*K+j)= U(j,c,Npop(1)+1,Npop(2)+1);
    end
    h((c-1)*(3*K+1)+3*K+1)= f(c,Npop(1)+1,Npop(2)+1);
end
h=h';%retyrns a column vector
%T,n,U,f,h



