%w=400000;m=127;
l=6000;l0=18000;
%beta=[10^(-7.5) 10^(-8) 10^(-9) 10^(-9.5) 10^(-10.5)];
beta=[1.0*10^(-8.5)]
alfa=1.158*l/l0;
w=400;m=80;
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
du=[ 1.2787e-002;
  1.0482e-002;
  7.3794e-003;
  5.1206e-003;
  3.7012e-003;
  2.8150e-003;
  2.2390e-003;
  1.8452e-003;
  1.5623e-003;
  1.3500e-003;
  1.1850e-003;
  1.0534e-003;
  9.4589e-004;
  8.5635e-004;
  7.8047e-004;
  7.1516e-004;
  6.5817e-004;
  6.0786e-004;
  5.6300e-004;
  5.2266e-004;
  4.8615e-004;
  4.5288e-004;
  4.2233e-004;
  3.9408e-004;
  3.6780e-004;
  3.4322e-004;
  3.2011e-004;
  2.9826e-004;
  2.7750e-004;
  2.5769e-004;
  2.3877e-004;
  2.2071e-004;
  2.0352e-004;
  1.8721e-004;
  1.7176e-004;
  1.5713e-004;
  1.4330e-004;
  1.3026e-004;
  1.1803e-004;
  1.0662e-004;
  9.6004e-005;
  8.6149e-005;
  7.7008e-005;
  6.8542e-005;
  6.0729e-005;
  5.3557e-005;
  4.7010e-005;
  4.1060e-005;
  3.5671e-005;
  3.0806e-005;
  2.6434e-005;
  2.2531e-005;
  1.9074e-005;
  1.6037e-005;
  1.3388e-005;
  1.1093e-005;
  9.1203e-006;
  7.4383e-006;
  6.0193e-006;
  4.8360e-006;
  3.8604e-006;
  3.0650e-006;
  2.4236e-006;
  1.9122e-006;
  1.5097e-006;
  1.1975e-006;
  9.5892e-007;
  7.7973e-007;
  6.4780e-007;
  5.5327e-007;
  4.8845e-007;
  4.4757e-007;
  4.2649e-007;
  4.2234e-007;
  4.3323e-007;
  4.5795e-007;
  4.9577e-007;
  5.4623e-007;
  6.0898e-007;
  6.8371e-007];
duc=du';
for l=1:k
   duu(l)=duc(l)*((f(l)^2/(f(l)^2+1)))^2;
   end

for g=1:1
for n=1:k
%hc3(n)=exp(-alfa*f(n)^.5);
hc3(n)=duu(n);
hx(n)=beta(g)*f(n)^1.5;
end
p1=100;
la=w/p1
n0=10^(-11);
for n=1:k
a(n)=hx(n)*(hx(n)+hc3(n));n1=n0*d*10^3;
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.0*10^(-13)') 
%gtext('beta=10^(-9.5)') 
%gtext('beta=10^(-10.5)') 
for  n=1:k; ca(n)=log2(1+p(n)*hc3(n)/(n1+hx(n)*p(n))); end;
%plot(f,ca)