lineardifferent.m

计算超短脉冲非线性传输的一个比较成熟的M文件

%有限差分法求解非线性传输的laplace算符项，单位：cm，s
%****************************************************
            %参数值
%****************************************************
format long;
clear all;
tic;
c=3e+10;
lamda0=0.00008;%the central wavelength
nb=1.0;%the refraction index of the medium.
n2=4.1e-16;%the nonlinear refraction index.
k0=2*pi/lamda0;%the wave number
k=nb*k0;
omega0=2*pi*c/lamda0;%frequency of the pulse.
w0=75e-4;% the beam waist'radius;
T_FWHM=70e-15;
tp=T_FWHM/sqrt(2*log(2));
R=500e-4;%the transverse length;
L0=400;%the  propagation length.
%****************************************************
               %步长定义
%****************************************************
n_t=2^10;
T=1.2e-12;
d_t=T/n_t;
kt=1:n_t;
t=(kt-n_t/2)*d_t;
n_r=2^8;
n_z=3000;
d_r=R/(n_r-1);
d_z=L0/(n_z-1);
kr=0:n_r-1;
r=kr.*d_r;
kz=0:d_z:L0;
d_f=1/T;
J=0:1:n_t-1;
F=(J-n_t/2)*d_f;
omega=2*pi*F;
%*****************************************************
                 %光束初始值
%*****************************************************
Pcr=lamda0^2/2/pi/n2;%光束的临界功率
Pin=1.17*Pcr;%输入功率为临界功率的1.17倍
E0=sqrt(2*Pin/pi/w0^2);
EE=zeros(n_r,n_t,n_z/100)*1e-30;
rho=zeros(n_r,n_t);
C=zeros(n_r,n_r);
%ez0=zeros(n_r,n_t)*1e-30;
%***************************************************
                  %非线性项参数
%***************************************************
k2=248e-30;%GVD系数
ee=1.60210e-19;%电子电荷C
me=9.1091e-31;%电子质量kg
yibixitong0=8.854e-12;%真空介电常数F/m
rhocritical=yibixitong0*me*omega0^2/ee/ee;%临界电子密度
deltak=1.2e-52;%多光子电离系数
rhoatom=6.68e22;%中性原子数目
Ui=6.5*ee;%电离势
h=1.05457266e-34;%约化普朗克常数
kk=(Ui-rem(Ui,h*omega0))/h/omega0;
betak=k*h*omega0*rhoatom*deltak;% K光子吸收非线性系数
%****************************************************
                    %线性laplace算符部分的计算
%****************************************************
%for n=1:n_r
 %   r=(n-1)*d_r;
  %  for m=1:n_t;
   %     t=-T/2+(m-1)*d_t;
        E=E0*exp(-r'.^2./w0^2)*exp(-t.^2./tp^2);
    %end
%end
A=sqrt(-1)/2/k*d_z/d_r^2/2;%前半个步长用后向差分
B=sqrt(-1)/2/k*d_z/4/d_r;
%C=-k*n2*(abs(E).^2)+k*r(n,i)./2./rhocritical-betak*i.*(abs(E(n,i)).^(2*kk-2))./2;
DN1=sparse(1:n_r,1:n_r,1-2*A,n_r,n_r);
DN2=sparse(2:n_r,1:n_r-1,A-B./(((1:n_r-1)*d_r)+1e-30),n_r,n_r);
DN3=sparse(1:n_r-1,2:n_r,A+B./(((0:n_r-2)*d_r)+1e-30),n_r,n_r);
HN10=DN1+DN2+DN3;
HN10(1,1)=1-4*A;
HN10(1,2)=4*A;
%HN10(n_r,n_r-1)=A-B;
%HN10(n_r-1,n_r)=0;%?

DM1=sparse(1:n_r,1:n_r,1+2*A,n_r,n_r);%后半个步长用前向差分 
DM2=sparse(2:n_r,1:n_r-1,-A+B./(((1:n_r-1)*d_r)+1e-30),n_r,n_r);
DM3=sparse(1:n_r-1,2:n_r,-A-B./(((0:n_r-2)*d_r)+1e-30),n_r,n_r);
HM10=DM1+DM2+DM3;
HM10(1,1)=1+4*A;
HM10(1,2)=-4*A;
%HM10(n_r,n_r-1)=-A+B;
%HM10(n_r-1,n_r)=0;
%IHM10=inv(HM10);
%E(n_r,:)=0;

%***************************************************
              %循环计算开始
%***************************************************
 for q=1:20;
     for m=2:n_t;
    K1=(rhoatom-rho(:,m-1)).*deltak.*(abs(E(:,m)).^2).^kk;
    K2=(rhoatom-rho(:,m-1)-K1.*d_t).*deltak.*(abs(E(:,m)).^2).^kk;
    rho(:,m)=rho(:,m-1)+d_t.*(K1+K2)./2;%rung-kuta法求电子密度矩阵
     end
     size(r)
     for mm=1:n_r;
%       C(n_r,n_r)=-k.*n2.*(abs(E).^2)+k.*r(mm,:)./2./rhocritical-betak.*i.*(abs(E).^(2*kk-2))./2;
       C(n_r,n_r)=-k.*n2.*(abs(E).^2)+k.*r(mm,:)./2./rhocritical;
     end
     DM1=DM1+C;
     DN1=DN1+C;
     HN10=DN1+DN2+DN3;
     HN10(1,1)=HN10(1,1)+C;
     HM10=DM1+DM2+DM3;
     HM10(1,1)=HM10+C;
     IHM10=inv(HM10);
     for j=1:n_r;
         E(j,:)=fftshift(fft(E(j,:),n_t));
         E_DL(j,:)=E(j,:).*exp(-k2*i*(i.*omega).^2.*d_z./2);
         E_TDL(j,:)=ifft(ifftshift(E_DL(j,:)));
     end
    
       for i=1:n_t;
         E_r(:,i)=IHM10*(HN10*E(:,i));
       end
     E=E_r+E_TDL-E;
   if rem(q,100)==0;
       q;
       time=toc;
   EE(:,:,q/100)=E;
   end
 end 
 %x=(-T/2:d_t:T/2);
 %y=(0:d_r:R);
 %[X,Y]=meshgrid(x,y);
 %I=real(E).^2+imag(E).^2;
surf(E)