clear all
close all
% zhou.gt@163.com
%光纖參數          超連續譜的話 N=2^15  detat=10^-15 步長取小一點

wavelength=1550*(10^-9);     %wavelength 單位m
c=3*(10^8);                  %光速 單位:m/s
w0=2*pi*c/wavelength;        %入射波長中心
%%%%%%%%%%% 光線參數%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha=9/4.343*(10^-3);            %attenuation coefficient 單位:/m 損耗
                       %這個參數可以改 它決定入射脈沖的形狀(陡峭度)
% D=-0.7*(10^-6);              %色散參量   單位s/m^2%-0.7
beta2=0.72*10^-27;             %二階色散參數


beta3=0.02*(10^-39);         %三階色散參數  單位s^3/m

gama=11*(10^-3);                    %非線性項 單位：/w/m
%%%%%%%%%%% 光線參數%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

TR=3*(10^-15);                %拉曼增益，單位：s   3fs 
                        

%%%%%%%%%%% 入射脈沖參數%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                      
TFWHM=6.89*(10^-13);   %脈沖半高寬度，單位為s  
T0=TFWHM/1.763;
%T0=2*(10^-13);              %1/e 點寬度0.20皮秒
% P_ave=26.12*10^-3;             %輸入光平均功率，單位：w
% P=2^0.5*P_ave/TFWHM/50/10^6;
P=50;    %此處可以改參數  輸入峰值脈沖功率 單位：w
%%%%%%%%%%% 入射脈沖參數%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
P=P^0.5;
N=2^12;              %取樣點數 有可能需要改 但必須是2的偶數次方
n=0:(N-1);  
n1=0:N/2-1;
n2=-N/2:1:-1;
h=1e-3;              %步長，單位:m 這個地方可以改的
Nh=2e2;            %模擬步長數 這個地方可以改的
NN=10;               %記錄波形數
              

Detat=3*10^-15;          %抽樣點時間間隔 這個地方可以改的
T=N*Detat;   %總時間長度，單位為s
t=(n-N/2)*Detat;    %時間序列
t1=n1*Detat;
t2=n2*Detat;
Detaw=2*pi/(T);     %頻譜間隔
w=(n-N/2)*Detaw;          %頻譜序列
wave=2*pi*c./(w+w0)*10^6;     %波長序列um
D_w=-alpha/2+j*beta2*(w.^2)/2-j*beta3*(w.^3)/6;
UU1=D_w(N/2+1:N);               %改變頻率顯示順序
UU2=D_w(1:N/2);
D_w=[UU1,UU2];

 U_in=P*sech(t/T0); %輸入時間脈沖的包絡,Hyperbolic-Secant Pulse
    

U_mh=ones(NN+1,N);
U_mhw=ones(NN+1,N);

U=U_in;
Uw=ifft(U,N)*N;%傅立葉變換
ABS_Uw=Uw.*conj(Uw);
U_mh(1,:)=abs(U).^2;
U_mhw(1,:)=fftshift(ABS_Uw);
% T_rms=sum(t1.*U_mh(1,(1:N/2)))./sum(U_mh(1,(1:N/2)));%用來計算均方差脈寬的 
% T2_rms=sum(t1.^2.*U_mh(1,(1:N/2)))./sum(U_mh(1,(1:N/2)));
% sigma_T(1)=2*sqrt(T2_rms-T_rms.^2);
for k=1:Nh 
    ABS_U=abs(U).^2;
    ABS_Uw=Uw.*conj(Uw);

   if(mod(k,Nh/NN)==0)
        U_mh(k*NN/Nh+1,:)=abs(U).^2;
        U_mhw(k*NN/Nh+1,:)=fftshift(ABS_Uw);
%         T_rms=sum(t1.*U_mh(k*NN/Nh+1,(1:N/2)))./sum(Umh(11,(1:N/2)));%用來計算均方差脈寬的
%         T2_rms=sum(t1.^2.*U_mh(k*NN/Nh+1,(1:N/2)))./sum(Umh(11,(1:N/2)));
%         sigma_T(k*NN/Nh+1)=2*sqrt(T2_rms-T_rms.^2);
    end
        Selfcraggedness=i*[ABS_U(1)*U(1)-ABS_U(N)*U(N),diff(ABS_U.*U)]./(w0*Detat*U);%這個地方不知道是點除還是直接除 
        ABS_UU=[ABS_U(1)-ABS_U(N),diff(ABS_U)]./Detat;
        Raman_effect=TR*ABS_UU;
        NL=i*gama*(ABS_U+Selfcraggedness-Raman_effect);
        U=ifft(exp(h*NL).*U)*N;   %傅立葉變換
        U=fft(exp(h*D_w).*U)/N;%反傅立葉變換
        Uw=ifft(U)*N;
    end
    Maxw=max(U_mhw(1,:));
    Max=max(U_mh(1,:));
    Max11=max(U_mhw(11,:));
 %圖形輸出         


 figure   %輸出波長譜
  plot(wave,10*log10(U_mhw(11,:)/Max11));
 title('出射脈沖波長譜');
 xlabel('\lambda({\mu}m)');
 ylabel('相對強度');
  figure
 
  waterfall(w*T0,[0:NN]*Nh*h/NN,U_mhw/Maxw)
  xlabel('(\omega-{\omega}_{0})T_0')
  ylabel('Z(m)')
  zlabel('相對強度')
  title('頻譜波形')

%   figure
% plot(t/T0, U_mh(11,:));
% xlabel('T/T0');
% ylabel('|U(z,t)|^2');
% title('出射脈沖');
% figure
%plot(sigma_T/sigma_T(1));
%ylabel('sigma_T');