%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^(-11.5)
alfa=1.158*l/l0;
w=500;m=100;
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
aa=[ 8.7068e-003;
  5.5734e-003;
  3.5481e-003;
  2.4058e-003;
  1.7429e-003;
  1.3329e-003;
  1.0624e-003;
  8.7391e-004;
  7.3623e-004;
  6.3200e-004;
  5.5072e-004;
  4.8567e-004;
  4.3241e-004;
  3.8789e-004;
  3.5003e-004;
  3.1731e-004;
  2.8869e-004;
  2.6339e-004;
  2.4081e-004;
  2.2054e-004;
  2.0222e-004;
  1.8556e-004;
  1.7029e-004;
  1.5619e-004;
  1.4310e-004;
  1.3088e-004;
  1.1943e-004;
  1.0870e-004;
  9.8627e-005;
  8.9168e-005;
  8.0281e-005;
  7.1948e-005;
  6.4167e-005;
  5.6941e-005;
  5.0271e-005;
  4.4144e-005;
  3.8540e-005;
  3.3437e-005;
  2.8817e-005;
  2.4665e-005;
  2.0964e-005;
  1.7687e-005;
  1.4805e-005;
  1.2287e-005;
  1.0107e-005;
  8.2392e-006;
  6.6562e-006;
  5.3302e-006;
  4.2318e-006;
  3.3326e-006;
  2.6062e-006;
  2.0280e-006;
  1.5748e-006;
  1.2251e-006;
  9.5949e-007;
  7.6114e-007;
  6.1600e-007;
  5.1266e-007;
  4.4206e-007;
  3.9727e-007;
  3.7316e-007;
  3.6617e-007;
  3.7396e-007;
  3.9508e-007;
  4.2860e-007;
  4.7389e-007;
  5.3031e-007;
  5.9710e-007;
  6.7337e-007;
  7.5804e-007;
  8.4987e-007;
  9.4741e-007;
  1.0490e-006;
  1.1530e-006;
  1.2577e-006;
  1.3613e-006;
  1.4625e-006;
  1.5600e-006;
  1.6524e-006;
  1.7388e-006;
  1.8182e-006;
  1.8900e-006;
  1.9539e-006;
  2.0096e-006;
  2.0570e-006;
  2.0960e-006;
  2.1268e-006;
  2.1497e-006;
  2.1649e-006;
  2.1730e-006;
  2.1742e-006;
  2.1692e-006;
  2.1582e-006;
  2.1420e-006;
  2.1209e-006;
  2.0954e-006;
  2.0661e-006;
  2.0333e-006;
  1.9976e-006;
  1.9593e-006];
du=aa';
for g=1:1
for n=1:k

hc3(n)=du(n);
hx(n)=beta(g)*f(n)^1.5;
end
p1=100.0;
la=w/p1
n0=10^(-8);
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.0*10^(-13)') 
%gtext('beta=10^(-9.5)') 
%gtext('beta=10^(-10.5)') 
