% Feb 2010, adapted from mvaproc_think
% modelCalc for one class by MVA
% code for analytic derivatives is commented out, to work with modelSens
% K = queues is global
% param(i) = demand for queue i, i = 1 to K
%param(K+1) = think time Z


% %Example Data input for an example in regressModel:
% nPotential = 4;
% parNamesCell = {'  D1  ', 'D2  ','  D3  ','  Z  '};
% parInteger = [0,0,0,0]';
% parInit = [1,1,1,10];
% parFit = [1,1,1,0];
% pControlInd = [1];
% TMAX = 10;
% mpotential = 11;
% dataNamesCell = {'  N  ', '  T1  ', '  T2  ','  T3  ',' n1  ',' n2 ',' n3 ',' U1 ',' U2 ',' U3 ',' f '};             
% 
% data = [1 1.21384 4.73153 2.64049 0.060692 0.2365765 0.1320245 0.0594781 0.264966 0.16847 0.05
% 2 1.60529 2.32118 4.72039 0.17337132 0.25068744 0.50980212 0.14541 0.204289 0.473156 0.108
% 3 1.71797 4.4539 4.26163 0.25425956 0.6591772 0.63072124 0.210235 0.432874 0.538983 0.148
% 4 2.82288 5.48302 6.825 0.4093176 0.7950379 0.989625 0.303433 0.586554 0.658857 0.145
% 5 3.49714 6.47212 8.42731 0.60850236 1.12614888 1.46635194 0.372242 0.499684 0.788018 0.174
% 6 3.00036 4.3094 13.6851 0.6150738 0.883427 2.8054455 0.389882 0.523607 0.875007 0.205
% 7 3.0938 5.55339 11.3586 0.7486996 1.34392038 2.7487812 0.439166 0.599545 0.870314 0.242
% 8 2.31285 7.90299 19.5202 0.4440672 1.51737408 3.7478784 0.354667 0.575443 0.925672 0.192
% 9 3.25927 5.07819 15.1546 0.85392874 1.33048578 3.9705052 0.53194 0.703916 0.990093 0.262
% 10 3.28088 7.04839 19.8235 0.81693912 1.75504911 4.9360515 0.437181 0.76528 0.974517 0.249 ];
% 
% dataVariances = [0.0109,0.0374,0.0745,0.0000737,0.0001531,0.0000872,.0787,.0304,.0792,0.0000154];
% maxIter = 20;
% convtest = 0.00001;
% %Example data input



function h=modelCalc(param,controlled)
%N customers, K queues
global K;
N = controlled(1);
Z= param(K+1);  %think time
D = zeros(1,K);   %row vector
for i = 1:K
    D(i)=param(i);
end
T(1,:)=D;
%dTdD(1,:,:)=eye(K);
f(1)=1.0/(sum(D)+Z);
%dfdD(1,:)=-f(1)^2;
Mean(1,:)=f(1)*D;
U(1,:) = f(1)*D;

%dMeandD(1,:,:)=f(1)^2*D'*ones(1,K) + f(1)*eye(K);
% i=1;   %initial values for derivs of f and U
% for k = 1:K
%       dfdD(i,k) = -(f(i)^2/i)*sum(dTdD(i,:,k));
%       for j = 1:K
%          dUdD(i,j,k) = dfdD(i,k) * D(j) + f(i)*(k==j);
%       end;
%  	end;
    
for i = 2:N
   denom = Z;
   for j = 1:K   %queue j
      T(i,j)= (1 + Mean(i-1,j))*D(j);
      denom = denom + T(i,j);
%       for k = 1:K   %wrto D(k)
%          dTdD(i,j,k)= (1+Mean(i-1,j))*(j==k) + D(j)*dMeandD(i-1,j,k);
%       end;
   end;
   
   f(i)=i/denom;
   U(i,:) = f(i)* D;
%    for k = 1:K   %wrto D(k)
%       dfdD(i,k) = -(f(i)^2/i)*sum(dTdD(i,:,k));
%       for j = 1:K
%          dUdD(i,j,k) = dfdD(i,k) * D(j) + f(i)*(k==j);
%       end;
%  	end;
   
   Mean(i,:) = f(i)*T(i,:);
%    for k = 1:K
%       dMeandD(i,:,k) = dfdD(i,k)*T(i,:) + f(i)*dTdD(i,:,k);
%    end
end;
for j = 1:K
   h(j) = T(N,j);
   h(K+j)= Mean(N,j);
   h(2*K+j)= U(N,j);
%    for k=1:K
%       H(j,k) = dTdD(N,j,k);
%       H(j,K+k) = dMeandD(N,j,k);
%       H(j,2*K+k) = dUdD(N,j,k);
%    end
%   H(j,3*K+1)=dfdD(N,j);
end
h(3*K+1) = f(N);
%H = H(1:K,:);
h = h';
%H = H';






