%%%%%%%%%%%%%%%%%%%%%%%%
% 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  問題3
%程序目的
%在kc處各放置一個高斯脈沖源hy
%在(kc-1)處放置另一個高斯脈沖源(-hy)
%%%%%%%%%%%%%%%%%%%%%%%%

clear %清除工作空間變量

%%%%%%%%%%%%%%%%%%%%%%%%
% 程序3 個初始參數
ke=200;  % ke是對模擬區域所取的節點（計算點）的數目
T=0;   %時間置零
nsteps=500;  % 循環的總次數  即模擬的總時間
%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%
%計劃在模擬區域的中心放置一個脈沖源
%此為脈沖源的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處各放置一個高斯脈沖源hy
    %在(kc-1)處放置另一個高斯脈沖源(-hy)
    pulse=exp(-0.5*power((t0-T)/spread,2));
    hy(kc)=pulse;
    hy(kc-1)=-hy(kc);
    
    %計算磁場部分
    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);
%      
%    }
% 
% }