fd1d_1_1_two_ex_source.m

程序目的 在kc处各放置一个高斯脉冲源hy 在(kc-1)处放置另一个高斯脉冲源(-hy)

%%%%%%%%%%%%%%%%%%%%%%%%
% ELECTROMAGNETIC SIMULATION USING THE FDTD METHOD
%     Dennis M.Sullivan 
% Electrical Engineering Department 
% University of  Idaho
% ISBN  0-7803-4747-1
% 书的配套程序的MATLAB版  本程序最后含有原C程序代码
%作者：徐飞
%2006年1月25日
%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%
%page 4 :PROBLEM SET 1.1  问题2
%程序目的
%在(kc-20)和(kc+20)处各放置一个高斯脉冲源ex
%%%%%%%%%%%%%%%%%%%%%%%%

clear %清除工作空间变量

%%%%%%%%%%%%%%%%%%%%%%%%
% 程序3 个初始参数
ke=200;  % ke是对模拟区域所取的节点（计算点）的数目
T=0;   %时间置零
nsteps=200;  % 循环的总次数  即模拟的总时间
%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%
%计划在模拟区域的中心放置一个脉冲源
%此为脉冲源的3个参数
kc=ke/2;  %计算模拟区域的中点
t0=40.0;    %脉冲源的中点 (我觉得是：时间中点，即脉冲出现最大值的时刻)
                   %参见吉林大学王涛硕士论文 17页
spread=12;  %脉冲源的宽度
%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%
%对ex和hy进行初始化，全部置0
ex=zeros(1,ke); 
hy=zeros(1,ke);
%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%
% 程序主体！！！！

    %%%%%%%%%%%%%%%%%%%%%%%%%%
    % %计算部分
for n=1:nsteps   % 循环的总次数  即模拟的总时间
    T=T+1;   %记录循环的总次数（其实和n是一样的）
    
    %计算电场部分
    for k=2:ke
        ex(k)=ex(k)+0.5*(hy(k-1)-hy(k));  %书上公式1.9a
    end

    %在(kc-20)和(kc+20)处各放置一个高斯脉冲源ex
    pulse=exp(-0.5*power((t0-T)/spread,2));
    ex(kc-20)=pulse;
    ex(kc+20)=pulse;
    
    %计算磁场部分
    for k=1:ke-1
        hy(k)=hy(k)+0.5*(ex(k)-ex(k+1));    %书上公式1.9b
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%
    % %绘图部分
    % 实时画出电场Ex和磁场Hy随时间的变化的情况
    subplot(2,1,1)
        plot(ex)
        axis([0 ke -2 2])  
        
    subplot(2,1,2)
        plot(hy)
        axis([0 ke -2 2]) 
     %%%%%%%%%%%%%%%%%%%%%%%%
    pause(0.001)

end

%最后加上网格和说明
subplot(2,1,1)  
    grid on   
    title('电场Ex')

subplot(2,1,2)  
    grid on    
    title('磁场Hy')


%%%%%%%%%%%%%%%%%%%%%%%
%下面是原书配套的c程序   fd1d_1.1.c
%此程序从网上（非作者网站）down来的
% 可能和书上的程序略有不同，但是编译运行无错
%%%%%%%%%%%%%%%%%%%%%%%




% /* FD1D_1.c. 1D FDTD simulation in free space 在自由空间模拟 */
% # include <math.h>
% # include <stdlib.h>
% # include <stdio.h>
% /*  */
% #define KE 200 /* KE is number of cells to be used KE是模拟所用单元的个数 */
% main ()
% {
%         float ex[KE],hy[KE];    /* 定义电场ex 和 磁场 hy*/
%         int n,k,kc,NSTEPS=1;  /* 定义 */
%         float T=0;        /* 定义时间 */
%         float t0,spread,pulse;
%         FILE *fp, *fopen();
% 
%         /* initialize 初始化*/
%         
%   for (k=1; k<KE; k++)
%   {    ex[k] = 0.;
%        hy[k] = 0.;
%   }
%        kc=KE/2;          /* center of problem space */
%        t0=40.0;          /* center of pulse */
%        spread=12;        /* width of pulse  */
%        NSTEPS = 1;
%   while (NSTEPS > 0)
%   {
%      printf("NSTEPS -->") ;  /* NSTEPS is number of times */
%      scanf ("%d",&NSTEPS);   /* main loop has executed */
%      printf("%d",NSTEPS);
%      n=0;
% 
%      for (n=1; n<=NSTEPS ; n++)
%      {
%          T=T+1;          /* T keeps track of total number */
%                          /* Main FDTD Loop  */
%                          /* Calulate the Ex field */
%                          
%          for (k=1; k<KE ; k++)
%          ex[k]=ex[k] + .5*(hy[k-1]-hy[k]);
%          
%                          /* Put a Gaussian pulse in the middle */
%                          
%          pulse = exp(-.5*(pow((t0-T)/spread, 2.0) ));
%          ex[kc] = pulse;
%          printf("%5.1f  %6.2f\n",t0-T,ex[kc]);
%         
%                          /* calculate the hy field */
%                          
%          for (k=0; k< KE-1;k++)
%          hy[k]=hy[k]+.5*(ex[k] - ex[k+1]);
%      }
%                          /* end of the main FDTD loop */
%                          /* At the end of the calculation, print out the Ex and Hy fields */
%                          
%      for (k=1; k<=KE; k++)
%      printf ( "%3d  ex[%d]= %6.2f  hy[%d]=%6.2f\n",k,k,ex[k],k,hy[k]);
%      
%                          /* write the E field out to a file "Ex" */
%                          
%      fp=fopen ("Ex","w");
%      fprintf(fp,"ex[k]\n");
%      for (k=1; k<=KE; k++)
%      fprintf (fp," ex[%d]=  %6.2f \n", k,ex[k]);
%      fclose(fp);
%                          /* Write the H field out to a hile "Hy" */
%                          
%      fp=fopen("Hy","w");
%      fprintf(fp,"hy[k]\n");
%      for (k=1; k<=KE; k++)
%      fprintf(fp," hy[%d]= %6.2f \n",k,hy[k]);
%      fclose(fp);
%      printf( "T = %5.0f\n", T);
%      
%    }
% 
% }