if 0
pmd             = 0.2 * ( 1e-12 / sqrt(1e3) ) ;   % 0.1 ps/sqrt(km) -> 1e-9 -> s/sqrt(m)
L               = 1000;     % m
corr_len        = 100;
num_wave_plates = 100;
wp_len          = L / num_wave_plates; %waveplate length
dgd             = pmd * sqrt(L); 
fdac            = 120e9;    % GHz

fsim            = fdac * 512 ; %oversampled simulation frequency
dt              = 1/fsim;   % sample time
nt              = 4;      % sampled signal length
omega           = 2*pi*[(0:nt/2-1),(-nt/2:-1)]/(dt*nt) ;  %angular frequency vector for optical signal of length nt
omega           = 2*pi*fsim ;  %angular frequency vector for optical signal of length nt

for rlz = 1:20000

    db0 = [];
    db1 = [];
    delta_beta = [];
    
    db0             = (rand(num_wave_plates,1)*2*pi - pi); % normal distr. between -pi <-> +pi

    db1             = sqrt(3*pi/8)*(dgd/fsim)/ num_wave_plates .* omega; % linear increasing delta beta 1
    
    %loop over waveplates
    for n_wp = 1:length(db0)
        delta_beta(n_wp,:)      = (db1+db0(n_wp))./corr_len;
    end
    
    dphi            = delta_beta*wp_len;       %phase shift due to propagation constant = beta * L

    if rlz == 1
        figure()
        plot(cumsum( delta_beta ));
    end
    
    dphi_end(rlz)   = sum( dphi );

end

figure;
histogram((dphi_end),100);

deltaT = mean(abs(dphi_end));

%H               = exp(-1j*dphi);       % Filter to apply phase shift in freq. domain
end


PMD         = 0.1 * ( 1e-12 / sqrt(1e3) ); %PMD s/sqrt(m)
L           = 10e3;                   %m
sigma_dgd   = PMD * sqrt(L);       % in sec.
disp(['variance of DGD in pico seconds: ',num2str(sigma_dgd*1e12)]);
disp(['variance of DGD in pico seconds: ',num2str(sigma_dgd*1e12)]);

dgd = [];
for i = 1:4
    dgd(i,:)       = DGD(sigma_dgd,10000);
end

rms =  PMD * sqrt(40e3) * 1e12 ;       % in sec.
t = [0:0.01:100];
z2 = sqrt(2/pi).*(t.^2)/(rms.^3).*exp(-t.^2/(2*rms^2)) ;
figure;
histogram(dgd*1e12,1000,"Normalization","pdf");
hold on
plot(t(1:200),z2(1:200),'r');
hold off

function z = DGD(sigma_dgd,N)

    %PMD = sqrt( (PMD*1e12)^2/3 );

   % PMD = PMD * 1/( 1e-12 / sqrt(1e3) );

    sigma_dgd = sigma_dgd*1e12;

    x = [0:N];

    almaga = randn(3,N)*sigma_dgd;

    y = almaga.^2;

    z0 = sqrt(sum(y,1));

   
    disp(mean(z0));
    disp(std(z0));
  % disp(var(z0));

    z1 = [];

    for l = 0:0.1:N
        z1 = [z1 sum(z0(1,:)>l & z0(1,:)<l+0.1) / N ];
    end

     z = z0 *1e-12;

    rms = sigma_dgd;
    t = [0:0.01:N];
    z2 = sqrt(2/pi).*(t.^2)/(rms.^3).*exp(-t.^2/(2*rms^2)) ;

    figure;
    histogram(z0,1000,"Normalization","pdf");
    hold on
    plot(t(1:200),z2(1:200),'r');

    return

end