clear
pack

f1=[0 2500 10000 20000 40000 100000 150000 200000 400000 500000 650000 750000 850000 1000000 1100000];

f=4000:4000:1100000;d=4000;
r1=[280 280.007 280.11 280.44  281.748 290.433 302.070 316.393 383.562 417.427 465.785 496.218 525.274 566.521 592.628];
l1=10^(-6).*[587.132 587.075 586.738 586.099 584.443 577.878 571.525 564.889 539.437 528.409 514.272 506.304 499.343 490.494 485.481];
%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)*43*ones(size(f));
g=zeros(size(f));z1=zeros(size(f));y=zeros(size(f));d1=1.68; d5=1.112;d7=.56*.3048; %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=[0  2500 10000 20000 40000  100000 150000 200000 400000  500000 650000 750000 850000 1000000 1100000];
f=4000:4000:1100000;%r11 este rez. per. de .5mm;
r11=[179 179.015 179.244 179.97 182.79 199.608 218.721 239.132 316.23 350.243 396.528 424.835 451.480 488.857 512.317];
l11=10^(-6)*[673.574 673.466 672.923 671.98 669.716 661.677 654.622 647.735 624.309 615.202 604.192 598.261 593.217 586.966 583.495];
r12=[ 113 113.028 113.442 114.737 119.523 143.115 164.938 185.689 255.575 284.753 323.844 347.542 369.758 400.816 420.264];
l12=10^(-6)*[699.258 697.943 693.361 687.008 674.593 642.718 622.050 605.496 563.845 551.323 537.868 531.212 525.813 519.402 515.956];
r63=spline(f11,r12,f);
l63=spline(f11,l12,f);
c63= 50*10^(-9).*ones(size(f));
g63=zeros(size(f));
d63=.5;
for k=1:length(f);
z63(k)=r63(k)+2*pi*f(k)*l63(k)*i;
y63(k)=g63(k)+2*pi*f(k)*c63(k)*i;
zc63(k)=sqrt(z63(k)/y63(k));
gamma63(k)=sqrt(z63(k)*y63(k));
B63(k)=zc63(k)*sinh(gamma63(k)*d63);
A63(k)=cosh(gamma63(k)*d63);
D63(k)=A63(k);
C63(k)=(1/zc63(k))*sinh(gamma63(k)*d63);
e63(k)=A63(k)*D63(k)-B63(k)*C63(k);
end





r2=spline(f11,r11,f);
l2=spline(f11,l11,f);
g2=zeros(size(f));
c2= 50*10^(-9).*ones(size(f));
%g=zeros(size(f));

z11=zeros(size(f));y=zeros(size(f)); d2=1.25;%portiunea cu d=.5mm;
d3=.560;%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));
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=150;hc3=ones(size(f));
for k=1:length(f);
m1=[1 150;0 1];
m2=[A(k) B(k);C(k) D(k)];
m63=[A63(k) B63(k);C63(k) D63(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*m7*m3*m9;
hc(k)= rr/(m5(1,1)*rr+m5(1,2));
hc3(k)=(abs(hc(k)))^2;
hc4=10*log10(hc3);
end
%semilogy(f,hc3)
plot(f,hc4)
hold on
betaa=10^(-13);
for n=1:length(f)
%hc(n)=exp(-alfa(o)*f(n)^.5);
hx(n)=betaa*f(n)^1.5;
end
semilogy(f,hx)


cap=0;
for n=1:k
cap= cap+d*log2(1+hc3(n)/(hx(n)));
end
n0=10^(-12.340);
nn=256; fs=2048000; kfext=7.74*10^(-21);t=7
ff=f(1:1,t:t+nn);
hc31=hc3(1:1,t:t+nn);
fext=ones(1,length(ff));
d0=3000+4000+1500;
c=1.*ones(1,length(ff));
for j=1:length(ff)
fext(j)= kfext*hc31(j)*d1*(ff(j))^2*(49^(.6))*1000*3.28;
end
f=4000:4000:1100000;f0=2.208*10^6;Kg=.1104
fh=1.104*10^6;
fl=4000;
fh1=f0/2;
fh=25875;
alfa=11.96;
alfa1=7.09;
m=length(f)
hpf=ones(size(f));
si=ones(size(f));
lpf=ones(size(f));

for i=1:length(f)
hpf(i)=(f(i)^alfa1+fl^alfa1)/(f(i)^alfa1+fh^alfa1);
si(i)=(sin(pi*f(i)/f0)/(pi*f(i)/f0))^2;




lpf(i)=fh1^alfa/(f(i)^alfa+fh1^alfa);
psd(i)=Kg*(2/f0)*si(i)*hpf(i)*lpf(i);
ps(i)=10*log10(1000*psd(i));
end
%plot(f,ps)
%axis([0 f0 -150 0]);grid
f1=4000:4000:40000;
hpf1=hpf(1:1,1:10);
psd1=psd(1:1,t:t+nn);
adslfext=fext.*psd1;
adslf=10*log10(1000.*adslfext);
plot(ff,adslf);grid
