% Vliv ohybu na standardni jednovidove telekomunikacni vlakno SMF-28
% Ing. Jan Sistek, Ph. D., Bc. Martin Valko
% katedra elektromagnetickeho pole
% CVUT v Praze, 2017

clc;
clear all;
close all;
precise=5e-5;                                       % presnost s jakou je urcena horni mez integralu
graf=4;                                             % 1 - utlum na vlnove delce; 
                                                    % 2 - utlum na polomeru ohybu;
                                                    % 3 - utlum na numericke aperture NA, kde parametrem je vlnová délka
                                                    % 4 - utlum na polomeru ohybu, kde parametrem je numerická apertura NA
                                                    % 5 - utlum na polomeru ohybu , kde parametrem je vlnova delka

% NASTAVENI VYCHOZICH HODNOT
% 1
if(graf==1)
R=8.5e-3;
lambda_start=1310e-9;                               % pocatecni vlnova delka
lambda_stop=1330e-9;                                % konecna vlnova delka
step=5e-9;                                          % krok vypoctu
NA=0.118404;                                        % numericka apertura testovaneho vlakna
end
% 2
if(graf==2)
lambda=1310e-9;					    % vlnova delka
R_start=7e-3;			  	        % pocatecni polomer ohybu
R_stop=10e-3;					    % konecny polomer ohybu
step=0.1e-3;					    % krok s jakym se bude menit na intervalu vyse polomer ohybu
NA=0.118404;                        % numericka apertura testovaneho vlakna
end
% 3
if(graf==3)
R=9e-3;						    % polomer ohybu
lambda_start=1550e-9;				    % pocatecni vlnova delka
lambda_stop=1565e-9;				    % konecna vlnova delka
stepl=5e-9;					    % krok zmeny vlnove delky na danem intervalu vyse
NA_start=0.03;					    % pocatecni numericka apertura
NA_stop=0.140;					    % konecna numericka apertura
step=0.001;					    % krok s jakym se bude menit numericka apertura NA na intervalu vyse
end                                                 
% 4
if(graf==4)
NA_v=[0.110 0.113 0.116 0.119 0.122 0.125 0.128];   % definuj pro jake numericke apertury NA chces vypocet provest
lambda=1310e-9;                                     % zadej pro jakou vlnovou delku chces vypocet provest
R_start=7.0e-3;                                     % pocatecni hodnota polomeru ohybu
R_stop=10e-3;                                       % konecna hodnota polomeru ohybu
step=0.1e-3;                                        % krok vypoctu
end
% 5
if(graf==5)
lambda_v=[1310 1470 1550 1570 1625].*10^-9;         % zadej pro jakou vlnovou delku chces vypocet provest
R_start=7e-3;                                       % pocatecni hodnota polomeru ohybu
R_stop=10e-3;                                       % konecna hodnota polomeru ohybu
step=0.05e-3;                                       % krok vypoctu
precise=5e-5;                                       % presnost s jakou je urcena horni mez integralu
NA=0.115804;                                        % numericka apertura testovaneho vlakna
end
% ---------------------------------------------------------------------------
% ---------------------------------------------------------------------------
if(graf==1)                                        % UTLUM V ZAVISLOSTI NA VLNOVE DELCE
    
%R=8.5e-3;                                            % [mm] definuj polomer ohybu
%lambda_start=1310e-9;                                % pocatecni vlnova delka
%lambda_stop=1330e-9;                                 % konecna vlnova delka
%step=5e-9;                                          % krok vypoctu
%precise=5e-5;                                        % presnost s jakou urcime veikost horni meze integralu (nizsi cislo = vyssi presnost)
lambda=lambda_start;                                 % promenna vlnove delky
%NA=0.118404;                                         % numericka apertura testovaneho vlakna                   

for(i=1:((lambda_stop-lambda_start)*10^9)/(step*10^9)+2) % upravuj 2 pokud chces rozsirit vlnovou delku

n1=1.4504;                                          % index lomu jadra
n2=sqrt(n1.^2-NA.^2);
%n2=1.4447;                                          %1.4446; n2<n1 !!! (indexy lomu pro namerene hodnoty)

if(n2<n1)                                           % podminka pro vedeni vlny
zeta=(1:precise:200).*1e4;                          % transformovana souradnice 'y' ve FT
pr=7;                                               % velikost ronbergovi tabulky (zvys, dokud se hodnoty posledni a predposledni nebudou rovnat)

[posledni, predposledni,loss,chi2,chi3,ch2,ch3,citatel]=ref7(R,lambda,zeta,pr,n2);
ch2s(1,i)=ch2;                                      % vektor data-pointeru pro reseni chi2
ch3s(1,i)=ch3;                                      % vektor data-pointeru pro reseni chi3
losss(1,i)=abs(loss);                               % vektor ztrat (pro jednotlive vlnove delky)
chi2s(i,:)=chi2;                                    % vektor chi2 (pro jednotlive vlnove delky)
chi3s(i,:)=chi3;                                    % vektor chi3 (pro jednotlive vlnove delky)
else
    fprintf('Spatne nastavene indexy lomu!')
    break                                           % prerus vypocet nepltati-li podminka n2<n1
end
lambdav(1,i)=lambda;                                % vektor vlnovych delek
lambda=lambda+step;                                 % nasledujici vlnova delka
end

figure(1) % zavislost utlumu na vlnove delce
plot(lambdav*10^9,(losss/1),'--r')
grid on
xlabel('\lambda [nm]')
ylabel('Loss [dB/m]')
title(['Zavislost utlumu na vlnove delce \lambda, R=' num2str(R*1e3),'mm', ', NA=',num2str(NA)])
%axis([1500 1600 0 200])
legend('Struktura s jednou vrstvou')
end

% ------------------------------------------------------------------------------------------------
% ------------------------------------------------------------------------------------------------
if(graf==2)                                         % UTLUM V ZAVISLOSTI NA POLOMERU OHYBU

%lambda=1310e-9;

%R_start=7.25e-3;
%R_stop=13e-3;
%step=0.1e-3;
%precise=5e-5;
R=R_start;
%NA=0.118404;                                         % numericka apertura testovaneho vlakna

for(i=1:((R_stop-R_start)*10^3)/(step*10^3)+2)      % upravuj 2 pokud chces rozsirit polomer ohybu

n1=1.4504;                                          % index lomu jadra
n2=sqrt(n1.^2-NA.^2);
%n2=1.4447;                                          % 1.4446; n2<n1 !!! (indexy lomu pro namerene hodnoty)

if(n2<n1)                                           % podminka pro vedeni vlny
zeta=(1:precise:200).*1e4;                          % transformovana souradnice 'y' ve FT
pr=7;                                               % velikost ronbergovi tabulky (zvys, dokud se hodnoty posledni a predposledni nebudou rovnat)

[~,~,loss,~,~,~,~,~]=ref7(R,lambda,zeta,pr,n2);
%ch2s(1,i)=ch2;                                      % vektor data-pointeru pro reseni chi2
%ch3s(1,i)=ch3;                                      % vektor data-pointeru pro reseni chi3
losss(1,i)=abs(loss);                               % vektor ztrat (pro jednotlive vlnove delky)
%chi2s(i,:)=chi2;                                    % vektor chi2 (pro jednotlive vlnove delky)
%chi3s(i,:)=chi3;                                    % vektor chi3 (pro jednotlive vlnove delky)
else
    fprintf('Spatne nastavene indexy lomu!')
    break                                           % prerus vypocet nepltati-li podminka n2<n1
end
Rv(1,i)=R;                                % vektor vlnovych delek
R=R+step;                                 % nasledujici vlnova delka
end

figure(1) % zavislost utlumu na polomeru ohybu
plot(Rv*10^3,losss,'--r')
grid on
xlabel('R [mm]')
ylabel('Ztraty ohybem [dB/m]')
title(['Zavislost utlumu na polomeru ohybu R'])
axis([R_start*1e3 R_stop*1e3 0 max(losss)])
legend(['Pro vlnovou delku \lambda=' num2str(lambda*1e9),'nm'])
end

% ------------------------------------------------------------------------------------------------
% ------------------------------------------------------------------------------------------------
if(graf==3)

%R=9e-3;
%lambda_start=1525e-9;
%lambda_stop=1565e-9;
%stepl=5e-9;
lambda=lambda_start;

%NA_start=0.03;
%NA_stop=0.140;
%step=0.001;
%precise=5e-5;
    

for(k=1:((lambda_stop-lambda_start).*1e9)/(stepl.*10^9)+1)
NA=NA_start;
for(i=1:((NA_stop-NA_start).*1e3)/(step.*10^3)+1)      % upravuj 1 pokud chces rozsirit numerickou aperturu

n1=1.4504;                                          % index lomu jadra
n2=sqrt(n1.^2-NA.^2);                               % n2<n1 !!! (indexy lomu pro namerene hodnoty)
n2v(1,i)=n2;
if(n2<n1)                                           % podminka pro vedeni vlny
zeta=(1:precise:200).*1e4;                          % transformovana souradnice 'y' ve FT
pr=7;                                               % velikost ronbergovi tabulky (zvys, dokud se hodnoty posledni a predposledni nebudou rovnat)

[~,~,loss,~,~,~,~,~]=ref7(R,lambda,zeta,pr,n2);
losss(k,i)=abs(loss);                               % vektor ztrat (pro jednotlive vlnove delky)
else
    fprintf('Spatne nastavene indexy lomu!')
    break                                           % prerus vypocet nepltati-li podminka n2<n1
end
NAv(k,i)=NA;                                % vektor vlnovych delek
NA=NA+step;                                 % nasledujici vlnova delka
end
lambdav(1,k)=lambda;
lambda=lambda+stepl;
end

close all
clc
for(i=1:k) 
figure(1) % zavislost utlumu na polomeru ohybu
hold on
col=jet(k);
plot(NAv(i,:),losss(i,:),'color',col(i,:),'LineStyle','--')
grid on
xlabel('NA [-]')
ylabel('Loss [dB/m]')
title(['Zavislost utlumu na numericke aperture NA, polomer ohybu R= ' num2str(R*1e3), ' mm'])
%axis([0.125 NA_stop 0 80])
%axis([NA_start NA_stop 0 150])
%legend(['Pro vlnovou delku \lambda= ' num2str(lambda*1e9),' nm'],'Location','SouthWest')
legend(['\lambda= ' num2str(lambdav(1,1)*1e9),' nm'], ['\lambda= ' num2str(lambdav(1,2)*1e9),' nm'], ['\lambda= ' num2str(lambdav(1,3)*1e9),' nm'],['\lambda= ' num2str(lambdav(1,4)*1e9),' nm'],['\lambda= ' num2str(lambdav(1,5)*1e9),' nm'],['\lambda= ' num2str(lambdav(1,6)*1e9),' nm'],['\lambda= ' num2str(lambdav(1,7)*1e9),' nm'],['\lambda= ' num2str(lambdav(1,8)*1e9),' nm'],['\lambda= ' num2str(lambdav(1,9)*1e9),' nm'],'Location','NorthEast')
end
end

% ------------------------------------------------------------------------------------------------
% ------------------------------------------------------------------------------------------------
if(graf==4)
    
%    NA_v=[0.09 0.11 0.13 0.14];                     % definuj pro jake numericke apertury NA chces vypocet provest
%    lambda=1600e-9;                                 % zadej pro jakou vlnovou delku chces vypocet provest


%R_start=7.0e-3;                                     % pocatecni hodnota polomeru ohybu
%R_stop=11e-3;                                       % konecna hodnota polomeru ohybu
%step=0.1e-3;                                        % krok vypoctu
%precise=5e-5;                                       % presnost s jakou je urcena horni mez integralu


for(j=1:length(NA_v))                               % cyklus pro vypocet zavislosti na NA
R=R_start;                                          % resetuje R na vychozi hodnotu po skonceni nasledujiciho cyklu
NA=NA_v(1,j);                                       % nastavi prislusnou numerickou aperturu definovanou vektorem vyse
for(i=1:((R_stop-R_start)*10^3)/(step*10^3)+2)      % upravuj 2 pokud chces rozsirit polomer ohybu, cyklus pro vypocet zvaislosti utlumu na polomeru ohybu

n1=1.4504;                                          % index lomu jadra
n2=sqrt(n1.^2-NA.^2);                               % 1.4446; n2<n1 !!! (indexy lomu pro namerene hodnoty)

if(n2<n1)                                           % podminka pro vedeni vlny
zeta=(1:precise:200).*1e4;                          % transformovana souradnice 'y' ve FT
pr=7;                                               % velikost ronbergovi tabulky (zvys, dokud se hodnoty posledni a predposledni nebudou rovnat)

[~,~,loss,~,~,~,~,~]=ref7(R,lambda,zeta,pr,n2);     % vypocet vysledneho utlumu ohybem pro prislusne parametry
losss(j,i)=abs(loss);                               % vektor ztrat (pro jednotlive vlnove delky)

else
    fprintf('Spatne nastavene indexy lomu!')
    break                                           % prerus vypocet nepltati-li podminka n2<n1
end
Rv(j,i)=R;                                          % vektor polomeru
R=R+step;                                           % nasledujici polomer

end
[maximum_o(j,:) pointer_o(j,:)]=max(losss(j,:));    % urceni maxima a jejiho ukazatele pro jednotlive faze simulace (tj. pro kazdou sadu hodnot prislusne NA)
R_o(j,:)=(Rv(j,pointer_o(j,:))).*1e3;               % urceni polomeru ohybu pro maximalni utlum
end

% vykresleni normovaneho grafu pro volbu graf==4
close all
check=0;                                            % pomocna kontrolni promenna zda je normovani provedeno korektne

%for(kk=1:length(NA_v))
%losss(kk,:)=losss(kk,:)./maximum_o(kk,:);           % normovani maximalni hodnotou
%[maximum_n(kk,:) pointer_n(kk,:)]=max(losss(kk,:)); % kontrola zda odpovidaji ukazatele, tedy zda bylo normovani provedeno pro spravne hodnoty
%end

for(kk=1:length(NA_v))
figure(1)                                           % zavislost utlumu na polomeru ohybu R, s parametrem NA
hold on
col=jet(length(NA_v));
plot(Rv(kk,:).*10^3,losss(kk,:),'color',col(kk,:),'LineStyle','--')
grid on
xlabel('R [mm]')
ylabel('Ztaty ohybem [ dB/m ]')
axis([7 10 0 1200])
legend(['NA= ' num2str(NA_v(1,1))],['NA= ' num2str(NA_v(1,2))],['NA= ' num2str(NA_v(1,3))],['NA= ' num2str(NA_v(1,4))],['NA= ' num2str(NA_v(1,5))],['NA= ' num2str(NA_v(1,6))],['NA= ' num2str(NA_v(1,7))])
title(['Zavislost utlumu na polomeru R, kde parametrem je numericka apertura NA, \lambda= ' num2str(lambda.*1e9), 'nm'])

%if(pointer_o==pointer_n)                            % odpovida-li hodnota stareho ukazatele (o-old) a noveho ukazatele(n-new)
%check=check+1;                                      % zvys hodnotu promenne check o 1
%else                                                % v opacnem pripade
%check=check-1;                                      % sniz o hodnotu 1
%end
%end
%if(check==length(NA_v))                             % jestlize vysledek scitani/odcitani v promenne check odpovida delce vektoru NA
%fprintf('Normovani je v poradku!')                  % je vse OK
%table=[R_o maximum_o];                              % a vytvor souhrnou tabulku polomeru ohybu a prislusnych maximanich utlumu
%else                                                % jinak
%fprintf('Chyba normovani!')                         % je vypocet spatny
%break                                               % a ukonci program
end
end

% ------------------------------------------------------------------------------------------------
% ------------------------------------------------------------------------------------------------
if(graf==5)
    
    
%lambda_v=[1310 1470 1550 1570 1625].*10^-9;       % zadej pro jakou vlnovou delku chces vypocet provest
%R_start=7e-3;                                     % pocatecni hodnota polomeru ohybu
%R_stop=10e-3;                                       % konecna hodnota polomeru ohybu
%step=0.05e-3;                                        % krok vypoctu
%precise=5e-5;                                       % presnost s jakou je urcena horni mez integralu
%NA=0.115804;                                         % numericka apertura testovaneho vlakna

for(j=1:length(lambda_v))                               % cyklus pro vypocet zavislosti na NA
R=R_start;                                          % resetuje R na vychozi hodnotu po skonceni nasledujiciho cyklu
lambda=lambda_v(1,j);                                       % nastavi prislusnou numerickou aperturu definovanou vektorem vyse
for(i=1:((R_stop-R_start)*10^3)/(step*10^3)+2)      % upravuj 2 pokud chces rozsirit polomer ohybu, cyklus pro vypocet zvaislosti utlumu na polomeru ohybu

n1=1.4504;                                          % index lomu jadra
n2=sqrt(n1.^2-NA.^2);
%n2=1.4447;                                          %1.4446; n2<n1 !!! (indexy lomu pro namerene hodnoty)

if(n2<n1)                                           % podminka pro vedeni vlny
zeta=(1:precise:200).*1e4;                          % transformovana souradnice 'y' ve FT
pr=7;                                               % velikost ronbergovi tabulky (zvys, dokud se hodnoty posledni a predposledni nebudou rovnat)

[~,~,loss,~,~,~,~,~]=ref7(R,lambda,zeta,pr,n2);     % vypocet vysledneho utlumu ohybem pro prislusne parametry
losss(j,i)=abs(loss);                               % vektor ztrat (pro jednotlive vlnove delky)

else
    fprintf('Spatne nastavene indexy lomu!')
    break                                           % prerus vypocet nepltati-li podminka n2<n1
end
Rv(j,i)=R;                                          % vektor polomeru
R=R+step;                                           % nasledujici polomer

end
[maximum_o(j,:) pointer_o(j,:)]=max(losss(j,:));    % urceni maxima a jejiho ukazatele pro jednotlive faze simulace (tj. pro kazdou sadu hodnot prislusne NA)
R_o(j,:)=(Rv(j,pointer_o(j,:))).*1e3;               % urceni polomeru ohybu pro maximalni utlum
end

% vykresleni normovaneho grafu pro volbu graf==5
close all
check=0;                                            % pomocna kontrolni promenna zda je normovani provedeno korektne

%for(kk=1:length(lambda_v))
%losss(kk,:)=losss(kk,:)./maximum_o(kk,:);           % normovani maximalni hodnotou
%[maximum_n(kk,:) pointer_n(kk,:)]=max(losss(kk,:)); % kontrola zda odpovidaji ukazatele, tedy zda bylo normovani provedeno pro spravne hodnoty
%end

for(kk=1:length(lambda_v))
figure(1)                                           % zavislost utlumu na polomeru ohybu R, s parametrem lambda
hold on
col=jet(length(lambda_v));
plot(Rv(kk,:).*10^3,losss(kk,:),'color',col(kk,:),'LineStyle','--')
grid on
xlabel('R [mm]')
ylabel('Loss [ dB/m ]')
legend(['\lambda= ' num2str(lambda_v(1,1)*10^9),' nm'],['\lambda= ' num2str(lambda_v(1,2)*10^9),' nm'],['\lambda= ' num2str(lambda_v(1,3)*10^9),' nm'],['\lambda= ' num2str(lambda_v(1,4)*10^9),' nm'],['\lambda= ' num2str(lambda_v(1,5)*10^9),' nm'])
title(['Zavislost utlumu na polomeru R, kde parametrem je vlnova delka \lambda, NA= 0.128'])
end
end

%% dodatecne simulace
figure(2) % zavislost chi2 na souradnici zeta, kde parametrem je vlnova delka
plot(zeta,chi2s)
hold on
plot(zeta(:,ch2s),ch2s,'ro')
grid on
xlabel('\zeta [-]')
ylabel('\chi_2 [-]')
title('Zavislost \chi_2 na souradnici \zeta')
%axis([0 10e5 -5e11 5e+11])
%legend(['\lambda=' num2str(lambdav(1,1)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,2)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,3)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,4)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,5)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,6)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,7)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,8)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,9)*1e9),'nm']);

%%
figure(3) % zavislost chi2 na souradnici zeta, kde parametrem je vlnova delka (ZOOM)
plot(zeta,chi2s)
hold on
plot(zeta(:,ch2s),ch2s,'ro')
grid on
xlabel('\zeta [-]')
ylabel('\chi_2 [-]')
title('Zavislost \chi_2 na souradnici \zeta')
%axis([5.47e5 5.65e5 -0.1e11 0.15e+11])
%legend(['\lambda=' num2str(lambdav(1,1)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,2)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,3)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,4)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,5)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,6)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,7)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,8)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,9)*1e9),'nm']);

%%
figure(4) % zavislost chi3 na souradnici zeta, kde parametrem je vlnova delka
plot(zeta,chi3s)
hold on
plot(zeta(:,ch3s),ch3s,'go')
grid on
xlabel('\zeta [-]')
ylabel('\chi_3 [-]')
title('Zavislost \chi_3 na souradnici \zeta')
%legend(['\lambda=' num2str(lambdav(1,1)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,2)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,3)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,4)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,5)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,6)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,7)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,8)*1e9),'nm'], ['\lambda=' num2str(lambdav(1,9)*1e9),'nm']);
%%
close all
citatel_real=real(citatel);
citatel_imag=imag(citatel);
figure(5)
plot(zeta, citatel_real,'-b')
hold on;
plot(zeta, citatel_imag,'--r')
hold on;
plot(zeta(1,ch2),0,'ro')
hold on;
plot(zeta(1,ch3),0,'go')
grid on;
xlabel('\zeta [-]')
ylabel('sqrt(\chi_2) . sqrt(\chi_3) [-]')
title('Zavislost citatele 3.integrandu na souradnici \zeta')
legend('Ralna cast citatele','Imaginarni cast citatele','horni mez pro \chi_2', 'horni mez pro \chi_3','Location','SouthWest')
