function[] = Turbulence(Lo,lo,ro,d,lambda,strength_of_turbulence,presentation_type,analyse,beamkind,mode,filename)
%Modeling atmpheric turbulences (scintillations) in propagation of free-space optical signals
%close all;clc;clear all;

%parameters
switch strength_of_turbulence
 case 0
      Cn_Quadrat = 0;
 case 1
      Cn_Quadrat = 10^-17;
 case 2
      Cn_Quadrat = 10^-16;
 case 3
      Cn_Quadrat = 10^-15;%refractive index structure constant
end

N = 2^7; 
L = 300; 
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)';
%Lo = 2000; % outer scale[mm]
Lo = Lo*10^3;%change Lo from m in mm 
%lo =7; % inner scale[mm]
km = 5.92/lo;
ko = 2*pi/Lo;
Phi_1D = 0.033*Cn_Quadrat.*exp(-(k_x.^2)./km^2)./(k_x.^2+ko^2).^(11/6);
a = max(Phi_1D);
b = k_x(Phi_1D >= a*10^-6);
omega_sx = 2*b(end);

N=2^7;
L=7*ro;
Fsx = N/L;
omega_sx = 2*pi*Fsx;
delta_kx = omega_sx/N;
k_x = (-omega_sx/2:delta_kx:omega_sx/2-delta_kx);


%Phi_1D = 0.033*Cn_Quadrat.*exp(-(k_x.^2)./km^2)./(k_x.^2+ko^2).^(11/6);

%ro = 30; % beam radius[mm] 
%lambda = 0.55e-3; % wavelength[mm]
lambda = lambda*10^-6;%change lambda from nm in mm
k0 = 2*pi/lambda; %wavenumber
%d = 1000000;                         % propagation distance in mm
distance = d;
d = d*10^6; %change propagation distance  from km in mm
delta_z = 500;%in mm
N_int = ceil(d/delta_z);
%alpha = lambda/(pi*ro^2);
epsilon_zero = 1;
x = linspace(-L/2,L/2,N)';

%%
% ---------- Impulserzeugung -----------------

%g_start = epsilon_zero.*exp(-(x/ro).^2);       % gaußian beam
%gx = g_start;
%E = zeros(N_int, N);
%E(1,:) = abs(gx);
%m = 2;
%%
%Implementierung von E(x,z+delta_z) >>SSFM
      
%    for z = linspace(0,d,N_int)
 %       white_noise_1D = randn(N,1);
        % n = 1000; autocorrelation_1D     = zeros(n,2*N-1);
        % Initialisierung von autokorrelation
        %for i = 1:n
        % generate colored_Noise/phasescreen
  %      colored_Noise_1D = ifft(ifftshift(fftshift(fft(white_noise_1D)).*sqrt(Phi_1D*delta_kx)));
        %  autocorrelation_1D(i,:) = xcorr(colored_Noise_1D,'biased');
        % end
        %delta_n = zeros(1,N);
   %     delta_n = real(colored_Noise_1D);
    %    S    = delta_n.*k0;
     %   Gx = fftshift(fft(gx)).*exp(1i.*(k_x.^2).*delta_z/(4*k0));
      %  gxp = (ifft(ifftshift(Gx))).*exp(-1i.*S.*delta_z);
       % Gxp = (fftshift(fft(gxp))).*exp(1i.*(k_x.^2).*delta_z/(4*k0));
       % gx = ifft(ifftshift(Gxp));
       % E(m,:) = abs(gx);
       % m = m+1;
   % end

%Ryy_1D = mean(autocorrelation_1D); %>>Mittelwert von xcorr(colored_Noise)
%hold on;
%psd_1D = abs(fftshift(fft(Ryy_1D))); 
%plot(linspace(-omega_sx/2,omega_sx/2,2*N-1),
%psd_1D,'r');

%hold on; 
%Ryy_1D_glatt = hamming(size(autocorrelation_1D,2)).'.*mean(autocorrelation_1D); 
%psd_1D_glatt = abs(fftshift(fft(Ryy_1D_glatt))); 
%plot(linspace(-omega_sx/2,omega_sx/2,2*N-1),
%psd_1D_glatt,'g'); 
%xlabel("k");
% legend("\Phi(k)", "Ryy(k)","Ryy_f(k)");
%grid on;


%pulse_evolution = surf(E);
%set(pulse_evolution,'LineStyle','none')
%shading interp;
%colormap jet;
%grid on;
%title(('evolution of 1D-Pulse according to the SSM');
%xlabel("x");
%ylabel("y");
%axis tight;
%m = 1;

     %E_zero_1D = zeros(N_int, N);
       % for z = linspace(0,d,N_int)
        %    E_zero_1D(m,:) = abs((epsilon_zero/(1+1i*alpha*z)).*exp(1i*k0*z-0.5*k0*alpha.*x.^2./(1+1i*alpha*z)));
         %   m = m+1;
       % end
 
%%
% 2 dimensionale Betrachtung
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));
%[Kx,Ky,Kz]= meshgrid(k_x,k_y ,k_z); 
%k = sqrt(Kx.^2+Ky.^2 + Kz.^2);%spatial frequency for 3D colored noise

Phi_2D = 0.033*Cn_Quadrat.*exp(-(k.^2)./km^2)./(k.^2+ko^2).^(11/6);

%figure(4);
%h = surf(Kx,Ky,Phi_2D);
%set(h,'LineStyle','none');
%shading interp;
%colormap jet;
%zlim auto;
%grid on;
%xlabel("kx");
%ylabel("ky");
%title("\Phi(k)");
%axis tight;

%Ryy_2D = squeeze(mean(correlation_2D));
%[p,q] = meshgrid(linspace(-omega_sx/2,omega_sx/2,2*length(Phi_1D)-1));

%figure(5); 
%psd_2D    = abs(fftshift(fft2(Ryy_2D))); % psd_2D == Phi_2D ? Wiener-Kintchine Theorem 
%b = surf(p, q, psd_2D);
%set(b,'LineStyle','none'); 
%shading interp;
%colormap jet; 
%grid on; 
%title('Ryy(k)');
% xlabel("kx");
%ylabel("ky");

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
           r = sqrt(X.^2+Y.^2);%radial coordinate
           g_start_2D = epsilon_zero.*besselj(0,kr.*r);
          end

gx = g_start_2D;

%%
%Implementierung von E(x,y,z+delta_z) >> SSFM

    data_memory(:,:,1) = abs(g_start_2D);
    position = zeros(1,10); 
    m = 2;
    n = 1;
    
  for z = 0:delta_z:d    
        %whiteNoise_2D = randn(N);
        %fourier_whiteNoise = fftshift(fft2(whiteNoise_2D));
        %fourier_coloredNoise = ifftshift(fourier_whiteNoise.*sqrt(Phi_2D*delta_kx*delta_ky));
        %colored_Noise_2D = ifft2(fourier_coloredNoise);
        colored_Noise_2D = ifft2(ifftshift(fftshift(fft2(randn(N))).*sqrt(Phi_2D)*Fsx^2));%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));
        if mod(n,ceil(N_int/10)) == 0 
        data_memory(:,:,m) = abs(gx1)./max(max(abs(gx1)));
        position(1,m) = z;
        m = m+1;
        end
        n = n+1;
        gx = gx1;
  end
     	
   
  switch analyse
     case 'evaluation'

	switch mode
		case 0
      
        g = figure(1);
        set(g,'visible','off');       
        set(0,'defaultaxesfontsize',6);
        set(0,'defaultaxesfontname',"Helvetica");
        set(0,'defaulttextfontsize',6);
        set(0,'defaulttextfontname',"Helvetica");
           
        xq = linspace(-L/2,L/2,2^9);
        yq = xq;
        [Xq,Yq] = meshgrid(xq,yq);
        V5 = interp2(X,Y,abs(g_start_2D),Xq,Yq);
        imagesc(xq,yq,V5);
        if beamkind == 0
        title('Gaussian beam','fontweight','bold');
        else
        title('Bessel beam','fontweight','bold');
        end
        shading interp;
        colormap('jet');
        xlabel("x [mm]");
        ylabel("y [mm]");
        caxis([0,1]);
        cb1 = colorbar;
        hL = ylabel(cb1,'^{E}/_{E_0}');
        set(hL,'Rotation',0);
        shading interp;
        colormap jet;
        print(filename,'-dpng', "-S660,490");


		case 1

        g = figure(1);
        set(g,'visible','off');
        set(0,'defaultaxesfontsize',6);
        set(0,'defaultaxesfontname',"Helvetica");
        set(0,'defaulttextfontsize',6);
        set(0,'defaulttextfontname',"Helvetica");


        xq = linspace(-L/2,L/2,2^9);
        yq = xq;
        [Xq,Yq] = meshgrid(xq,yq);
        V5 = interp2(X,Y,abs(g_start_2D),Xq,Yq);
        if beamkind == 0
        title('Gaussian beam','fontweight','bold');
        else
        title('Bessel beam','fontweight','bold');
        end
        surf(xq,yq,V5,'EdgeColor','none');
        xlabel("x [mm]");
        ylabel("y [mm]");
        caxis([0,1]);
        cb2 = colorbar;
        hL = ylabel(cb2,'^{E}/_{E_0}');
        set(hL,'Rotation',0);
        axis tight;
        shading interp;
        colormap jet;
        print(filename,'-dpng',"-S660,490");



		case 2


        g = figure(1);
        set(g,'visible','off');
        set(0,'defaultaxesfontsize',6);
        set(0,'defaultaxesfontname',"Helvetica");
        set(0,'defaulttextfontsize',6);
        set(0,'defaulttextfontname',"Helvetica");

        subplot(2,1,1);
        xq = linspace(-L/2,L/2,2^9);
        yq = xq;
        [Xq,Yq] = meshgrid(xq,yq);
        V5 = interp2(X,Y,abs(g_start_2D),Xq,Yq);
        imagesc(xq,yq,V5);
        if beamkind == 0
        title('Gaussian beam','fontweight','bold');
        else
        title('Bessel beam','fontweight','bold');
        end
        shading interp;
        colormap('jet');
        xlabel("x [mm]");
        ylabel("y [mm]");
        caxis([0,1]);
        cb1 = colorbar;
        hL = ylabel(cb1,'^{E}/_{E_0}');
        set(hL,'Rotation',0);
        shading interp;
        colormap jet;
        print(filename,'-dpng', "-S545,545");
        subplot(2,1,2);
        surf(xq,yq,V5,'EdgeColor','none');
        xlabel("x [mm]");
        ylabel("y [mm]");
        caxis([0,1]);
        cb2 = colorbar;
        hL = ylabel(cb2,'^{E}/_{E_0}');
        set(hL,'Rotation',0);
        axis tight;
        shading interp;
        colormap jet;
        print(filename,'-dpng',"-S545,545");



	end

                  
     case 'simulation'
       switch presentation_type

         case 'dynamic'
                       set(0,"defaultaxesfontsize",8);
                       set(0,"defaultaxesfontname","Helvetica");
                       set(0,"defaulttextfontsize",8);
                       set(0,"defaulttextfontname","Helvetica");
                       xq = linspace(-L/2,L/2,2^9);
                       yq = xq;
                       [Xq,Yq] = meshgrid(xq,yq);                     
                    
                       for l = 1:length(position)
			switch mode
				case 0                  

                       beamMagnitude = data_memory(:,:,l);
                       V1 = interp2(X,Y,beamMagnitude,Xq,Yq);
                       h0 = imagesc(xq,yq,V1);
                       title(['z = ',num2str(position(l)/10^6),'km']);
                       xlabel("x [mm]");
                       ylabel("y [mm]");
                       caxis([0,1]);
                       cb1 = colorbar;
                       hL = ylabel(cb1,'^{E}/_{E_0}');
                       set(hL,'Rotation',0);
                       axis tight;
                       shading interp;
                       colormap jet;
                       new_filename = [filename, '_' num2str((l-1))];
                       print(new_filename, '-dpng',"-S660,490");  
                       delete(h0);
                      
				case 1

                       beamMagnitude = data_memory(:,:,l);
	               V2 = interp2(X,Y,beamMagnitude.^2,Xq,Yq);
                       h1 = surf(xq,yq,V2,'EdgeColor','none');
                       xlabel("x [mm]");
                       ylabel("y [mm]");
                       caxis([0,1]);
                       cb2 = colorbar;
                       hL = ylabel(cb2,'^{I}/_{I_0}');
                       set(hL,'Rotation',0);
                       axis tight;
                       shading interp;
                       colormap jet;
                       hold on;
                       h2 = plot3(xq,L/2*ones(size(yq)),V2(:,end/2).^2);
                       h3 = plot3(L/2*ones(size(xq)),yq,V2(end/2,:).^2);
                       set(gca, 'children', flipud(get(gca,'children')));
                       zlim([0,1]);
                       hold off;
                       new_filename = [filename, '_' num2str((l-1))];
                       print(new_filename, '-dpng',"-S660,490");
                       delete(h1);
                       delete(h2);
                       delete(h3);


				case 2

                       subplot(2,1,1);
                       beamMagnitude = data_memory(:,:,l);
                       V1 = interp2(X,Y,beamMagnitude,Xq,Yq);
                       h0 = imagesc(xq,yq,V1);
                       title(['z = ',num2str(position(l)/10^6),'km']);
                       xlabel("x [mm]");
                       ylabel("y [mm]");
                       caxis([0,1]);
                       cb1 = colorbar;
                       hL = ylabel(cb1,'^{E}/_{E_0}');
                       set(hL,'Rotation',0);
                       axis tight;
                       shading interp;
                       colormap jet;

                       subplot(2,1,2);
                       V2 = interp2(X,Y,beamMagnitude.^2,Xq,Yq);
                       h1 = surf(xq,yq,V2,'EdgeColor','none');
                       xlabel("x [mm]");
                       ylabel("y [mm]");
                       caxis([0,1]);
                       cb2 = colorbar;
                       hL = ylabel(cb2,'^{I}/_{I_0}');
                       set(hL,'Rotation',0);
                       axis tight;
                       shading interp;
                       colormap jet;
                       hold on;
                       h2 = plot3(xq,L/2*ones(size(yq)),V2(:,end/2).^2);
                       h3 = plot3(L/2*ones(size(xq)),yq,V2(end/2,:).^2);
                       set(gca, 'children', flipud(get(gca,'children')));
                       zlim([0,1]);
                       hold off;
                       new_filename = [filename, '_' num2str((l-1))];
                       print(new_filename, '-dpng',"-S700,700");
                       delete(h0);
                       delete(h1);
                       delete(h2);
                       delete(h3);


			end
                       end                               
                       
                       
        case 'static'
            if mode == 0
     g = figure(1);
     set(g,'visible','off');
     set(0, 'defaultaxesfontsize', 6);
     set(0, 'defaultaxesfontname', "Helvetica");
     set(0, 'defaulttextfontsize', 10);
     set(0, 'defaulttextfontname', "Helvetica");
  

     subplot(2,1,1);
     xq = linspace(-L/2,L/2,2^9);
     yq = xq;
     [Xq,Yq] = meshgrid(xq,yq);
     V3 = interp2(X,Y,abs(g_start_2D),Xq,Yq);
     imagesc(xq,yq,V3);
     title('source plane','fontweight','bold');
     shading interp;
     colormap('jet');
     xlabel("x [mm]");
     ylabel("y [mm]");
     caxis([0,1]);
     h1 = colorbar;
     hL = ylabel(h1,'^{E}/_{E_0}');
     set(hL,'Rotation',0);     
     shading interp;
     colormap jet;
     print(filename, '-dpng', "-S510,510");
 
     subplot(2,1,2);
     V4 = interp2(X,Y,abs(gx)/max(max(abs(gx))), Xq,Yq);
     imagesc(xq,yq,V4);
     title('observation plane','fontweight','bold');
     xlabel("x [mm]");
     ylabel("y [mm]");
     caxis([0,1]);
     h2 = colorbar;
     hL = ylabel(h2,'^{E}/_{E_0}');
     set(hL,'Rotation',0);
     shading interp;
     colormap jet;
     print(filename, '-dpng',"-S510,510");
           
            else
   g=figure(1);
     set(g,'visible','off');
     set(0, 'defaultaxesfontsize', 8);
     set(0, 'defaultaxesfontname', "Helvetica");
     set(0, 'defaulttextfontsize', 8);
     set(0, 'defaulttextfontname', "Helvetica");
     
     xq = linspace(-L/2,L/2,2^9);
     yq = xq;
     [Xq,Yq] = meshgrid(xq,yq);  
     V4 = interp2(X,Y,abs(gx)/max(max(abs(gx))), Xq,Yq);
     imagesc(xq,yq,V4);
     title('observation plane','fontweight','bold');
     xlabel("x [mm]");
     ylabel("y [mm]");
     caxis([0,1]);
     h2 = colorbar;
     hL = ylabel(h2,'^{E}/_{E_0}');
     set(hL,'Rotation',0);
     shading interp;
     colormap jet;
     print(filename, '-dpng',"-S660,490");
            end    
        end
  end

