w=400;m=160;l=9000;l0=18000;beta=10^(-10);
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=100;la=1/14.5
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(n);
end
dif=s-p1
x=1:1:k;
plot(x,p)
cap=0
for n=1:k
cap= cap+d*1000*log2(1+p(n)*hc(n)/(n1+hx(n)*p(n)));
end