tlzscannomatrix.m

由于非线性光学中计算非线性折射系数的数值模拟

%constant
clc;
clear;
h=6.626176e-34;                              %普朗克常数 
c=2.99792458e8;                              %光速
Na=6.022045e23;
Hd=50e-6;                                    %样品厚度2mm
l=Hd/30;                                      %切片
%z-scan parameter
Tl=0.4;
alfa_1=-log(Tl)/Hd;
alfa_2=0e-10;                                 %非线性吸收系数
Re=1.42;                                      %线性折射率
beita_2=0e-17;                               %非线性折射率                            
lamda=532e-9;                                %激光波长
k=2*pi/lamda;                                %波矢量
omiga=6e-6;                                 %焦点处束腰半径
tau=8e-9;                                    %脉冲宽度    
f=pi*omiga^2/lamda;
w=0.55e-6;		                             %单脉冲能量
E0=sqrt(4*sqrt(log(2))*w/pi/sqrt(pi)/omiga^2/tau); %电场振幅
d=0.4;
dR=1e-3/10;
R=[0:dR:1e-3]';
% thermal parameter
cp=1.71;     %specific heat        J/(g.K)
% K=0.188;     %thermal conductivity  W/(m.K)
p=0.867e6;    %denstity        g/m^3
Dt=0.97e-3;  %thermal expansion coefficient
% D=K/cp/p;   %thermal diffusion    m^2/s
v=1170/sqrt(1.35);      %sound velocity (m/s)
es=(Re*Re-1)*(Re*Re+2)/3;
bb=es*Dt/2/Re;
% program parameter
dt=tau/30;                                   %时间积分步长
s=1;
 zstep=[-1000 -800 -600 -400 -300 -250 -200 -180 -160 -140 -110 -90 -70 -50 -30 -15 0 15 30 50 70 90 110 140 160 180 200 250 300 400 600 800]'/10;
for  z=0;%[-1000 -800 -600 -400 -300 -250 -200 -180 -160 -140 -110 -90 -70 -50 -30 -15 0 15 30 50 70 90 110 140 160 180 200 250 300 400 600 800]*1e-4,  
  wz=omiga*sqrt(1+z^2/f/f);                  %光腰
  dr=wz/30;
  ac1(:,:)=-log(Tl)/Hd*ones(4*round(wz/dr)+3,3*tau/dt+3);
  dn1(:,:)=zeros(4*round(wz/dr)+3,3*tau/dt+3);
  a=[0:1:4*round(wz/dr)]';
  b=[0:1:3*round(tau/dt)];                           %从-4~+4脉宽
  r=(a+0.5)*dr;
  r1=(a+1)*dr;
  r2=a*dr;
  Ez_in=zeros(4*round(wz/dr)+3,3*round(tau/dt)+3);
  Iz_in=zeros(4*round(wz/dr)+3,3*round(tau/dt)+3);
  Ez_in(3:4*round(wz/dr)+3,3:3*round(tau/dt)+3)=E0*GSBM(lamda,omiga,z,r)*(exp(-(b*dt-1.5*tau).*(b*dt-1.5*tau)*2/tau/tau));%E(r,z)
  for n=1:1:1;%Hd/l,  
      Ez_out=zeros(4*round(wz/dr)+3,3*round(tau/dt)+3);
      Iz_out=zeros(4*round(wz/dr)+3,3*round(tau/dt)+3);
      dn0=zeros(4*round(wz/dr)+3,3*round(tau/dt)+3);
      T=zeros(4*round(wz/dr)+3,3*round(tau/dt)+3);
      Q=zeros(4*round(wz/dr)+3,3*round(tau/dt)+3);
  for a1=3:1:4*round(wz/dr)+3,
    Iz_in(a1,:)=(abs(Ez_in(a1,:))).^2;                   %initial laser intensity
  ac2(a1,:)=alfa_1+alfa_2*Iz_in(a1,:);
  dc(a1,:)=beita_2*Iz_in(a1,:);
    Q(a1,:)=ac2(a1,:).*Iz_in(a1,:);
     for b1=3:1:3*round(tau/dt)+3,
        T(a1,b1)=sum(Q(a1,1:b1),2)*dt/cp/p;
     end %b1
     end %a1
         for a2=3:1:4*round(wz/dr)+3,   
          for b2=3:1:3*round(tau/dt)+3, 
     x1(a2,b2)=[dn0(a2-1,b2)-2*dn0(a2-1,b2-1)+dn0(a2-1,b2-2)]/dt/dt/v/v; 
     x2(a2,b2)=(-2*r(a2-2)*dn0(a2-1,b2)+r2(a2-2)*dn0(a2-2,b2)+...
         2*(r1(a2-2)*dn0(a2,b2-1)-2*r(a2-2)*dn0(a2-1,b2-1)+r2(a2-2)*dn0(a2-2,b2-1))+...
         r1(a2-2)*dn0(a2,b2-2)-2*r(a2-2)*dn0(a2-1,b2-2)+r2(a2-2)*dn0(a2-2,b2-2))/4/r(a2-2)/dr/dr;
     x3(a2,b2)=(r1(a2-2)*T(a2,b2)-2*r(a2-2)*T(a2-1,b2)+r2(a2-2)*T(a2-2,b2)+...
         2*(r1(a2-2)*T(a2,b2-1)-2*r(a2-2)*T(a2-1,b2-1)+r2(a2-2)*T(a2-2,b2-1))+...
         r1(a2-2)*T(a2,b2-2)-2*r(a2-2)*T(a2-1,b2-2)+r2(a2-2)*T(a2-2,b2-2))*bb/4/r(a2-2)/dr/dr;
     dn0(a2,b2)=(x1(a2,b2)-x2(a2,b2)-x3(a2,b2))*4*r(a2-2)*dr*dr/r1(a2-2);
        end %a2
    end %b2
    dn2=dc+dn0;break
%     for b3=3:1:3*round(tau/dt)+3,
%         for a3=3:1:4*round(wz/dr)+3,
%             Y1=2*j*k*Re*(Ez_out(a3-1,b3)-Ez_in(a3-1,b3))/l-...
%                 k*k*(((dn2(a3-1,b3)+Re)^2-Re*Re)*Ez_out(a3-1,b3)+((dn1(a3-1,b3)+Re)^2-Re*Re)*Ez_in(a3-1,b3))/2+...
%                 j*k*Re*(ac2(a3-1,b3)*Ez_out(a3-1,b3)+ac1(a3-1,b3)*Ez_in(a3-1,b3))/2;
%             Y2=(-2*r(a3-2)*Ez_out(a3-1,b3)+r2(a3-2)*Ez_out(a3-2,b3)+...
%                 r1(a3-2)*Ez_in(a3,b3)-2*r(a3-2)*Ez_in(a3-1,b3)+r2(a3-2)*Ez_in(a3-2,b3))/2/r(a3-2)/dr/dr;
%             Ez_out(a3,b3)=(Y1-Y2)*2*r(a3-2)*dr*dr/r1(a3-2);
%          end %a3
%         end  %b3
%         Ez_in=Ez_out;
%         ac1=ac2;
%         dn1=dn2;
     end    %n
    E=Ez_out(3:4*round(wz/dr)+3,3:3*round(tau/dt)+3);
    Iz_out=(abs(E)).^2;
    Pz_out(s)=r'*sum(Iz_out,2)*dt*dr*2*pi;
    hk=exp(-i*pi*R.^2/lamda/(d-z-Hd))*(exp(-i*pi*r'.^2/lamda/(d-z-Hd)).*r');
    hf=besselj(0,2*pi*R*r'/lamda/(d-z-Hd)).*hk*exp(-i*k*(d-z-Hd)/Re);   
    Ez_s=i*hf*E*2*pi*dr/lamda/(d-z-Hd);               %小孔距离焦点0.8米
    Iz_r=(abs(Ez_s)).^2;
    Pz_r(s)=R'*sum(Iz_r,2)*dt*dR*2*pi;
    s=s+1; 
   end %[-Z-Z]
   %从折射中去掉总吸收
  T4=Pz_r./Pz_out;
  %对总的吸收进行归一化
   All=Pz_out'/Pz_out(1);

   %对除去总吸收的折射进行归一化
   S=T4'/T4(1);
save S S -ascii;
save zstep zstep -ascii;
save All All -ascii;
plot(zstep,S,zstep,All)