clear
pack

f1=[1000 4000 10000 20000 40000 100000 150000 200000 400000 500000 1000000];

f=4000:4000:1060000;d=4000;
r1=[268 268 268 269  271 282 295 312 390 425 603];
l1=10^(-6).*[680 680 678 675 669 650 642 635 619 608 574];
%c5=polyfit(f1,r1,3);

%c15=polyfit(f1,l1,3);
%r=polyval(c5,f);
%l=polyval(c15,f);
r=spline(f1,r1,f);
l=spline(f1,l1,f);
c=10^(-9)*45.5.*ones(size(f));
g=zeros(size(f));z1=zeros(size(f));y=zeros(size(f));d1=3*.3048; d5=6*.3048;d7=1.5*.03048; %ultima portiune;
for k=1:length(f);
z1(k)=r(k)+2*pi*f(k)*l(k)*i;
y(k)=g(k)+2*pi*f(k)*c(k)*i;
zc(k)=sqrt(z1(k)/y(k));
gamma(k)=sqrt(z1(k)*y(k));
alfa(k)=real(gamma(k));
beta(k)=imag(gamma(k));
 B(k)=zc(k)*sinh(gamma(k)*d1);
B5(k)=zc(k)*sinh(gamma(k)*d5);
B7(k)=zc(k)*sinh(gamma(k)*d7);





A(k)=cosh(gamma(k)*d1);
A5(k)=cosh(gamma(k)*d5);
A7(k)=cosh(gamma(k)*d7);
D7(k)=A7(k);
D5(k)=A5(k);
D(k)=A(k);
C(k)=(1/zc(k))*sinh(gamma(k)*d1);
C5(k)=(1/zc(k))*sinh(gamma(k)*d5);
C7(k)=(1/zc(k))*sinh(gamma(k)*d7);


e(k)=A(k)*D(k)-B(k)*C(k);

end
f11=[1000  4000 10000 20000 40000  100000 150000 200000 400000  500000 1000000];
f=4000:4000:1060000;
r11=[172 172 172 173 175 190 207 227 302 334 470];
%r11=(1/1.6093).*[277.2 277.5 278  286.8 308.4 337.2 431.6 541.7];
l11=10^(-3)*[.679 .678 .678 .675 .667 .646 .637 .629 .603 .592 .559];
%l11=(1/1.6093)*10^(-3).*[.986 .984 .982 .958 .935 .920 .888 .857];
g11=(1/1.6093)*10^(-6).*[.115 .4 .853 1.53 3.258 6.32 8.993 11.95 18.44  25.633 50.0];
%c5=polyfit(f1,r1,3);

%c15=polyfit(f1,l1,3);
%r=polyval(c5,f);
%l=polyval(c15,f);
r2=spline(f11,r11,f);
l2=spline(f11,l11,f);
g2=spline(f11,g11,f);
c2=(1/1.6093)*10^(-9)*83.*ones(size(f));
%g=zeros(size(f));

z11=zeros(size(f));y=zeros(size(f)); d2=.0002;%portiunea cu d=.5mm;
d3=.3048*1.5;%prima derivatie;
d4=.3048*1.5;%a 2-a derivatie;
d6=.3048*1.5;% a 3-a derivatie;
for k=1:length(f);
z11(k)=r2(k)+2*pi*f(k)*l2(k)*i;
y2(k)=g2(k)+2*pi*f(k)*c2(k)*i;
zc2(k)=sqrt(z11(k)/y2(k));
gamma2(k)=sqrt(z11(k)*y2(k));
%for k=1:length(f);
B2(k)=zc2(k)*sinh(gamma2(k)*d2);
A2(k)=cosh(gamma2(k)*d2);
D2(k)=A2(k);
C2(k)=(1/zc2(k))*sinh(gamma2(k)*d2);
e2(k)=A2(k)*D2(k)-B2(k)*C2(k);
end
for k=1:length(f);
zg(k)=zc(k)*cosh(gamma(k)*d3)/sinh(gamma(k)*d3);
zg2(k)=zc(k)*cosh(gamma(k)*d4)/sinh(gamma(k)*d4);
zg3(k)=zc(k)*cosh(gamma(k)*d6)/sinh(gamma(k)*d6);



mod(k)=abs(zg(k));fi(k)=angle(zg(k));xr(k)=mod(k)*cos(fi(k));
yi(k)=mod(k)*sin(fi(k));
end
m1=ones(2,2);
m2=ones(2,2);
m3=ones(2,2);
m4=ones(2,2);
m5=ones(2,2);
m6=ones(2,2);
m7=ones(2,2);
m9=ones(2,2);
m10=ones(2,2);
rr=105;hc3=ones(size(f));
for k=1:length(f);
m1=[1 105;0 1];
m2=[A(k) B(k);C(k) D(k)];
m3=[1 0;1/zg(k) 1];
m4=[A2(k) B2(k);C2(k) D(k)];
m6=[1 0;1/zg2(k) 1];
m10=[1 0;1/zg3(k) 1];
m7=[A5(k) B5(k);C5(k) D5(k)];
m9=[A7(k) B7(k);C7(k) D7(k)];

m5=m1*m2*m3*m7*m6*m9*m10;
hc(k)= rr/(m5(1,1)*rr+m5(1,2));
hc3(k)=(abs(hc(k)))^2;
end
semilogy(f,hc3)
betaa=10^(-13);
for n=1:k
%hc(n)=exp(-alfa(o)*f(n)^.5);
hx(n)=betaa*f(n)^1.5;
end
cap=0;
for n=1:k
cap= cap+d*log2(1+hc3(n)/(hx(n)));
end
n0=10^(-13.8);
nn=256; fs=2048000; kfext=2*10^(-21.49);t=10
ff=f(1:1,t:nn);
hc31=hc3(1:1,t:nn);
fext=ones(1,length(ff));
d=3000+6000+1500;
c=1.*ones(1,length(ff));
for j=1:length(ff)
fext(j)= kfext*hc31(j)*d*(ff(j))^2*(49^(.6));
end
ic=0;icmax=14;
ga=1;
l=length(ff);
snr=zeros(1,l);b=zeros(1,l); bb=zeros(1,l);
for j=1:length(ff)
snr(j)= c(j)*hc31(j)/(n0*fs/512+c(j)*fext(j));
end
 uc=l;
btot=1;
btar=432;

q=(14.5-10*log10(3))/10;gaa=10^q;
while (btot~=btar)&(ic<=icmax)


for i=1:l
uc=l;
b(i)= log2(1+snr(i)/(gaa*ga));
bb(i)= round(b(i));
diff (i)=b(i)-bb(i);
if bb(i)>=10
bb(i)=10;
end
if (bb(i)<=1)
bb(i)=0;
uc=uc-1;
end
end

btot= sum(bb); ga=ga*2^((btot-btar)/uc);

ic=ic+1;
end
btot
cond=btot-btar
 while cond <=-1
[m,i]=max(diff);

bb(i)=bb(i)+1;
%c(i)=c(i)*(2^bb(i)-1)/(2^(bb(i)-1)-1);
diff=b-bb;

btot=sum(bb);
cond=btot-btar;
end
while cond >=1
[m,i]=min(diff);
if bb(i)>=2
bb(i)=bb(i)-1;
c(i)=c(i)*(2^bb(i)-1)/(2^(bb(i)+1)-1);
end
diff=b-bb;

btot=sum(bb);
cond=btot-btar;


end
for i=1:l


if bb(i)<=1
c(i)=0;
end
end
cc=sum(c);
p=100
aa=p/cc;
dd=aa.*c;
for j=1:l
snr1(j)=dd(j)*hc31(j)/(n0*fs/512+dd(j)*fext(j));
end
for i=1:l
%uc=nn;
b1(i)= log2(1+snr(i)/(gaa*ga));
bb1(i)= round(b(i));
end
