function [alfa_x2,vysledek_p,Loss,Chi2,Chi3,ch2,ch3,citatel]=ref7(R,lambda,zeta,pr,n2)
%[alfa_x2]=ref7(R,lambda);
%Vypocet 2*alfa ohnuteho vlakna typu SMF-28 podle [7]:
%Faustini, L., Martini, G.: Bend Loss in Single-Mode Fibers. IEEE Journal
%of Ligthwave technology, vol. 15, no. 4, Apr. 1997, pp. 671-679.

%Definovani velicin
a=4e-6;     b=62.5e-6;  c=125e-6; %nevyuzito (polomer 1.ochranne vrstvy)
n1=1.4504;   
n3=1.4786;

%Konstanty sireni
k0=2*pi./lambda;     k1=k0.*n1;   k2=k0.*n2;   k3=k0.*n3;

%!! Pri integraci je zeta vektor => vsechny funkce vraceji take vektor !!
%   R_ je aktualne zpracovavana hondota (vektoru) R.
%   Funkce jsou vnorene (nested) a prebiraji hodnoty k0, k1, k2, beta atd.

%Definovani X2,3(x,zeta)    (vztah 2.1)
function [vysl]= X2(x,zeta)
    vysl=(R_./2./k2^2).^(2/3)*(beta.^2+zeta.^2-k2.^2*(1+2.*x./R_));
end

function [vysl]= X3(x,zeta)
    vysl=(R_./2./k3^2).^(2/3)*(beta.^2+zeta.^2-k3.^2.*(1+2.*x./R_));
end

%Definovani chi2,3 (zeta)   (vztah 13.1)
%chi_q je vlastne jen hranata zavorka z X_q pro x=b
function [vysl]=chi2(zeta)
    vysl=-beta.^2-zeta.^2+k2.^2.*(1+2.*b./R_);
end
function [vysl]=chi3(zeta)
    vysl=-beta.^2-zeta.^2+k3.^2.*(1+2.*b./R_);
end

%Definovani theta (zeta) (vztah 13.2)
%exponent 3/2 kompenzuje 2/3 z X_q
function [vysl]=theta(zeta)
    vysl=2/3.*R_./2./k2.^2*(-beta.^2-zeta.^2+k2.^2.*(1+2.*b./R_)).^(3/2)+pi()/4;
end

%Integrand ve vztahu (14), pro prehlednost je rozlozen na 3 soucinitele
function [vysl]=integrand(zeta)
    pomocne=sqrt(gamma.^2+zeta.^2);
    cinitel1=exp(-a.*pomocne)./pomocne;
    cinitel2=airy(0,X2(0,zeta))./airy(2,X2(a,zeta));
    cinitel3=(sqrt(chi2(zeta)).*sqrt(chi3(zeta)))./((chi2(zeta).*cos(theta(zeta)).^2+chi3(zeta).*sin(theta(zeta)).^2));
    vysl=cinitel1.*cinitel2.*cinitel3;
end

% Zadefinovani disperzni rovnice
function [u,v,norm_beta]=disper_LP01(V)
    
    f=inline('abs((besselj(0,u))/(u.*besselj(1,u))-(besselk(0,sqrt(V.^2-u.^2)))/(sqrt(V.^2-u.^2)*besselk(1,sqrt(V.^2-u.^2))))','u','V'); %definice disperzni rovnice --> dostaneme fci f(u,v)
    u=fminbnd(@(u)f(u,V),0,V); % vypocet v (znam-li u) na intervalu <0,V>
    v=sqrt(V.^2-u.^2);
    zkouska=inline('abs((besselj(0,u))/(u.*besselj(1,u))-(besselk(0,v))/(v.*besselk(1,v)))','u','v'); %definice disperzni rovnice --> dostaneme fci f(u,v)
    %fprintf('---------------------------------------------------------------------------------------- \n')
    fprintf('Chyba reseni disperzni rovnice: %i.\n', zkouska(u,v))
    norm_beta=(v./V)^2;
    
end
%                   ******************
%                   * Hlavni vypocet *
%                   ******************
%        Nasledujci hodnoty jsou behem integrace pevne
NA=sqrt(n1^2-n2^2);                 %numericka apertura
V=2*pi()*a./lambda*NA;               %norm. frekvence, neplest s 'v'
for (i=1:length(lambda))
[u(1,i),v(1,i),norm_beta(1,i)]=disper_LP01(V(1,i));     %pro LP01 nalezne z disperzni rovnice u,v a normovane beta
end
nef=sqrt(norm_beta*NA^2+n2^2);      %ef. index lomu
beta=k0.*nef;
kappa=sqrt(k1.^2-beta.^2);    %take rovno u/a
gamma=sqrt(beta.^2-k2.^2);    %take rovno v/a

%kontrolni vypis, normalne neni treba
% clc
 %fprintf('SMF-28, ciselne vysledky\n')
 fprintf('lambda= %4.0f\n ',lambda*1e9)
 %fprintf('a= %4.0f\t b=%4.0f\t c=%4.0f\n',a*1e6,b*1e6,c*1e6)
 %fprintf('V= %4.2f\n',V)
 fprintf('NA= %5.3f\n',NA)
 %fprintf('b= %5.2f\t u=%5.2f\t v=%5.2f\n',norm_beta,u,v)
 %fprintf('n1= %5.3f\t n2=%5.3f\t nef= %5.3f\n',n1,n2,nef)
 %fprintf('k0= %5.2e\t k1= %5.2e\t beta= %5.2e\t kappa= %5.2e\t gamma= %5.2e\n',k0,k1,beta,kappa,gamma)

%cyklus pres vechny prvky (vektoru) polomeru R

Chi2=(-((beta^2)+(zeta.^2)-(k2^2)*(1+2*b/R)));         % vztah (13.1) v [ref7]
Chi3=(-((beta^2)+(zeta.^2)-(k3^2)*(1+2*b/R)));         % vztah (13.1) v [ref7]
citatel=sqrt(Chi2).*sqrt(Chi3);                        % citatel tretiho integrandu vztahu (14) v [ref7]

for(i=1:length(Chi2)-1)                                % hledani horni meze integralu, tj. chi2=0
if(Chi2(1,i)>0 & Chi2(1,i+1)<0)
    horni_mez=i;
    ch2=horni_mez;
    horni_mez=zeta(1,horni_mez);
    nula2=Chi2(1,i);
end
end

for(i=1:length(Chi3)-1)                                 % hledani horni meze integralu, tj. chi2=0
if(Chi3(1,i)>0 & Chi3(1,i+1)<0)
    ch3=i;
    %ch3=zeta(1,ch3);
    nula3=Chi3(1,i);
end
end


for prvek=1:max(size(R))    
    R_=R(prvek);                                        %aktualne zpracovavany polomer, na nem zavisi vsechny funkce
    %integral=quad(@integrand,0,limit_integrace*gamma);  %namisto nekonecna je horni mez=1*gamma, pro 10*gamma jiz hlasi problemy
    for (i=1:length(lambda))
    integral=romberg(@integrand,0,horni_mez,pr); %romberguvo metoda integrovani
    end
    vysledek=2*kappa.^2./beta./V.^2./besselk(1,v).^2.*integral; %vztah (14)
    vysledek_p=vysledek(pr,pr-1); %kontrola posledni a predposledni hodnoty rombergrovy tabulky
    vysledek=vysledek(pr,pr); %pri rombergove integrace vybere posledni
    %clen tabulky
    
   % koef=10*log10(exp(1));                              %prepocet Np/m na dB/m
    koef=20/log(10); 
    fprintf('Pro R= %4.1f mm je 2*alfa= %5.2f Np/m, tj. %5.2f dB/m\n',R_*1e3,vysledek,vysledek*koef);
    alfa_x2(prvek)=vysledek;
    %Loss(1,prvek)=vysledek.*koef; %odkomentuj pro srovnani
    Loss=vysledek.*koef;%odkomentuj pro main
end

end %funkce ref7
