%Simulace sireni impulsu optickym vlaknem s PMD
close all; clear all; clc;
i=sqrt(-1);c=3e8;       %konstanty
%vstupni parametry--------------------------------------------------------
P0=0.001;       %vstupni spickovy vykon gauss pulsu, W
alfadB=0;       %utlum vlakna, dB/km
Dcd=17;         %disperze, ps/nm*km
Dpmd=3;         %koef PMD, ps/sqrt(km)
lambda=1550;    %vlnova delka, nm
n2=3.2e-16;     %nelinearni koef indexu lomu vlakna, cm^2/W
Aef=50;         %efektivni plocha jadra vlakna, um^2
z=50;           %delka vlakna, km
dz=0.25;        %element vzdalenosti v rozmezi 0.1 - 1 km
T0=50;          %sirka pulsu v ps, kde A=A0*1/e
C=0;            %koef rozmitaneho gauss pulsu
m=1;            %rad gauss pulsu
N=1024;         %pocet bodu casoveho vektoru
okno=400;       %sirka casoveho okna, ps
pp=5;           %pocet vystupnich pulsu
%-------------------------------------------------------------------------
beta2=-Dcd*1e-24*((lambda)^2)/(2*pi*c);%disperzni koef 2.r, s^2/m
gamma=(2*pi*n2*1e17)/(lambda*Aef);  %nelinearni koef vlakna, 1/W*m
d=sqrt(3*pi/8)*Dpmd*1e-15/sqrt(dz); %koef rozdilu grup rychlosti
T0=(T0/2)*1e-12;
h=z/dz;                                             %pocet kroku vypoctu
dz=dz*1000;                                         %km -> m
alfa=(2*alfadB/4.343)/1000;                         %utlum vlakna 1/m
t=-(okno/2)*1e-12:(okno*1e-12)/(N-1):(okno/2)*1e-12;%casovy vektor
dt=(okno*1e-12)/(N-1);                              %casovy element
for kk=1:pp,
A=sqrt(P0)*exp(-(1+i*C)*(t/T0).^(2*m));             %gauss puls
Ax=A/2;
Ay=A/2;
if  kk==1,
    figure(1);
    plot(t*1e12,abs(A)/sqrt(P0),'b');
    xlabel('time [ps]');
    ylabel('normalized power');
    grid on;hold on;
    l=length(A);
    dw=2*pi*1/(dt*l);
    w=(-1*l/2:1:l/2-1)*dw;
    w=fftshift(w);
end;
spAx=fft(Ax);
spAy=fft(Ay);
for k=1:h,
    spAx=spAx.*exp(-(alfa/2)*(dz/2)+i*w*(d/2)*(dz/2)+...
         i*w.^2*(beta2/2)*(dz/2));
    spAy=spAy.*exp(-(alfa/2)*(dz/2)-i*w*(d/2)*(dz/2)+...
         i*w.^2*(beta2/2)*(dz/2));
    Ax=ifft(spAx);
    Ay=ifft(spAy);
    Ax=Ax.*exp(i*gamma*((abs(Ax).^2)+(2/3)*(abs(Ay).^2))*dz);
    Ay=Ay.*exp(i*gamma*((abs(Ay).^2)+(2/3)*(abs(Ax).^2))*dz);
    spAx=fft(Ax);
    spAy=fft(Ay);
    spAx=spAx.*exp(-(alfa/2)*(dz/2)+i*w*(d/2)*(dz/2)+...
         i*w.^2*(beta2/2)*(dz/2));
    spAy=spAy.*exp(-(alfa/2)*(dz/2)-i*w*(d/2)*(dz/2)+...
         i*w.^2*(beta2/2)*(dz/2));
    Ax=ifft(spAx);
    Ay=ifft(spAy);
   if k>1,          %navrat do vychozi pozice vektoru x,y
       Ax1=Ax;
       Ay1=Ay;
       Ax=[cos(theta) -exp(i*fi)*sin(theta)]*[Ax1;Ay1];
       Ay=[exp(-i*fi)*sin(theta) cos(theta)]*[Ax1;Ay1];
   end;
   if k<h,          %nahodna elipticka rotace vektoru x,y
           theta=2*pi*rand(1);
           fi=2*pi*rand(1);
           Ax1=Ax;
           Ay1=Ay;
           Ax=[cos(theta) exp(i*fi)*sin(theta)]*[Ax1;Ay1];
           Ay=[-exp(-i*fi)*sin(theta) cos(theta)]*[Ax1;Ay1];           
    end;
    spAx=fft(Ax);
    spAy=fft(Ay);    
end;
Ax=ifft(spAx);
Ay=ifft(spAy);
A=abs(Ax)+abs(Ay);
plot(t*1e12,A/sqrt(P0),'r');
end;
pp=num2str(pp);
pp=[pp 'x output'];
legend('input',pp);