# copyright Barry L. Nelson and Michael R. Taaffe 2000

with(linalg):

PhPhInfMCNMoments := proc(A, lambda, R, B, mu, K, Pr, nodelist, lower, upper)

#
# this procedure returns a plot of the mean and variance of
# the number in system an at a list of nodes
# for a [Ph_t/Ph_t/Infinity]^K network with R different arrival processes 
# 
# initial conditions for the queue are empty and idle
#
# this version is designed to work with Maple V release 5
#              
# last update 7/28/00
#
# Inputs:
# A = Markov chain transition matrix for the arrival process
#     [vector of matrices of dimension R]
# lambda = arrival rate vector
#     [vector of vectors of dimension R]
# R = number of distinct Ph arrival processes
# B = Markov chain transition matrix for the service process
#     [vector of matrices of dimension K]
# mu = service rate vector
#     [vector of vectors of dimension K]
# Pr = Markov routing probabilities
# nodelist = list of nodes to be plotted
# lower, upper = range of plots

	local r, deq, initial, temp, mplots;
	if R = 1 then
	   RETURN(PhPhInfNetMoments(A, lambda, B, mu, K, Pr, nodelist, lower, upper));
	else
	   mplots := vector(R);
	   temp := PhPhInfNetMoments(A[1], lambda[1], B, mu, K, Pr[1], nodelist, lower, upper);
	   mplots[1] := temp;
	   for r from 2 to R do
	      temp := PhPhInfNetMoments(A[r], lambda[r], B, mu, K, Pr[r], nodelist, lower, upper);
	      mplots[r] := temp;
	   od;
	   RETURN(mplots);
	fi;
	end;