%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^(-14)]
alfa=1.158*l/l0;
w=50;m=100;
d=w/m;
f=d/2:d:w-d/2;
k=max(size(f));p=ones(1,k);

aa=[5.4318e-001;
  3.4746e-001;
  2.5545e-001;
  1.9894e-001;
  1.6026e-001;
  1.3210e-001;
  1.1074e-001;
  9.4066e-002;
  8.0749e-002;
  6.9926e-002;
  6.1003e-002;
  5.3558e-002;
  4.7283e-002;
  4.1948e-002;
  3.7378e-002;
  3.3436e-002;
  3.0016e-002;
  2.7033e-002;
  2.4418e-002;
  2.2116e-002;
  2.0082e-002;
  1.8277e-002;
  1.6670e-002;
  1.5235e-002;
  1.3951e-002;
  1.2797e-002;
  1.1758e-002;
  1.0821e-002;
  9.9736e-003;
  9.2055e-003;
  8.5080e-003;
  7.8734e-003;
  7.2951e-003;
  6.7671e-003;
  6.2843e-003;
  5.8422e-003;
  5.4368e-003;
  5.0644e-003;
  4.7220e-003;
  4.4067e-003;
  4.1160e-003;
  3.8477e-003;
  3.5998e-003;
  3.3705e-003;
  3.1582e-003;
  2.9614e-003;
  2.7789e-003;
  2.6093e-003;
  2.4518e-003;
  2.3052e-003;
  2.1687e-003;
  2.0415e-003;
  1.9229e-003;
  1.8122e-003;
  1.7089e-003;
  1.6123e-003;
  1.5219e-003;
  1.4374e-003;
  1.3582e-003;
  1.2840e-003;
  1.2144e-003;
  1.1491e-003;
  1.0878e-003;
  1.0303e-003;
  9.7615e-004;
  9.2527e-004;
  8.7740e-004;
  8.3234e-004;
  7.8990e-004;
  7.4991e-004;
  7.1221e-004;
  6.7665e-004;
  6.4309e-004;
  6.1142e-004;
  5.8150e-004;
  5.5323e-004;
  5.2651e-004;
  5.0124e-004;
  4.7734e-004;
  4.5472e-004;
  4.3330e-004;
  4.1301e-004;
  3.9379e-004;
  3.7557e-004;
  3.5829e-004;
  3.4191e-004;
  3.2636e-004;
  3.1161e-004;
  2.9760e-004;
  2.8429e-004;
  2.7165e-004;
  2.5963e-004;
  2.4821e-004;
  2.3735e-004;
  2.2701e-004;
  2.1718e-004;
  2.0782e-004;
  1.9891e-004;
  1.9043e-004;
  1.8234e-004];

 du=aa';

%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
   %hc3(n)=exp(-alfa)*f(n)^.5);
   hc3(n)=du(n);
hx(n)=beta(g)*f(n)^1.5;
end
p1=100;
la=w/p1
n0=10^(-11);
for n=1:100
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^(-14)') 
%gtext('beta=10^(-9.5)') 
%gtext('beta=10^(-10.5)') 
