Files
imdd_silas/Theory/Optical/PMD_maxwellshape/dgd_considerations.m
2026-03-25 10:57:48 +01:00

117 lines
2.6 KiB
Matlab

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