clc;
clear all
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓ Import data and Values of variables ↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓%

LANDA=fliplr(importdata('C:\Users\msi\Desktop\data.txt')');

Landa=1e-9*LANDA(1,:);
phase_value=LANDA(2,:);

length_Landa=length(Landa);
clear LANDA 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓ Main ↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓%

landa_bragg=1550e-9;

neff=1.45;

L=2e-3;   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓ Apodization ↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓%

S=0.5;

AC=0.0001;

DC=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓ Chirp ↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓%

Delta=0.0001e-6; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓ Calculation ↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓%

seg=101;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Lambda_B=landa_bragg/(2*neff);

z=L :-L/seg: 0;

length_z=length(z);

for j=1:length_z-1

Z(j)=(z(j)+z(j+1))/2;

end

length_Z=length(Z);
clear z length_z
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Apodization= {['Uniform(without Apodization)'],['Gaussian'],['Hyperbolic tangent']};
type = menu('Apodization type? :','Uniform(without Apodization)','Gaussian','Hyperbolic tangent');
Chirp= {['Uniform(without chirp)'],['Linear'],['Quadratic'],['Square root'],['Cubic root']};
Type = menu('Chirp type? :','Uniform(without chirp)','Linear','Quadratic','Square root','Cubic root');
%↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓ Apodization functions ↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓%

if type==1

       A(1:length_Z)=1;                                                                                   % Uniform 
      dA(1:length_Z)=0;                                                                                   % Derivative 
% figure(1)
% plot(Z,A)
% title('Uniform Apodization Function')
% xlabel('Length(mm)')
% ylabel('Value')

end

%---------------------------------------------------------------------------------------------------------------------%

if type==2
   for j=1:length_Z
      A(j)=exp((-log(2))*(2*(Z(j)-0.5*L)/(S*L))^2);                                                       % Gaussian 
     dA(j)=(-8*log(2)*(Z(j)-0.5*L)/(L*S)^2)*exp((-log(2))*(2*(Z(j)-0.5*L)/(S*L))^2);                      % Derivative
   end
% figure(2)
% plot(Z,A)
% title('Gaussian Apodization Function')
% xlabel('Length(mm)')
% ylabel('Value')

end

%---------------------------------------------------------------------------------------------------------------------%

if type==3
    for j=1:length_Z
        A(j)=(tanh((S*Z(j))/L))*tanh(S*(1-((Z(j))/L)))+1-(tanh(0.5*S))^2;                                 % Hyperbolic tangent 
       dA(j)=(S/L)*(((sech(S*Z(j)/L))^2)*tanh(S*(1-Z(j)/L))-(tanh(S*Z(j)/L))*(sech(S*(1-Z(j)/L)))^2) ;    % Derivative
    end
% 
% figure(3)
% plot(Z,A)
% title('Hyperbolic tangent Apodization function')
% xlabel('Length(mm)')
% ylabel('Value')

end

clear Apodization type
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓ Chirp functions ↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓%

if Type==1
   Lambda(1:length_Z)=Lambda_B;                                                                           % Uniform
  dLambda(1:length_Z)=0;                                                                                  % Derivative
end

%---------------------------------------------------------------------------------------------------------------------%

if Type==2
    for j=1:length_Z
       Lambda(j)=Lambda_B-Delta*((Z(j)-0.5*L)/L);                                                         % Linear
      dLambda(j)=-Delta/L;                                                                                % Derivative
    end
end

%---------------------------------------------------------------------------------------------------------------------%

if Type==3
    for j=1:length_Z
       Lambda(j)=Lambda_B-Delata*((Z(j)/L)^2-0.25);                                                       % Quadratic
      dLambda(j)=-2*Z(j)*Delta/L^2;                                                                       % Derivative
    end
end

%---------------------------------------------------------------------------------------------------------------------%

if Type==4
    for j=1:length_Z
       Lambda(j)=Lambda_B-Delta*(sqrt(Z(j)/L)-sqrt(0.5));                                                 % Square root
      dLambda(j)=-Delta/(2*sqrt(Z(j)*L));                                                                 % Derivative
    end
end

%---------------------------------------------------------------------------------------------------------------------%

if Type==5
    for j=1:length_Z
       Lambda(j)=Lambda_B-Delta*((Z(j)/L)^(1/3)-0.5^(1/3));                                               % Cubic root
      dLambda(j)=-Delta/(3*(Z(j)^(2/3))*L^(1/3));                                                         % Derivative
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓ K=1/2 dpi/dz ↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓%

for p=1:length_Landa
    for j=1:length_Z
      K(p,j)=(-4*pi*(neff+DC*A(j)+AC*A(j))*Z(j)/landa_bragg^2)*2*(((DC*dA(j)+AC*dA(j))*Lambda(j)+neff+DC*A(j)+AC*A(j))*dLambda(j));  %
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓ Delta & Sigma & Sigma-hat & kapa & Gamma & Transfer matrix ↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓%
step=L/(seg);

for p=1:length_Landa

    Tt=[1 0;0 1];

    for j=1:length_Z
    delta(p,j)=2*pi*neff*(((1)/(Landa(p)))-((1)/(landa_bragg)));  
    sigma(p,j)=(2*pi/Landa(p))*(DC*A(j));  
    sigma_hat(p,j)=delta(p,j)+sigma(p,j)+K(p,j); 
    kapa(p,j)=(pi/Landa(p))*AC*A(j); 
    Gamma(p,j)=sqrt(((kapa(p,j))^2)-((sigma_hat(p,j))^2));

    t11(p,j)=cosh(Gamma(p,j)*step)-sqrt(-1)*(sigma_hat(p,j)/Gamma(p,j))*sinh(Gamma(p,j)*step);

    t12(p,j)=-sqrt(-1)*(kapa(p,j)/Gamma(p,j))*sinh(Gamma(p,j)*step);

    t21(p,j)=sqrt(-1)*(kapa(p,j)/Gamma(p,j))*sinh(Gamma(p,j)*step);

    t22(p,j)=cosh(Gamma(p,j)*step)+sqrt(-1)*(sigma_hat(p,j)/Gamma(p,j))*sinh(Gamma(p,j)*step);

    T1(:,:,j) =[t11(p,j) t12(p,j);t21(p,j) t22(p,j)];
    T=T1;
    Tt=Tt*T(:,:,j) ;
    end

    T=Tt*[1;0];

    TT=T(2,1)/T(1,1);

    a(p)=-phase(TT);

end


plot(Landa,a)


compare=a-phase_value;