function [newQN] = addStage(qn, inout, S, rates, P)
% ADDSTAGE add a stage to the queuing network to represent a random environment
% qn: queuing network to be modified
% inout: rates for arriving (first row) and exiting (second row) the current 
%        stage from any existing stage. This is a matrix 2x qn.E.
% S: number of servers in the new stage
% rates: server rates in the new stage
% P: routing probabilities in the new stage
% newQN: modified queuing network

	if ~isa(qn, 'CMCQNCSRE')
		qn = CMCQNCSRE(qn.M, qn.K, 1, qn.N, 0, {qn.S}, {qn.rates}, qn.sched, ...
			{qn.P}, qn.NK, qn.classMatch, qn.refNodes, {'noReset'}, qn.nodeNames, ...
			qn.classNames);
	end
	
	% calculate new E
	newE = qn.E + 1;
	
	% calculate new ENV
	newENV = qn.ENV;
	newENV(1:qn.E, newE) = inout(1,:)'; 
	newENV(newE, 1:qn.E) = inout(2,:);
	for e=1:newE
		newENV(e,e) = newENV(e,e)-sum(newENV(e,:));
	end
	
	% calculate new S
	newS = qn.S;
	newS{newE} = S;
	
	% calculate new rates
	newRates = qn.rates;
	newRates{newE} = rates;
	
	% calculate new P
	newP = qn.P;
	newP{newE} = P;

	% calculate new QN
	newQN = CMCQNCSRE(qn.M, qn.K, newE, qn.N, newENV, newS, newRates, qn.sched, ...
			newP, qn.NK, qn.classMatch, qn.refNodes);
	newQN.nodeNames = qn.nodeNames;
	newQN.classNames = qn.classNames;
end