function [newQN] = splitStation(qn, x, proportions, cpu)
% SPLITSTATION This function splits a station that contains replicas.
% qn: queuing network
% x: index of the station that will be split
% proportions: vector of numbers which indicate the proportions in which
%					the station should be divided. The sum of the vector is 1. Its
%				   cardinality is the number of new stations that will be added to
%				   the QN minus one.
% cpu: vector of number which indicates the number of cpus of the given
% proportion
% newQN: new queuing network calculated by this function

	resNum = 1:length(proportions);

	% remove zero proportions
	sel = proportions == 0;
	proportions(sel) = [];
	resNum(sel) = [];
	cpu(sel) = [];

	lproportions = length(proportions);
	
	% calculate new M
	newM = qn.M + lproportions - 1;
	
	% calculate new S
	if isa(qn,'CMCQNCSRE')
		newS = cell(qn.E, 1);
		for e=1:qn.E
			newS{e} = [qn.S{e}(1:(x-1)); reshape(cpu, lproportions, 1); ... 
				qn.S{e}((x+1):qn.M)];
		end
	else
		newS = [qn.S(1:(x-1)); reshape(cpu, lproportions, 1); qn.S((x+1):qn.M)];
	end
	
	% calculate new rates
	if isa(qn,'CMCQNCSRE')
		newRates = cell(qn.E, 1);
		for e=1:qn.E
			newRates{e} = [qn.rates{e}(1:(x-1),:); repmat(qn.rates{e}(x,:), lproportions, 1); ...
				qn.rates{e}((x+1):qn.M,:)];
			for ii=1:lproportions
				newRates{e}(x+ii-1,:) = newRates{e}(x+ii-1,:) .* proportions(ii) ./ cpu(ii);
			end
		end
	else
		newRates = [qn.rates(1:(x-1),:); repmat(qn.rates(x,:), lproportions, 1); ...
			qn.rates((x+1):qn.M,:)];
		for ii=1:lproportions
			newRates(x+ii-1,:) = newRates(x+ii-1,:) .* proportions(ii) ./ cpu(ii);
		end
	end

	% calculate new sched
	newSched = [qn.sched(1:(x-1)); repmat(qn.sched(x), lproportions, 1); ...
		qn.sched((x+1):qn.M)];
	
	% calculate new P
	if isa(qn,'CMCQNCSRE')
		newP = cell(qn.E, 1);
		for e=1:qn.E
			newP{e} = zeros(qn.K*newM, qn.K*newM);
			for kk=1:qn.K
				Pk = qn.P{e}(kk:qn.K:end, kk:qn.K:end);
				newPk = [Pk(1:(x-1),:); repmat(Pk(x,:), lproportions, 1); Pk((x+1):qn.M,:)];
				newPk = [newPk(:,1:(x-1)) repmat(newPk(:,x), 1, lproportions) ... 
					newPk(:,(x+1):qn.M)];
				for ii=1:lproportions
					newPk(:, x + ii - 1) = newPk(:, x + ii - 1)*proportions(ii);
				end
				newP{e}(kk:qn.K:end, kk:qn.K:end) = newPk;
			end
		end
	else
		newP = zeros(qn.K*newM, qn.K*newM);
		for kk=1:qn.K
			Pk = qn.P(kk:qn.K:end, kk:qn.K:end);
			newPk = [Pk(1:(x-1),:); repmat(Pk(x,:), lproportions, 1); Pk((x+1):qn.M,:)];
			newPk = [newPk(:,1:(x-1)) repmat(newPk(:,x), 1, lproportions) ... 
				newPk(:,(x+1):qn.M)];
			for ii=1:lproportions
				newPk(:, x + ii - 1) = newPk(:, x + ii - 1)*proportions(ii);
			end
			newP(kk:qn.K:end, kk:qn.K:end) = newPk;
		end
	end
	
	% calculate new refNodes
	newRefNodes = qn.refNodes;
	newRefNodes(newRefNodes>x) = newRefNodes(newRefNodes>x) + lproportions - 1;
	
	% calculate new nodeNames
	newNodeNames = [qn.nodeNames(1:(x-1)); ... 
		repmat(qn.nodeNames(x), lproportions, 1); qn.nodeNames((x+1):qn.M)];
	for ii=1:lproportions
		newNodeNames(x + ii - 1) = strcat(newNodeNames(x + ii - 1), ... 
			['_r' num2str(resNum(ii))]);
	end
	
	% create new queuing network
	if isa(qn,'CMCQNCSRE')
		newQN = CMCQNCSRE(newM, qn.K, qn.E, qn.N, qn.ENV, newS, newRates, newSched, ...
			newP, qn.NK, qn.classMatch, newRefNodes, qn.resetRules, ...
			newNodeNames, qn.classNames, qn.adhocResetRules);
	else
		newQN = CMCQNCS(newM, qn.K, qn.N, newS, newRates, newSched, newP, qn.NK, ...
			qn.classMatch, newRefNodes, newNodeNames, qn.classNames);
	end
end
