function [avgvisits] = avgVisits(qn)
% AVGVISITS calculate the average number of visits of each server (row) for each
% class (column). This function assumed no class switching.
% qn: queuing network
% avgvisits: number of average visits calculated by this function

	persistent envCTMC solCTMC;

	if isa(qn,'CMCQNCSRE')
		
		if ~exist('envCTMC', 'var') || ~isequal(envCTMC, qn.ENV)
			envCTMC = qn.ENV;
			solCTMC = ctmc_solve(qn.ENV);
		end
		
		P = qn.P;
		E = qn.E;
		w = solCTMC;
		S = [qn.S{1}];
	else
		P{1} = qn.P;
		E = 1;
		w = 1;
		S = [qn.S];
	end
	
	avgvisits = zeros(qn.M, qn.K);

%	classMatch = qn.classMatch;
	
	for e=1:E
%		THIS CODE DOES NOT SUPPORT CLASS SWITCHING
		for ii=1:qn.K
			p = P{e}(ii:qn.K:end, ii:qn.K:end);
			p1 = p(1,2:qn.M);
			p2 = p(2:qn.M,2:qn.M);
			avgvisits(1,ii) = avgvisits(1,ii) + w(e);
			avgvisits(2:end,ii) = avgvisits(2:end,ii) + w(e).*(p1/(eye(qn.M-1)-p2))';
			avgvisits(2:end,ii) = avgvisits(2:end,ii)./S(2:end);
		end

%     THIS COMMENTED CODE SUPPORTS CLASS SWITCHING, BUT IS NOT MAINTAINED
% 		for ii=1:qn.K
% 			
% 			eventualVisitProb = reshape( myCQN.P((refNode-1)*myCQN.K+idxOrigClasses(k),idxNR)*eventualVisit, myCQN.K , myCQN.M-sum(delayRefNodes>0) )'; %probability of eventual visit 
% 			eventualVisitProb = eventualVisitProb(:,classMatch(k,:)==1);
% 			avgvisits(:,ii) = eventualVisitProb;
% 		end
		
		
	end
end