clear;
tr=300;tl=-300; %ps
Tmax=tr-tl;% window size
N=2^14;% FFT points
dt=(tr-tl)/N; %序列采样间隔
c=3*10^5;%light velocity nm/ps
% % % % % % % % % % % % % % % %
t=dt*(-N/2:1:N/2-1);% temporal grid
omega=2*pi/Tmax*(-N/2:1:N/2-1);%2pi*k/N*dt
df=1/Tmax;
f=(-N/2:N/2-1)*df;
wavelength0=1060;% center wavelength nm
frequence0=c/wavelength0;
frequence=f+frequence0;
lambda=c./frequence;% nm
% % % % % % % % % % % % % % % % AM
mf =10; %Ghz modulation window is 50ps, which means frequency is 20GHz.
T=1000/mf;%单个脉冲占据的时间窗口 ps
b = 0;%相对偏置 2*Vbias/Vpi
m = 1;% 调制深度 2*Vm/Vpi Vm为微波调制电压的幅度
omega_intense = 2*pi*mf/1000;% modulation frequensy is 20GHz.
nn = min(2*round(500/mf*(N/Tmax)),N);
piecewise = [zeros(1,(N-nn)/2),ones(1,nn),zeros(1,(N-nn)/2)];
T_intense = 0.5*(1+sin(pi/2*(b...
+m*cos(omega_intense*t)))).*piecewise;
% PM 1
b_phase=0;
m_phase=10;% 相位调制深度 V/Vpi，每个调制器最大调制5个pi
%% 不同调制深度
MR1=[];
EpR1=[];
swidthR1=[];
TwidthR1=[];
for loop=1:1
omega_phase = 2*pi*mf/1000;%
T_phase = exp(1i*(pi*(b_phase...
+m_phase*cos(omega_phase*t)))).*piecewise;
%%fiber
%% 不同色散
L1=3;
MR=[];
EpR=[];
swidthR=[];
TwidthR=[];
for loop1=1:1
L2=1.5;
L3=4;
%% YDF
g0=2; %small signal gain parameter (1/m)
Psat=50/1000; %yb fiber saturable power (w)
beta2=0.023;%ps^2/m D=-38.5
%%%%%%%%%%%%%%%
dispersion1=exp((1i/2*beta2*omega.^2)*L1);
dispersion2=exp((1i/2*beta2*omega.^2)*L3);
%input
% uu=ones(1,length(t));
% uu=uu+i*uu;
aa=unifrnd(0,1,[1,N]);
bb=unifrnd(0,1,[1,N]);
uu=0.1*aa+0.1*1i*bb;
Ep=trapz(t,abs(uu).^2);
Puu=Ep/T;
answer=[];
%% 腔内循环
for loop2=1:3000;
u1=uu;
%% AM
uu=uu.*T_intense;
%% PM
uu=uu*0.5;
uu=uu.*T_phase;
%% fiber1
Utemp=dispersion1.*fftshift(ifft(fftshift(uu)));
uu=fftshift(fft(fftshift(Utemp)));
%% iso
uu=uu*0.63;%2db 腔内损耗和增益的平衡 只要总损耗一样 放在这里还是最后结果都是差不多的。虽然有一点点不同，是因为YDF的增益和他入射初始光强有关，所以放大不同
%% YDF
for loop4=1:10
Ep=trapz(t,abs(uu).^2);%pulse enegy (pj)
Puu=Ep/T;%脉冲平均功率 (W)
g=g0/(1+Puu/Psat);
Utemp2=exp(((1i/2*(beta2+1i*g*2.2230e-004)*omega.^2)+g/2)*L2/10).*fftshift(ifft(fftshift(uu)));%GVD fft并把前后部分对调
uu=fftshift(fft(fftshift(Utemp2)));
end
%% OC
uu=uu.*sqrt(0.9);
uc=uu;
uc2=uu.*sqrt(0.1);
%% fiber2
Utemp2=dispersion2.*fftshift(ifft(fftshift(uu)));
uu=fftshift(fft(fftshift(Utemp2)));
%% lost
% uu=uu*0.63;% 3db PL and IM
if mod(loop,10)==0
answer=[answer,uc'];
end
% if max(abs(uu).^2-abs(u1).^2)/max(abs(uu).^2)<=10^(-8)
% break
% end
if sqrt(sum((abs(uu).^2-abs(u1).^2).^2))...
/sum(abs(uu).^2) <= 10^(-9)
break
end
size=fix(loop/10);
end
% L1=L1+1;
X=ones(N,size);
Y=ones(N,size);
for n=1:(size)
X(:,n) = t;
end
for n=1:(size)
Y(:,n) = n';
end
%% change the output point
uu=uc;
% uu=uc2;
figure(1)
plot3(X,Y,1000*abs(answer).^2,'b')
xlabel(' t ps');
ylabel(' loop/10 ')
zlabel('power mw')
figure(2);
plot(t,1000*abs(uu).^2,'-r');
xlabel('t ps');
ylabel('power mw');
axis([-100 100 0 inf]);
figure(3)
plot(lambda,10*log10(abs(fftshift(ifft(uu))).^2*10.^3));
xlabel('\lambda nm');
ylabel('Spectral Power dbm ');
axis([1040 1080 -80 -20]);
figure(4)
plot(lambda,1000*abs(fftshift(ifft(uu))).^2);
xlabel('\lambda nm');
ylabel('Spectral Power mw ');
axis([1050 1070 0 inf]);
figure(5);
%plot(t(1:N-1),unwrap(-diff(phase(uu))/dt));%相邻两点差值大于2*pi则翻转至，差值在2*pi内
ph=phase(uu);
plot(t(1:N-1),-diff(phase(uu))/dt);
axis([-50 50 -inf inf]);
xlabel('t ps');
ylabel('chirp ');
figure(6);
plot(t,T_intense);
xlabel('t ps');
ylabel('modulator T');
%% 显示参数
%峰值功率 mw
M=max(1000*abs(uu).^2)
MR=[MR,M];
%单脉冲能量 pj
Ep=trapz(t,abs(uu).^2)
EpR=[EpR,Ep];
%3dB脉宽 ps
Tm=max(abs(uu).^2)/2;
point=find(abs(uu).^2>=Tm);
Twidth=t(max(point))-t(min(point))
TwidthR=[TwidthR,Twidth];
%20dB 带宽 nm
% Sm=max(10*log10(abs(fftshift(ifft(u))).^2*10.^3))-20;
% spoint=find(10*log10(abs(fftshift(ifft(u))).^2*10.^3)>=Sm);
% Swidth=lambda(min(spoint))-lambda(max(spoint))
%3dB频宽 nm
sm2=max(abs(fftshift(ifft(uu))).^2)/2;
point=find(abs(fftshift(ifft(uu))).^2>=sm2);
swidth=lambda(min(point))-lambda(max(point))
swidthR=[swidthR,swidth];
end
%% 记录二维数据
% MR=MR';
% EpR=EpR';
% TwidthR=TwidthR';
% swidthR=swidthR';
%% 记录三维数据
% MR1=[MR1,MR];
% EpR1=[EpR1,EpR];
% TwidthR1=[TwidthR1,TwidthR];
% swidthR1=[swidthR1,swidthR];
% m_phase=m_phase+0.5;
end
%% 画等高图
% xx=[0.5:0.5:10];%%20元素 x轴 调制深度
% yy=[6.5:35.5]*-38.5/1000;% 30元素 y轴 色散D
% [C,h]=contourf(bb,b,swidthR1,60);
% xlabel('modulate depth pi');
% ylabel('dispersion ps/nm ');
