% Author: 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]= wongpowerallo(channel,suballo,N,K,Ptot,noise,gamma)
global oops
Nk = sum(suballo');
H = channel.*suballo/noise;
[Hsorted,I] = sort(H');
Hsorted = Hsorted';
n = [2:N];
k = [2:K];
V = zeros(1,K);
W = ones(1,K);
for k = 1 : K,
id = find(Hsorted(k,:));
Hk1 = Hsorted(k,id(1));
if k == 1,
H11 = Hk1;
end;
if length(id) > 1,
for n = id(2) : N,
V(k) = V(k) + (Hsorted(k,n) - Hk1)/(Hk1*Hsorted(k,n));
W(k) = W(k)*(Hsorted(k,n)/Hk1)^(1/(Nk(k)));
a(k) = -Nk(1)/Nk(k)*Hk1*W(k)/(H11*W(1));
b(k) = Nk(1)/(H11*W(1))*(W(k) - W(1) + H11*V(1)*W(1)/Nk(1) - Hk1*V(k)*W(k)/Nk(k));
end;
else,
V(k) = 0;
W(k) = 1;
a(k) = -Nk(1)/Nk(k)*Hk1*W(k)/(H11*W(1));
b(k) = Nk(1)/(H11*W(1))*(W(k) - W(1) + H11*V(1)*W(1)/Nk(1) - Hk1*V(k)*W(k)/Nk(k));
end;
end;
a(1) = 0;
b(1) = Ptot;
A = diag(a) + not([0; ones(K-1,1)]*[0 ones(1,K-1)]);
% p = (A\b')'
p = zeros(length(a),1);
bcut = b(2:length(b));
acut = a(2:length(a));
p(1) = (Ptot-sum(bcut./acut))/(1-sum(1./acut));
p(2:length(a)) = ((bcut)-p(1))./acut;
if any(p < zeros(size(p)) | ~isreal(p)),
oops = oops + 1;
p = 1/K*ones(1,K);
end;
