clc
clear all
tic
nn=input('nn=')
theta1=3*pi/8; %input('theta1=')
detphi=0*pi/50; %input('detphi')
% Given Parameters
gama=3; % 1/W/km
theta2=theta1+pi/2-0.05;
NUM=1024;
T=30; %ps
t=linspace(-T,T,NUM+1);
%qu=0.5*exp(-t(1:NUM).^2/100000/2); %quasi-CW,pulse width 1ns=1000ps
qu=0.1*rand(1,NUM);
%qu=0.1*randn(1,1024);
%qu=0.1*imnoise(1:1024,'gaussian');
%qu=wgn(1,1024,0.01); % Gaussian noise
Dt=2*T/NUM;
w=linspace(-pi./Dt,pi./Dt,NUM+1);
w1=fftshift(w(1:NUM));
tt=linspace(-T,T,NUM);
% the units are W,m,ps
h=0.1*1e-3; % calculation step unit km
for m=1:nn % nn number of round trips
u=qu*cos(theta1);
v=qu*sin(theta1)*exp(-i*detphi);
%%% propogation in standard fiber,5m long
z=5*1e-3;
Lb=2*1e-3;%12/5*1e-3; % according to the paper it's ~10m long
beta=pi/Lb;
deltaj=beta*1.554*1e-9/(2*pi*3e-7);
disp2=23; % ps^2/km
for m1=1:fix(z/h)
D1=exp((-deltaj*(i*w1)+i*0.5*disp2*(i*w1).^2)*h/2);
u1=fft(u).*D1;
u2=ifft(u1);
D2=exp((deltaj*(i*w1)+i*0.5*disp2*(i*w1).^2)*h/2);
v1=fft(v).*D2;
v2=ifft(v1);
N1=i*gama*((abs(u)).^2+2*(abs(v)).^2/3+conj(u).*v.^2./(u+eps+i*eps)/3)+i*beta;
N2=i*gama*((abs(v)).^2+2*(abs(u)).^2/3+conj(v).*u.^2./(v+eps+i*eps)/3)-i*beta;
u3=u2.*exp(N1.*h);
v3=v2.*exp(N2.*h);
D4=exp((-deltaj*(i*w1)+i*0.5*disp2*(i*w1).^2)*h/2);
u4=fft(u3).*D4;
u5=ifft(u4);
D5=exp((deltaj*(i*w1)+i*0.5*disp2*(i*w1).^2)*h/2);
v4=fft(v3).*D5;
v5=ifft(v4);
u=u5;
v=v5;
end
%% propogation in erbium-doped fiber,4m long
z=4*1e-3;
Lb=2*1e-3;%12/5*1e-3; % according to the table 1,L/Lb=5;
beta=pi/Lb;
deltaj=beta*1.554*1e-9/(2*pi*3e-7);
disp2=15; % ps^2/km
omegag=2*pi*4; % 4THz*2pi,change to the unit ps
g0=1600; %? unit? 80
Es=300;%1000;
for m2=1:fix(z/h)2
a=(abs(u)).^2+(abs(v)).^2;
gw=sum((abs(u)).^2)+sum((abs(v)).^2)*Dt;
%gw=san1(a,T,NUM);
gk=(1+gw/Es);
gp=g0./(1+(w(1:NUM)/(omegag)).^6)/gk; % according to equation (10)
D1=exp((-deltaj*(i*w1)+i*0.5*disp2*(i*w1).^2)*h/2);
u1=fft(u).*D1;
u2=ifft(u1);
D2=exp((+deltaj*(i*w1)+i*0.5*disp2*(i*w1).^2)*h/2);
v1=fft(v).*D2;
v2=ifft(v1);
N1=i*gama*((abs(u)).^2+2*(abs(v)).^2/3+conj(u).*v.^2./(u+eps+i*eps)/3)+i*beta+gp/2;
N2=i*gama*((abs(v)).^2+2*(abs(u)).^2/3+conj(v).*u.^2./(v+eps+i*eps)/3)-i*beta+gp/2;
u3=u2.*exp(N1.*h);
v3=v2.*exp(N2.*h);
D4=exp((-deltaj*(i*w1)+i*0.5*disp2*(i*w1).^2)*h/2);
u4=fft(u3).*D4;
u5=ifft(u4);
D5=exp((deltaj*(i*w1)+i*0.5*disp2*(i*w1).^2)*h/2);
v4=fft(v3).*D5;
v5=ifft(v4);
u=u5;
v=v5;
end
%%% propogation in standard single-mode fiber,3m long
z=3*1e-3;
Lb=2*1e-3;%12/5*1e-3Lb=12/5*1e-3; % 10m long
beta=pi/Lb;
deltaj=beta*1.554*1e-9/(2*pi*3e-7);
disp2=23;
for m3=1:fix(z/h)
D1=exp((-deltaj*(i*w1)+i*0.5*disp2*(i*w1).^2)*h/2);
u1=fft(u).*D1;
u2=ifft(u1);
D2=exp((deltaj*(i*w1)+i*0.5*disp2*(i*w1).^2)*h/2);
v1=fft(v).*D2;
v2=ifft(v1);
N1=i*gama*((abs(u)).^2+2*(abs(v)).^2/3+conj(u).*v.^2./(u+eps+i*eps)/3)+i*beta;
N2=i*gama*((abs(v)).^2+2*(abs(u)).^2/3+conj(v).*u.^2./(v+eps+i*eps)/3)-i*beta;
u3=u2.*exp(N1.*h);
v3=v2.*exp(N2.*h);
D4=exp((-deltaj*(i*w1)+i*0.5*disp2*(i*w1).^2)*h/2);
u4=fft(u3).*D4;
u5=ifft(u4);
D5=exp((deltaj*(i*w1)+i*0.5*disp2*(i*w1).^2)*h/2);
v4=fft(v3).*D5;
v5=ifft(v4);
u=u5;
v=v5;
end
qu=cos(theta2).*u+sin(theta2).*v;
P1=(abs(u)).^2+(abs(v)).^2+2*abs(u).*abs(v).*cos(angle(u)-angle(v));
P2=(abs(u)).^2+(abs(v)).^2-2*abs(u).*abs(v).*cos(angle(u)-angle(v));
qu2=P1;
qu3=P2;
end
figure(1)
plot(w(1:NUM),qu2);
figure(2)
plot(t(1:NUM),qu3);
figure(3)
plot(w(1:NUM),conj(qu).*qu);
toc
%plot(w(1:NUM),abs(fftshift(fft(qu))).^2/max(abs(fftshift(fft(qu))).^2),'b')