% Author: Zukang Shen and Ian C. Wong  
% Copyright (C) 2004   Ian C. Wong
% 
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
% 
% You may reach the author at wongic@mail.utexas.edu.
% Or visit his website at www.ece.utexas.edu/~iwong
function [p]= shenpowerallo(channel,suballo,N,K,Ptot,noise,gamma)
global shenoops
global rootiter
global rootfnct

n = sum(suballo');


for user=1:K
   used=[];
	for j=1:N
   	if suballo(user,j)>0
      	used=[used channel(user,j)];
   	end
	end
	usedlength=length(used);
	for i=1:length(used)
   	[minvalue minindex]=min(used);
   	usedordered(user,i)=minvalue;
   	used=[used(1:minindex-1) used(minindex+1:length(used))];
	end
   usedordered(user,:)=usedordered(user,:)/noise;	
   
   h(user)=usedordered(user,1);
   vtemp=0;
   for m=2:usedlength
   	vtemp=vtemp+(usedordered(user,m)-usedordered(user,1))/(usedordered(user,1)*usedordered(user,m));   
   end
   v(user)=vtemp;
   
   ytemp=1;
   for m=2:usedlength
   	ytemp=ytemp*usedordered(user,m)/usedordered(user,1);   
   end
   y(user)=ytemp;
   
end

coef(1)=1;
power(1)=1;
for i=2:K
    coef(i)=n(i)/(h(i)*y(i)^(1/n(i)))*(h(1)*y(1)^(1/n(1))/n(1))^(n(1)*gamma(i)/(n(i)*gamma(1)));
    power(i)=n(1)*gamma(i)/(n(i)*gamma(1));
end

for i=1:101
   guessx=0.01*(i-1)*Ptot;
   guess(i)=guessx;
   for j=2:K
   	guess(i)=guess(i)+coef(j)*guessx^power(j);   
   end
   
   guess(i)=guess(i)-Ptot;   
end
[guessvalue,guessindex]=min(abs(guess));

[p(1) f e o]=fzero(@rootfinding,0.01*(guessindex-1)*Ptot,[],Ptot,coef,power);

if isnan(p(1))==1
   for i=1:K      
      p(i)=Ptot/N*n(i);
   end
else 
   rootiter = rootiter + o.iterations;
   rootfnct = rootfnct + o.funcCount;
   for i=2:K
      p(i)=coef(i)*p(1)^power(i);
   end
end