cIc tic 1=3*pi/8 %input( theta=') detphi=0*pi/50: %input( detphi') Gi P theta2=thetal+pi/2-0 NUM=1024 t=linspace(-T, T, NUM+1) %qu=0. 5kexp(t(1: NUM). 2/100000/2): %quasi-CW, pulse width Ins=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) wl=fftshift(w(1: NUM)) tt=linspace(T, T, NUM) the units are W, m, ps h=0. 1*le-3: calculation step unit ki for m=l: nn nn number of round trips etal) v=qu*sin(theta1)*exp (-i*detphi) %% propogation in standard fiber, 5m long Lb=2*le-3: %12/5*1e-3: according to the paper it's 10m long beta=pi/Lb
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) ps for m1=1: fix(z/h) DI=exp ((-deltaj(i*w1)+i*0. 5*disp2*(i=*w1).2)*h/2) 1=fft(u).*DI 2=ifft(ul) D2=exp((deltaj* (i*w1)+i=0. 5*disp2*(i*w1).2)*h/2) vl=fft(v).*D2 N1=i*gama*((abs(u)).2+2*(abs(v)).2/3+conj(u).*V.2./(u+eps+i*eps)/3)+i*eta; N2=i*gama*((abs(v)).2+2*(abs(u)).2/3+conj(v).*u. 2./(v+eps+i*eps)/3)-isbeta u3=12.*exp(N1.*h) 3=V2.*exp(N2.*h) D4=exp((-delta j*(i=sw1)+i=0. 5*disp2*i*w1).2)*h/2) u4=fft(u3).*D4 D5=exp((deltaj* (i*w1)+i=0. 5*disp2*(i*w1).2)*h/2) 4=fft(v3).*D5; (v4) % propogation in erbium-doped fiber, 4m long b=2*1e-3: %12/5*1e-3: according to the table 1, L/Lb=5 beta=pi/Lb deltaj=beta*l 554*le-9/(2*pi*3e-7) omega o the unit g0=1600 Es=300:‰1000
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;
or m2=1: fix(z/h)2 a=(abs(u)).2+(abs(v)).2 w=sum((abs(u)).2)+sum((abs(v)).2)*Dt %gw=sanl(a, T, NUM) gk=(1+gw/Es) gp=gO /(1+(w(1: NUM)/(omega)).6)/gk: according to equation(10) Dl=exp(- deltas*(i*W1)+i*0.5*dsp2*(*1).2)*h/2); ul=fft(u).**DI 2=ifft(ul) D2=exp(+ deltas(i*W1)+1*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./(uteps+i*eps)/3)+i*beta+gp 12=i*gama*((abs(v)). 2+2*(abs(u)).2/3+conj(v).*u. 2. /(v+eps+isps)/3)-i*betatgp u3=u12.*eXp(N1.*h) v3=v2. *exp (N2 *h) D4=exp((- deltas*(i*W1)+1*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) v5=ifft(v4) u=u5 %% propogation in standard single-mode fiber, 3m long Lb=2*1e-3;%12/5*1e-3Lb=12/5*1e-3;%10m1ong beta=pi/ Lb delta j=beta*1. 554*le-9/(2*pi*3e-7)
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) DI=exp((-deltaj*(i=w1)+i*0. 5*disp2* (i*w1).2)*h/2) ul=fft(u).*D1 u2=ifft(ul) D2=exp((delta j*(i*w1)+i=0. 5*disp2*(i*w1).2)*h/2) vl=fft(v).*D2 v2=ifft(v1) N1=i*gam*(abs(u).2+2*(abs(v).2/3+conj(u).*V.2./(ueps+i*eps)/3)+i*beta; N2=i*gama*((abs(v)).2+2*(abs(u).2/3+conj().*u. 2./(v+eps+i*eps)/3)-i*beta u3=12.*exp(N1.*h) v3=v2. *exp (N2 *h) D4=exp(- deltas*(i*W1)+i*0.5*disp2*(*W1).2)*h/2); u4=fft(u3).*4 u5=ifft(u4 5=exp((delta*(i=*w1)+i*0. 5*disp2*(i*w1).2)*h/2) 4=fft(v3).*D5; v5=ifft(v4) u=u5 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=Pl 3=P2 plot(w(1: NUM), qu2) plot(t(1: NUM), qu3) plot(w(1: NUM), conj (qu). *qu %plot(w(1: NUM), abs(fftshift(fft(qu))). 2/max(abs(fftshift(fft(qu))).2), b')
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')