oned_hyex_sourceh_dng.asv

一维双负介质的drude等效模型的防真程序．可以从输出图象中看出双负介质的基本特性

   %1D FDTD 双负介质
   % 定义计算参数
    nz=200;    %总计算格数
	TIME=1000;    %总计算时间
    hyi=zeros(1,nz);exi=zeros(1,nz+1);    %网格剖分
    kyi=zeros(1,nz);jxi=zeros(1,nz);     % 等效电磁流
     d=10;       %PML总层数
     m=2;        %PML中参数衰减速率阶数
     C=3.0e8;       %光速
     R0=1.0e-5;     %设定的波吸收率
     pi=3.1415926;
     ee=8.851e-12;
     mu=4*3.1415926e-7;
     f=30*10^9;
     dx=0.0005;dy=0.0005;dz=0.0005;dt=1*1e-12;      %网格剖分与时间剖分
     dez=zeros(1,nz+1);dhz=zeros(1,nz);       %网格中的介质参数
     
  
	c1=(m+1)*C*ee/2.0/(d*dz)*log(1/R0);		     %最大煤质电参数
	m1=mu*c1/ee;     %最大煤质磁参数
    
    %DNG频率设定
    wpe=5*1e11;wpm=5*1e11;
    w0=2*pi*f;
    gaie=1e8;gaim=1e8;
   
    
    %计算各网格点的介质参数
	 	   for i=1:nz+1 
			   if  i<=d
				   dez(i)=c1*(d+1-i)*(d+1-i)/d/d;
			   elseif i>=nz-d+2
				   dez(i)=c1*(i+d-nz-1)*(i+d-nz-1)/d/d;
			   else
				   dez(i)=0;
               end
           end
		   
		   for i=1:nz
			   if i<=d
				   dhz(i)=m1*(d-i+1.0/2.0)*(d-i+1.0/2.0)/d/d;
               elseif i>=nz-d+1
				   dhz(i)=m1*(i+d-nz-1.0/2.0)*(i+d-nz-1.0/2.0)/d/d;
               else
				   dhz(i)=0;
               end
           end
           
    %场值迭代
  %   for (t=0;t<=TIME;t++)    
    for t=0:TIME
   
%	   for(i=0;i<=nx-1;i++)
           for i=1:nz                                     %更新Hy
	   			if i==80
		        hyi(i)=(1-dhz(i)*dt/(2*mu))/(1+dhz(i)*dt/(2*mu))*hyi(i)+dt/(mu*(1+dhz(i)*dt/(2*mu)))*(exi(i)-exi(i+1))/dz+sin(2*pi*f*t*dt);
            elseif i>nz/2&i<nz*3/4
				hyi(i)=(1-dhz(i)*dt/(2*mu))/(1+dhz(i)*dt/(2*mu))*hyi(i)+dt/(mu*(1+dhz(i)*dt/(2*mu)))*(exi(i)-exi(i+1)-kyi(i)*dz)/dz;
            else                
				hyi(i)=(1-dhz(i)*dt/(2*mu))/(1+dhz(i)*dt/(2*mu))*hyi(i)+dt/(mu*(1+dhz(i)*dt/(2*mu)))*(exi(i)-exi(i+1))/dz;
               end
           end
           
           for i=nz/2+1:nz*3/4-1                                    %更新ky
	   			kyi(i)=(1-0.5*gaim*dt)/(1+0.5*gaim*dt)*kyi(i)+mu*wpm^2*dt/(1+0.5*gaim*dt)*hyi(i);
            end

           for i=2:nz                                 %更新ex
             if i>nz/2&i<nz*3/4
				exi(i)=(1-dez(i)*dt/2/ee)/(1+dez(i)*dt/2/ee)*exi(i)+dt/ee/(1+dez(i)*dt/2/ee)*(hyi(i-1)-hyi(i)-0.5*(jxi(i)+jxi(i-1))*dz)/dz; 
            else 
                exi(i)=(1-dez(i)*dt/2/ee)/(1+dez(i)*dt/2/ee)*exi(i)+dt/ee/(1+dez(i)*dt/2/ee)*(hyi(i-1)-hyi(i))/dz;  
            end
           end
           
          for i=nz/2+1:nz*3/4-1                                    %更新ky
	   			jxi(i)=(1-0.5*gaie*dt)/(1+0.5*gaie*dt)*jxi(i)+ee*wpe^2*dt/(1+0.5*gaie*dt)*0.5*(exi(i)+exi(i+1));
            end
            
           % % -------------    PLOT    -------------------
    figure(1);clf;
    set(gcf,'Renderer','OpenGL');     %使用windows系统自带的画笔速度更快

    h=plot(hyi);
    grid on;
      set(h,'LineWidth',3);
      set(gca,'FontName','TimesNewRoman');
      set(gca,'FontSize',14);
      h = xlabel('x向坐标');
      ylabel('磁场强度');
      hold on
      set(gca,'YLim',[-1.5 1.5]);
      hold on; %plot([N+1 N+1],1.1*get(gca,'YLim'),'-.'); 
      % METHOD NAMES
      x=[nz/2,nz/2,nz,nz];y=[1.5,-1.5,-1.5,1.5];
      fill(x,y,[0.8,0.8,0.5]);
      h=text(0.75*nz,.75*min(get(gca,'YLim')),'DNG');
      set(h,'FontName','Times');
      set(h,'FontSize',30);
      title(['时间步为:',num2str(t),'步']);
      
     % h=text(1.25*N,.75*min(get(gca,'YLim')),'DWG');
     % set(h,'FontName','Times');
     % set(h,'FontSize',20);
      hold off;
      drawnow;
    
    end