w=500000;m=100;%l=9000;l0=18000;
%beta=[10^(-7.5) 10^(-8) 10^(-9) 10^(-9.5) 10^(-10.5)];
beta=[1.0*10^(-13)]
%alfa=1.4053*l/l0;

d=w/m;
f=d/2:d:w-d/2;
k=max(size(f));p=ones(1,k);

%l=linspace(.000985,.000872,127);
%c=(.083*10^(-6)).*ones(1,k);hc=ones(1,k);

%for t=1:k
%r(t)=438.29996+33.3*(.00001*f(t))+5.012*(.00001*f(t))^2-.42*(.00001*f(t))^3+.012*(.00001*f(t))^4-.001*(.00001*f(t))^5;

%gamma(t)=sqrt((r(t)+i*2*pi*f(t)*l(t))*(i*2*pi*f(t)*c(t)));
%end
%alfa=unwrap(real(gamma));
%for t=1:k

%hc(t)=exp(-alfa(t)*2*9000*(.3048/1609.3));
%hc1(t)=exp(-1.4*.5*f(t)^.5);

%end
for g=1:1
for n=1:k
%hc(n)=exp(-alfa*f(n)^.5);
hx(n)=beta(g)*f(n)^1.5;
end
p1=50;
la=w/p1
n0=10^(-11);
for n=1:k
a(n)=hx(n)*(hx(n)+hc3(n));n1=n0*d;
b(n)=n1*(2*hx(n)+hc3(n));c(n)=n1*(n1-d*(1/la)*hc3(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;
while dif<0
e=.02;
la=(1-e)*la;
for n=1:k
a(n)=hx(n)*(hx(n)+hc3(n));n1=n0*d;
b(n)=n1*(2*hx(n)+hc3(n));c(n)=n1*(n1-(1/la)*d*hc3(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;
%plot(f1,u)
hold on
semilogy(f1,p2)
cap(g)=0;
for n=1:k
cap(g)= cap(g)+d*log2(1+p(n)*hc3(n)/(n1+hx(n)*p(n)));
end
end
xlabel ('frecventa in khz');
ylabel ('densitatea spectrala per canal in mW')
title('Distributia puterii pentru o capacitate maxima')
%gtext('beta=10^(-7.5)') 
%gtext('beta=10^(-8)') 
gtext('beta=1.4*10^(-9)') 
%gtext('beta=10^(-9.5)') 
%gtext('beta=10^(-10.5)')