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;
k=max(size(f));
for n=1:k
hc(n)=exp(-alfa*f(n)^.5);
hx(n)=beta*f(n)^1.5;
end
p1=10^(-2.3);la=w/p1
%la=12647
n0=10^(-11);
for n=1:k
a(n)=hx(n)*(hx(n)+hc(n));n1=n0*d*1000;
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);
end
s=0;
for t=1:k
s=s+p(t);
end
dif=s-p1;
e=.08
while dif<0
la=(1-e)*la;
for n=1:k
a(n)=hx(n)*(hx(n)+hc(n));n1=n0*d*1000;
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);
end
q=zeros(1,k);y=q>=c;u=y.*p;
s=0;
for t=1:k
s=s+u(t);
end
dif=s-p1;
end
x=1:1:k;f1=d.*x;p2=(1/1).*u;
%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^1.7 mW,beta=10^(-9),model f^.5,m=1000 canale ')
cap=0
for n=1:k
cap= cap+d*1000*log2(1+p(n)*hc(n)/(n1+hx(n)*p(n)));
end