gaussiansource.m

时域有限差分中平面光源及高斯光源的设置.对光源设置有一定意义

  clear
 clc

 %************************************************************
 %fundamental constants
 %*************************************************************
cc=2.99792458e8 ;            
muz=4.0*pi*1.0e-7;          
epsz=1.0/(cc*cc*muz);       

    % parameters of wave    
lambda=780e-9; 
freq=cc./lambda;                 
omega=2.0*pi*freq
%***************************************************************
% grid parameters
%***************************************************************
   
   %  total grid
 ie=120 ; 
 je=120  ;
  ib=ie+1;
 jb=je+1;
  % space increment and time step and total number of time steps
 dx=20e-9 ;   
 dt=dx/(2.0*cc) ;
 nmax=100000;    
 
 
 
 
%************************************************************
%  material parameters and coefficients
%************************************************************
  
  
   %material parameters
  eps=1.0 ;
  sig=0.0; 
  mur=1.0;
  sim=0.0;
  
  %coefficients
   
 
   
       eaf=dt*sig/(2.0*epsz*eps);
       ca=(1.0-eaf)/(1.0+eaf);
       cb=dt/epsz/eps/dx/(1.0+eaf);
       haf=dt*sim/(2.0*muz*mur);
       da=(1.0-haf)/(1.0+haf);
       db=dt/muz/mur/dx/(1.0+haf);
  
 
 %***********************************************************
 % specification the coefficients
 %***********************************************************
  
  %initialize the all coefficients 
  caex(1:ib,1:je)=ca;
  cbex(1:ib,1:je)=cb;
  
  caey(1:ie,1:jb)=ca;
  cbey(1:ie,1:jb)=cb;
  
  dahz(1:ib,1:jb)=da;
  dbhz(1:ib,1:jb)=db;
  
  %add the special material's coefficients in!
      %position of the special material
     
  
%*************************************************************  
%  initialization of field arrays
%************************************************************ 
 ex=zeros(ib,je);  %fields in main grid
 ey=zeros(ie,jb);
 hz=zeros(ib,jb);

 
 %*************************************************************
 % wave  source
 %*************************************************************
  slength=20:100 ;
  rwaist=55;
  scenter=slength(1)+(length(slength)-1)/2  ;         
  
  
  source=zeros(ib,jb);
  source(10,slength)=exp(-(slength-scenter).^2/rwaist^2);
 
 %*************************************************************
 % initialization of visualizing
 %***********************************************************
 
  pcolor(hz');
  shading flat;
  %caxis ([-1000 1000]);
  caxis auto;
  axis([1 ib 1 jb]);
  colorbar;
  axis image;
 % axis off;
  opengl neverselect;
  title(['ex at time step=0']);
  
  
   
  
  %****************************************************************
  % begain time stepping!!
  %****************************************************************
  
  
  
   for n=1:nmax
        
       %add the source
        
              if n<pi/omega/dt
                  switchf=0.5.*(1-cos(n*dt*omega));
                  else
                  switchf=1;
              end              
            
           % hz=hz+exp(-j*(omega.*n.*dt))*source*switchf;
       
      %update electric fields in main   grid
        ex(:,:)=caex(:,:).*ex(:,:)+...
            cbex(:,:).*(hz(:,2:jb)-hz(:,1:je));
        
        ey(:,:)=caey(:,:).*ey(:,:)+...
            cbey(:,:).*(hz(1:ie,:)-hz(2:ib,:));
        
        hzabc(:,:)=hz(:,:);  %for the boundary
        
        hz(2:ie,2:je)=dahz(2:ie,2:je).*hz(2:ie,2:je)-...
            dbhz(2:ie,2:je).*((ey(2:ie,2:je)-ey(1:ie-1,2:je))...
            -(ex(2:ie,2:je)-ex(2:ie,1:je-1)));
        
        
                  
        %update electric fields on boundary 
              
          hz(1,2:je)=hzabc(2,2:je)+((cc.*dt-dx)/(cc.*dt+dx))*(hz(2,2:je)-hz(1,2:je))...
              +cc.^2.*epsz*eps(1)*dt/(2*(cc.*dt+dx))*(ex(1,2:je)-ex(1,1:je-1)+ex(2,2:je)-ex(2,1:je-1));
        
         hz(ib,2:je)=hzabc(ie,2:je)+((cc.*dt-dx)/(cc.*dt+dx))*(hz(ie,2:je)-hz(ib,2:je))...
              +cc.^2.*epsz*eps(1)*dt/(2*(cc.*dt+dx))*(ex(ib,2:je)-ex(ib,1:je-1)+ex(ie,2:je)-ex(ie,1:je-1));
          hz(2:ie,1)=hzabc(2:ie,2)+((cc.*dt-dx)/(cc.*dt+dx))*(hz(2:ie,2)-hz(2:ie,1))...
              -cc.^2.*epsz*eps(1)*dt/(2*(cc.*dt+dx))*(ey(2:ie,1)-ey(1:ie-1,1)+ey(2:ie,2)-ey(1:ie-1,2));
         hz(2:ie,jb)=hzabc(2:ie,je)+((cc.*dt-dx)/(cc.*dt+dx))*(hz(2:ie,je)-hz(2:ie,jb))...
              -cc.^2.*epsz*eps(1)*dt/(2*(cc.*dt+dx))*(ey(2:ie,jb)-ey(1:ie-1,jb)+ey(2:ie,je)-ey(1:ie-1,je));
         
           %update electric fields in corner
          
            hz(1,1)=hzabc(2,2)+((cc.*dt-sqrt(2)*dx)/(cc.*dt+sqrt(2)*dx))*(hz(2,2)-hz(1,1)); 
           
            hz(1,jb)=hzabc(2,je)+((cc.*dt-sqrt(2)*dx)/(cc.*dt+sqrt(2)*dx))*(hz(2,je)-hz(1,jb));
            
            hz(ib,jb)=hzabc(ie,je)+((cc.*dt-sqrt(2)*dx)/(cc.*dt+sqrt(2)*dx))*(hz(ie,je)-hz(ib,jb));
           
            hz(ib,1)=hzabc(ie,2)+((cc.*dt-sqrt(2)*dx)/(cc.*dt+sqrt(2)*dx))*(hz(ie,2)-hz(ib,1));
        
            
            
            
            hz=hz+exp(-j*(omega.*n.*dt))*source*switchf;
          %visualize fields                         
             
          if mod(n,10)==0;
             
            fprintf('第 %d 时间步已完成\n',n);               
            timestep=int2str(n);
%              hzabs=abs(hz');
              hzabs=abs(hz');
              hzabs=abs(ex');
             pcolor(hzabs);
         
             shading flat;
             caxis auto;
             axis([1 ie 1 jb]);
             colorbar;
             axis image;
            % axis off;
             title(['ex at time step=',timestep]);
             
             disp('正在计算......');
             
           end
             
            pause(0.000001)
        
         
   end
    
   %*******************************************************
   % other work ,for example ,save datas
   %*********************************************************
   
   %*********************************************************
   % end
   %**********************************************************