%% Scintillation index
function[] = testScintillation(Lo,lo,ro,d,lambda,strength_of_turbulence,beamkind,filename)
%clear all;clc;close all;

%% Einstellungen
N = 2^6;  
L = 700; 
Fsx = N/L;
omega_sx = 2*pi*Fsx;
delta_kx = omega_sx/N; %increment in spatial frequencies
k_x = (-omega_sx/2:delta_kx:omega_sx/2-omega_sx/N)';

switch strength_of_turbulence
       case 1
            Cn_Quadrat = 10^-17;
       case 2
            Cn_Quadrat = 10^-16;
       case 3 
            Cn_Quadrat = 10^-15;
end
Lo = Lo*10^3;%change Lo from m in mm
km = 5.92/lo;
ko = 2*pi/Lo;
lambda = lambda*10^-6;%change lambda from nm in mm
k0 = 2*pi/lambda; 
d0 = 0;
distance = d;
d = d*10^6;%change d from km in mm
delta_z = 1000;%[mm]
N_int = ceil(d/delta_z);
alpha = lambda/(pi*ro^2);
epsilon_zero = 1;
L = ro;
x = (-L/2:L/N:L/2-L/N)';
k_y = k_x;
[Kx,Ky] = meshgrid(k_x,k_y);
k = sqrt(Kx.^2 + Ky.^2); %spatial frequency
A =  -(k.^2/(1i*2*k0));
Phi_2D = 0.033*Cn_Quadrat.*exp(-(k.^2)./km^2)./(k.^2+ko^2).^(11/6);
y = x;
[X,Y] = meshgrid(x,y);
delta_ky = delta_kx;
D = exp(A.*delta_z.*0.5);
%%
% -------Impulserzeugung --------
 switch beamkind
        case 0 %gauss
     g_start_2D = epsilon_zero.*exp(-(X.^2+Y.^2)./ro^2);
        case 1 % bessel
     kr = 2.405/ro;%radial wave vector
     kz = sqrt(k0^2-kr^2);
     r = sqrt(X.^2+Y.^2);%radial coordinate
     g_start_2D = epsilon_zero.*besselj(0,kr.*r);
 end
   
   gx = g_start_2D;
   averageNumber = 100;
   m = 2;
   p = length(gx)/2;
   q = 2;
   SI_array(1,1) = 0;

  for s = d0:delta_z:d
    I = zeros(size(gx));
    I_square = zeros(size(gx));
    colored_Noise_2D = ifft2(ifftshift(fftshift(fft2(randn(N))).*sqrt(Phi_2D*delta_kx*delta_ky)));
    delta_n_2D = real(colored_Noise_2D);
    S = delta_n_2D*k0;
    R = exp(-1i.*S.*delta_z);
    gx1 = ifft2(ifftshift((fftshift(fft2((ifft2(ifftshift(fftshift(fft2(gx)).*D))).*R))).*D));    

   for val = 1:averageNumber
        I = I + (abs(gx1)).^2;%Irradiance or Intensity[W/m^2]
        I_square = I_square + (abs(gx1)).^4;%I^2
   end
        I_interior = I(p-q:p+q, p-q:p+q);
        I_square_interior = I_square(p-q:p+q, p-q:p+q);
        SI = (mean(mean(I_square_interior/averageNumber))-(mean(mean(I_interior/averageNumber)))^2)/((mean(mean(I_interior/averageNumber)))^2);%scintillation index
        SI_array(1,m) = SI;  
        m = m+1
        gx = gx1;
  end  

        df = 1/length(SI_array);
        f = distance*(0:df:1-df);
        set(0, "defaultaxesfontsize", 9);
        set(0, "defaultaxesfontname", "Helvetica");
        set(0, "defaulttextfontsize", 9);
        set(0, "defaulttextfontname", "Helvetica");
        set(0, "defaultaxeslinewidth", 1);
        plot(f,SI_array,'LineWidth',3);
        xlabel('propagation distance z [km]');
        ylabel('\sigma^2');
        %title('scintillation index');          
        print(filename,'-dpng',"-S930,350");
end
