w=100000;m=1000;l=600;l0=18000;beta=10^(-9);
alfa=1.158*l/l0;d=w/m;
f=d/2:d:w-d/2;r=44; pa=21.9908
%fa=f(r):d:W-d;
k=max(size(f));hc=ones(1,k);
for n=r:k
hc(n)=exp(-alfa*f(n)^.5);
%hx(n)=beta*f(n)^1.5;
end
p1a=10^(2.7)-pa;%la=(w-f(r))/(2*p1a)
la=98.3
n0=10^(-10.99);
n1=n0*d*1000;p=zeros(1,k);
for n=r:k
%a(n)=hx(n)*(hx(n)+hc(n));

%b(n)=n1*(2*hx(n)+hc(n));c(n)=n1*(n1-(1/la)*hc(n));
%p(n)=(1/(2*a(n)))*(-b(n)+(b(n)^2-4*a(n)*c(n))^.5);
p(n)=d/(2*la)-n1/(4*hc(n));
end
s=0;
for t=r:k
s=s+p(t);
end
dif=s-p1a;
e=.0005
while dif<0
la=(1-e)*la;
n1=n0*d*1000;
for n=r:k
%a(n)=hx(n)*(hx(n)+hc(n));

%b(n)=n1*(2*hx(n)+hc(n));c(n)=n1*(n1-(1/la)*hc(n));
%p(n)=(1/(2*a(n)))*(-b(n)+(b(n)^2-4*a(n)*c(n))^.5);
p(n)=d/(2*la)-n1/(4*hc(n));
end
%q=zeros(1,k);y=p>=q;u=y.*p;
s=0;
q=zeros(1,k);
for t=r:k
qq=[p(t) 0];q(t)=max(qq);
s=s+q(t);
end
dif=s-p1a;
end
%[i,j]=min(y)
x=1:1:k;f1=d.*x;p2=(1/1).*q;
%s=0;
%for t=1:k
%s=s+p(t);
%end
%dif=s-p1;
%end
%x=1:1:k;
plot(f1,p2)

xlabel ('frecventa in khz');
ylabel ('densitatea spectrala per canal in mW')
title('Distributia puterii pentru o capacitate maxima')
gtext('p=10^2.7 mW,awgn,model f^.5,m=1000 canale ')
cap=0
for n=r:k
cap= cap+d*500*log2(1+2*p(n)*hc(n)/(n1));
%dr(n)=(hc(n)-2*hx(n))/(hx(n)^2);
end
%v=2*p2/n1-dr;
%yy=v>=q;
%[z,r]=max(yy)
%dd=p2(1:1,1:r-1);
%pa=sum(pa)