%Simulace sireni dvou impulsu optickym vlaknem s PMD
%impulsy na blizkych vlnovych delkach -> chromaticka disperze, nelinearni
%koef vlakna a utlum vlakna stejny pro oba pulsy 
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=1;          %disperze, ps/nm*km
d1=-20e-16;     %koeficienty rozdilu grup rychlosti
d2=5e-16;
d3=30e-16;
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=12.5;        %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=200;       %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
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 1
B=sqrt(P0)*exp(-(1+i*C)*(t/T0).^(2*m));             %gauss puls 2
Ax=A/2;
Ay=A/2;
Bx=B/2;
By=B/2;
if  kk==1,
    subplot(1,2,1);
    plot(t*1e12,abs(A)/sqrt(P0),'b');
    xlabel('time [ps]');
    ylabel('normalized power');
    title('lambda 1');
    grid on;hold on;
    subplot(1,2,2);
    plot(t*1e12,abs(B)/sqrt(P0),'b');
    xlabel('time [ps]');
    ylabel('normalized power');
    title('lambda 2');
    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);
spBx=fft(Bx);
spBy=fft(By);
for k=1:h,
    spAx=spAx.*exp(-(alfa/2)*(dz/2)-i*w*(d1/2)*(dz/2)+...
         i*w.^2*(beta2/2)*(dz/2));
    spAy=spAy.*exp(-(alfa/2)*(dz/2)+i*w*(d1/2)*(dz/2)+...
         i*w.^2*(beta2/2)*(dz/2));
    spBx=spBx.*exp(-(alfa/2)*(dz/2)-i*w*(d2+d1/2)*(dz/2)+...
         i*w.^2*(beta2/2)*(dz/2));
    spBy=spBy.*exp(-(alfa/2)*(dz/2)-i*w*(d3+d1/2)*(dz/2)+...
         i*w.^2*(beta2/2)*(dz/2));
    Ax=ifft(spAx);
    Ay=ifft(spAy);
    Bx=ifft(spBx);
    By=ifft(spBy);
    Ax=Ax.*exp(i*gamma*((abs(Ax).^2)+2*(abs(Bx).^2)+...
       (2/3)*(abs(Ay).^2)+(2/3)*(abs(By).^2))*dz);
    Ay=Ay.*exp(i*gamma*((abs(Ay).^2)+2*(abs(By).^2)+...
       (2/3)*(abs(Ax).^2)+(2/3)*(abs(Bx).^2))*dz);
    Bx=Bx.*exp(i*gamma*((abs(Bx).^2)+2*(abs(Ax).^2)+...
       (2/3)*(abs(By).^2)+(2/3)*(abs(Ay).^2))*dz);
    By=By.*exp(i*gamma*((abs(By).^2)+2*(abs(Ay).^2)+...
       (2/3)*(abs(Bx).^2)+(2/3)*(abs(Ax).^2))*dz);
    spAx=fft(Ax);
    spAy=fft(Ay);
    spBx=fft(Bx);
    spBy=fft(By);
    spAx=spAx.*exp(-(alfa/2)*(dz/2)-i*w*(d1/2)*(dz/2)+...
         i*w.^2*(beta2/2)*(dz/2));
    spAy=spAy.*exp(-(alfa/2)*(dz/2)+i*w*(d1/2)*(dz/2)+...
         i*w.^2*(beta2/2)*(dz/2));
    spBx=spBx.*exp(-(alfa/2)*(dz/2)-i*w*(d2+d1/2)*(dz/2)+...
         i*w.^2*(beta2/2)*(dz/2));
    spBy=spBy.*exp(-(alfa/2)*(dz/2)-i*w*(d3+d1/2)*(dz/2)+...
         i*w.^2*(beta2/2)*(dz/2));
    Ax=ifft(spAx);
    Ay=ifft(spAy);
    Bx=ifft(spBx);
    By=ifft(spBy);    
   if k>1,          %navrat do vychozi pozice vektoru x,y
       Ax1=Ax;
       Ay1=Ay;
       Bx1=Bx;
       By1=By;       
       Ax=[cos(theta(1)) -exp(i*fi(1))*sin(theta(1))]*[Ax1;Ay1];
       Ay=[exp(-i*fi(1))*sin(theta(1)) cos(theta(1))]*[Ax1;Ay1];
       Bx=[cos(theta(2)) -exp(i*fi(2))*sin(theta(2))]*[Bx1;By1];
       By=[exp(-i*fi(2))*sin(theta(2)) cos(theta(2))]*[Bx1;By1];
   end;
   if k<h,          %nahodna elipticka rotace vektoru x,y
       theta=2*pi*rand(1,2);
       fi=2*pi*rand(1,2);
       Ax1=Ax;
       Ay1=Ay;
       Bx1=Bx;
       By1=By;
       Ax=[cos(theta(1)) exp(i*fi(1))*sin(theta(1))]*[Ax1;Ay1];
       Ay=[-exp(-i*fi(1))*sin(theta(1)) cos(theta(1))]*[Ax1;Ay1];           
       Bx=[cos(theta(2)) exp(i*fi(2))*sin(theta(2))]*[Bx1;By1];
       By=[-exp(-i*fi(2))*sin(theta(2)) cos(theta(2))]*[Bx1;By1];
   end;
    spAx=fft(Ax);
    spAy=fft(Ay);
    spBx=fft(Bx);
    spBy=fft(By);
end;
Ax=ifft(spAx);
Ay=ifft(spAy);
Bx=ifft(spBx);
By=ifft(spBy);
A=abs(Ax)+abs(Ay);
B=abs(Bx)+abs(By);
subplot(1,2,1);
plot(t*1e12,A/sqrt(P0),'r');
subplot(1,2,2);
plot(t*1e12,B/sqrt(P0),'r');
end;
pp=num2str(pp);
pp=[pp 'x output'];
for ll=1:2,
    subplot(1,2,ll);
    legend('input',pp);
end;