fdtd3d.c

matlab编写的三维FDTD程序 pml边界条件 实现电磁波传播数值模拟和仿真

  
    ======================================================================== */

for(t=0;t<fd1.timnum;t++) {

   zend=zcent+10+1.58*t/pow(fd1.zperm,0.35);
   if(zend>zm) zend=zm;


/* =========================================================================

   SPECIFICATION OF ELECTROMAGNETIC ENERGY INPUT INTO GRID DURING EACH TIMESTEP
   
   Eamp = GAUSSIAN VOLTAGE PULSE 

   ========================================================================= */

   /* USING A GAUSSIAN PULSE */

   Eamp=Vo*exp(-(t-tcval)*(t-tcval)*(fd1.dt*fd1.dt)/(timvar*timvar))/(fd1.dx);
   
                  
/* =========================================================================

   STORE OLD VALUES OF E AND H NEAR GRID EDGES AND SIDES FOR USE IN ABC
   EQUATIONS 

   ========================================================================= */

    for(y=0;y<=ym;y++) {
       for(z=0;z<=zend;z++) {
          Ezxolda[0][y][z]=Ez[0][y][z];
          Eyxolda[0][y][z]=Ey[0][y][z];
          Ezxmaxa[0][y][z]=Ez[xm][y][z];
          Eyxmaxa[0][y][z]=Ey[xm][y][z];
          Ezxolda[1][y][z]=Ez[1][y][z];
          Eyxolda[1][y][z]=Ey[1][y][z];
          Ezxmaxa[1][y][z]=Ez[xm-1][y][z];
          Eyxmaxa[1][y][z]=Ey[xm-1][y][z];
          Hzxolda[0][y][z]=Hz[1][y][z];
          Hyxolda[0][y][z]=Hy[1][y][z];
          Hzxmaxa[0][y][z]=Hz[xm][y][z];
          Hyxmaxa[0][y][z]=Hy[xm][y][z];
          Hzxolda[1][y][z]=Hz[2][y][z];
          Hyxolda[1][y][z]=Hy[2][y][z];
          Hzxmaxa[1][y][z]=Hz[xm-1][y][z];
          Hyxmaxa[1][y][z]=Hy[xm-1][y][z];
       }
    }
    for(x=0;x<=xm;x++) {
      for(z=0;z<=zend;z++) {
         Ezyolda[x][0][z]=Ez[x][0][z];
         Exyolda[x][0][z]=Ex[x][0][z];
         Ezymaxa[x][0][z]=Ez[x][ym][z];
         Exymaxa[x][0][z]=Ex[x][ym][z];
         Ezyolda[x][1][z]=Ez[x][1][z];
         Exyolda[x][1][z]=Ex[x][1][z];
         Ezymaxa[x][1][z]=Ez[x][ym-1][z];
         Exymaxa[x][1][z]=Ex[x][ym-1][z];
         Hzyolda[x][0][z]=Hz[x][1][z];
         Hxyolda[x][0][z]=Hx[x][1][z];
         Hzymaxa[x][0][z]=Hz[x][ym][z];
         Hxymaxa[x][0][z]=Hx[x][ym][z];
         Hzyolda[x][1][z]=Hz[x][2][z];
         Hxyolda[x][1][z]=Hx[x][2][z];
         Hzymaxa[x][1][z]=Hz[x][ym-1][z];
         Hxymaxa[x][1][z]=Hx[x][ym-1][z];
      }
   }
   for(x=0;x<=xm;x++) {
      for(y=0;y<=ym;y++) {
         Eyzolda[x][y][0]=Ey[x][y][0];
         Eyzolda[x][y][1]=Ey[x][y][1];
         Exzolda[x][y][0]=Ex[x][y][0];
         Exzolda[x][y][1]=Ex[x][y][1];
         Eyzmaxa[x][y][0]=Ey[x][y][zm];
         Eyzmaxa[x][y][1]=Ey[x][y][zm-1];
         Exzmaxa[x][y][0]=Ex[x][y][zm];
         Exzmaxa[x][y][1]=Ex[x][y][zm-1];
         Hyzolda[x][y][0]=Hy[x][y][1];
         Hyzolda[x][y][1]=Hy[x][y][2];
         Hxzolda[x][y][0]=Hx[x][y][1];
         Hxzolda[x][y][1]=Hx[x][y][2];
         Hyzmaxa[x][y][0]=Hy[x][y][zm];
         Hyzmaxa[x][y][1]=Hy[x][y][zm-1];
         Hxzmaxa[x][y][0]=Hx[x][y][zm];
         Hxzmaxa[x][y][1]=Hx[x][y][zm-1];
      }
    }

/* ========================================================================

    CALCULATE CURRENT AND VOLTAGE AT CENTER OF ANTENNAS FOR OUTPUT TO FILE(S) 

   ========================================================================  */

   Iyt=(Hx[xc][yc][zcent-5]-Hx[xc][yc+1][zcent-5])*fd1.dy+
        (Hy[xc+1][yc][zcent-5]-Hy[xc][yc][zcent-5])*fd1.dx;
   v=Ey[xc][yc+1][zcent-5]*fd1.dx; 
   Zx=0.0;
   if(rec1flag==1) {
      Iyr=(Hx[xcere][yc][zcere-5]-Hx[xcere][yc+1][zcere-5])*fd1.dy+
        (Hy[xcere+1][yc][zcere-5]-Hy[xcere][yc][zcere-5])*fd1.dx;
   }

/* =========================================================================

   RESTORE H-FIELD VALUES AT JUNCTION BETWEEN COAX EXTENDING ABOVE GRID AND 
      PORTION OF COAX WITHIN GRID 

   ========================================================================= */

   if(tra1flag>0) {
      Hx[xc1][yc1][1]  =Hxc1[0][zcm+1];
      Hx[xc1][yc1+1][1]=Hxc1[1][zcm+1];
      Hx[xc2][yc2][1]  =Hxc2[0][zcm+1];
      Hx[xc2][yc2+1][1]=Hxc2[1][zcm+1];

      Hy[xc1][yc1][1]    =Hyc1[0][zcm+1];
      Hy[xc1+1][yc1][1]  =Hyc1[1][zcm+1];
      Hy[xc2][yc2][1]    =Hyc2[0][zcm+1];
      Hy[xc2+1][yc2][1]  =Hyc2[1][zcm+1];
   }
   if(tra2flag>0) {
      Hx[xc3][yc3][1]  =Hxc3[0][zcm+1];
      Hx[xc3][yc3+1][1]=Hxc3[1][zcm+1];
      Hx[xc4][yc4][1]  =Hxc4[0][zcm+1];
      Hx[xc4][yc4+1][1]=Hxc4[1][zcm+1];

      Hy[xc3][yc3][1]    =Hyc3[0][zcm+1];
      Hy[xc3+1][yc3][1]  =Hyc3[1][zcm+1];
      Hy[xc4][yc4][1]    =Hyc4[0][zcm+1];
      Hy[xc4+1][yc4][1]  =Hyc4[1][zcm+1];
   }
   if(rec1flag>0) {
      Hx[xcr1][ycr1][1]  =Hxcr1[0][zcm+1];
      Hx[xcr1][ycr1+1][1]=Hxcr1[1][zcm+1];
      Hx[xcr2][ycr2][1]  =Hxcr2[0][zcm+1];
      Hx[xcr2][ycr2+1][1]=Hxcr2[1][zcm+1];
  
      Hy[xcr1][ycr1][1]    =Hycr1[0][zcm+1];
      Hy[xcr1+1][ycr1][1]  =Hycr1[1][zcm+1];
      Hy[xcr2][ycr2][1]    =Hycr2[0][zcm+1];
      Hy[xcr2+1][ycr2][1]  =Hycr2[1][zcm+1];
   }
 /*
  if(rec2flag>0) {
     Hx[xcr3][ycr3][1]  =Hxcr3[0][zcm+1];
     Hx[xcr3][ycr3+1][1]=Hxcr3[1][zcm+1];
     Hx[xcr4][ycr4][1]  =Hxcr4[0][zcm+1];
     Hx[xcr4][ycr4+1][1]=Hxcr4[1][zcm+1];

     Hy[xcr3][ycr3][1]    =Hycr3[0][zcm+1];
     Hy[xcr3+1][ycr3][1]  =Hycr3[1][zcm+1];
     Hy[xcr4][ycr4][1]    =Hycr4[0][zcm+1];
     Hy[xcr4+1][ycr4][1]  =Hycr4[1][zcm+1];
  }
 */

/* =========================================================================

    ADVANCE E-FIELD VALUES IN COAX ABOVE GRID  
 
   ========================================================================= */

    for(z=0;z<=zcm;z++) {
      if(tra1flag>0) {
         Exc1[0][z] = eoldmulc*Exci1uolda[0][z]+enmulc*(Hzc1[0][1][z]-Hzc1[0][0][z]+Hyc1[0][z]-Hyc1[0][z+1]);
         Exc1[1][z] = eoldmulc*Exci1uolda[1][z]+enmulc*(Hzc1[1][1][z]-Hzc1[1][0][z]+Hyc1[1][z]-Hyc1[1][z+1]);
         Exci1uolda[0][z]=Exc1[0][z];
         Exci1uolda[1][z]=Exc1[1][z];

         Eyc1[0][z] = eoldmulc*Eyci1uolda[0][z]+enmulc*(Hxc1[0][z+1]-Hxc1[0][z]+Hzc1[0][0][z]-Hzc1[1][0][z]);
         Eyc1[1][z] = eoldmulc*Eyci1uolda[1][z]+enmulc*(Hxc1[1][z+1]-Hxc1[1][z]+Hzc1[0][1][z]-Hzc1[1][1][z]);
         Eyci1uolda[0][z]=Eyc1[0][z];
         Eyci1uolda[1][z]=Eyc1[1][z];

         Exc2[0][z] = eoldmulc*Exci2uolda[0][z]+enmulc*(Hzc2[0][1][z]-Hzc2[0][0][z]+Hyc2[0][z]-Hyc2[0][z+1]);
         Exc2[1][z] = eoldmulc*Exci2uolda[1][z]+enmulc*(Hzc2[1][1][z]-Hzc2[1][0][z]+Hyc2[1][z]-Hyc2[1][z+1]);
         Exci2uolda[0][z]=Exc2[0][z];
         Exci2uolda[1][z]=Exc2[1][z];

         Eyc2[0][z] = eoldmulc*Eyci2uolda[0][z]+enmulc*(Hxc2[0][z+1]-Hxc2[0][z]+Hzc2[0][0][z]-Hzc2[1][0][z]);
         Eyc2[1][z] = eoldmulc*Eyci2uolda[1][z]+enmulc*(Hxc2[1][z+1]-Hxc2[1][z]+Hzc2[0][1][z]-Hzc2[1][1][z]);
         Eyci2uolda[0][z]=Eyc2[0][z];
         Eyci2uolda[1][z]=Eyc2[1][z];
      }
      if(tra2flag>0) {
         Exc3[0][z] = eoldmulc*Exci3uolda[0][z]+enmulc*(Hzc3[0][1][z]-Hzc3[0][0][z]+Hyc3[0][z]-Hyc3[0][z+1]);
         Exc3[1][z] = eoldmulc*Exci3uolda[1][z]+enmulc*(Hzc3[1][1][z]-Hzc3[1][0][z]+Hyc3[1][z]-Hyc3[1][z+1]);
         Exci3uolda[0][z]=Exc3[0][z];
         Exci3uolda[1][z]=Exc3[1][z];

         Eyc3[0][z] = eoldmulc*Eyci3uolda[0][z]+enmulc*(Hxc3[0][z+1]-Hxc3[0][z]+Hzc3[0][0][z]-Hzc3[1][0][z]);
         Eyc3[1][z] = eoldmulc*Eyci3uolda[1][z]+enmulc*(Hxc3[1][z+1]-Hxc3[1][z]+Hzc3[0][1][z]-Hzc3[1][1][z]);
         Eyci3uolda[0][z]=Eyc3[0][z];
         Eyci3uolda[1][z]=Eyc3[1][z];

         Exc4[0][z] = eoldmulc*Exci4uolda[0][z]+enmulc*(Hzc4[0][1][z]-Hzc4[0][0][z]+Hyc4[0][z]-Hyc4[0][z+1]);
         Exc4[1][z] = eoldmulc*Exci4uolda[1][z]+enmulc*(Hzc4[1][1][z]-Hzc4[1][0][z]+Hyc4[1][z]-Hyc4[1][z+1]);
         Exci4uolda[0][z]=Exc4[0][z];
         Exci4uolda[1][z]=Exc4[1][z];

         Eyc4[0][z] = eoldmulc*Eyci4uolda[0][z]+enmulc*(Hxc4[0][z+1]-Hxc4[0][z]+Hzc4[0][0][z]-Hzc4[1][0][z]);
         Eyc4[1][z] = eoldmulc*Eyci4uolda[1][z]+enmulc*(Hxc4[1][z+1]-Hxc4[1][z]+Hzc4[0][1][z]-Hzc4[1][1][z]);
         Eyci4uolda[0][z]=Eyc4[0][z];
         Eyci4uolda[1][z]=Eyc4[1][z];
      }
      if(rec1flag>0) {
         Excr1[0][z] = eoldmulc*Excri1uolda[0][z]+enmulc*(Hzcr1[0][1][z]-Hzcr1[0][0][z]+Hycr1[0][z]-Hycr1[0][z+1]);
         Excr1[1][z] = eoldmulc*Excri1uolda[1][z]+enmulc*(Hzcr1[1][1][z]-Hzcr1[1][0][z]+Hycr1[1][z]-Hycr1[1][z+1]);
         Excri1uolda[0][z]=Excr1[0][z];
         Excri1uolda[1][z]=Excr1[1][z];

         Eycr1[0][z] = eoldmulc*Eycri1uolda[0][z]+enmulc*(Hxcr1[0][z+1]-Hxcr1[0][z]+Hzcr1[0][0][z]-Hzcr1[1][0][z]);
         Eycr1[1][z] = eoldmulc*Eycri1uolda[1][z]+enmulc*(Hxcr1[1][z+1]-Hxcr1[1][z]+Hzcr1[0][1][z]-Hzcr1[1][1][z]);
         Eycri1uolda[0][z]=Eycr1[0][z];
         Eycri1uolda[1][z]=Eycr1[1][z];

         Excr2[0][z] = eoldmulc*Excri2uolda[0][z]+enmulc*(Hzcr2[0][1][z]-Hzcr2[0][0][z]+Hycr2[0][z]-Hycr2[0][z+1]);
         Excr2[1][z] = eoldmulc*Excri2uolda[1][z]+enmulc*(Hzcr2[1][1][z]-Hzcr2[1][0][z]+Hycr2[1][z]-Hycr2[1][z+1]);
         Excri2uolda[0][z]=Excr2[0][z];
         Excri2uolda[1][z]=Excr2[1][z];

         Eycr2[0][z] = eoldmulc*Eycri2uolda[0][z]+enmulc*(Hxcr2[0][z+1]-Hxcr2[0][z]+Hzcr2[0][0][z]-Hzcr2[1][0][z]);
         Eycr2[1][z] = eoldmulc*Eycri2uolda[1][z]+enmulc*(Hxcr2[1][z+1]-Hxcr2[1][z]+Hzcr2[0][1][z]-Hzcr2[1][1][z]);
         Eycri2uolda[0][z]=Eycr2[0][z];
         Eycri2uolda[1][z]=Eycr2[1][z];
       }
      /*
      if(rec2flag>0) {
         Excr3[0][z] = eoldmulc*Excri3uolda[0][z]+enmulc*(Hzcr3[0][1][z]-Hzcr3[0][0][z]+Hycr3[0][z]-Hycr3[0][z+1]);
         Excr3[1][z] = eoldmulc*Excri3uolda[1][z]+enmulc*(Hzcr3[1][1][z]-Hzcr3[1][0][z]+Hycr3[1][z]-Hycr3[1][z+1]);
         Excri3uolda[0][z]=Excr3[0][z];
         Excri3uolda[1][z]=Excr3[1][z];

         Eycr3[0][z] = eoldmulc*Eycri3uolda[0][z]+enmulc*(Hxcr3[0][z+1]-Hxcr3[0][z]+Hzcr3[0][0][z]-Hzcr3[1][0][z]);
         Eycr3[1][z] = eoldmulc*Eycri3uolda[1][z]+enmulc*(Hxcr3[1][z+1]-Hxcr3[1][z]+Hzcr3[0][1][z]-Hzcr3[1][1][z]);
         Eycri3uolda[0][z]=Eycr3[0][z];
         Eycri3uolda[1][z]=Eycr3[1][z];

         Excr4[0][z] = eoldmulc*Excri4uolda[0][z]+enmulc*(Hzcr4[0][1][z]-Hzcr4[0][0][z]+Hycr4[0][z]-Hycr4[0][z+1]);
         Excr4[1][z] = eoldmulc*Excri4uolda[1][z]+enmulc*(Hzcr4[1][1][z]-Hzcr4[1][0][z]+Hycr4[1][z]-Hycr4[1][z+1]);
         Excri4uolda[0][z]=Excr4[0][z];
         Excri4uolda[1][z]=Excr4[1][z];

         Eycr4[0][z] = eoldmulc*Eycri4uolda[0][z]+enmulc*(Hxcr4[0][z+1]-Hxcr4[0][z]+Hzcr4[0][0][z]-Hzcr4[1][0][z]);
         Eycr4[1][z] = eoldmulc*Eycri4uolda[1][z]+enmulc*(Hxcr4[1][z+1]-Hxcr4[1][z]+Hzcr4[0][1][z]-Hzcr4[1][1][z]);
         Eycri4uolda[0][z]=Eycr4[0][z];
         Eycri4uolda[1][z]=Eycr4[1][z];
      }
     */
  }
  

/* =========================================================================

   FORCE E-FIELD TOP OF TRANSMIT COAX CABLE TO CORRESPOND WITH VOLTAGE SPIKE 

   ========================================================================= */

   if(Eamp>0.01){ 
      if(tra1flag==1) {
         Eyc1[1][1]=Eamp;
         Eyc1[0][1]=(-1.)*Eamp;
         Exc1[1][1]=Eamp;
         Exc1[0][1]=(-1.)*Eamp;

         Eyc2[1][1]=(-1.)*Eamp;
         Eyc2[0][1]=(1.)*Eamp;
         Exc2[1][1]=(-1.)*Eamp;
         Exc2[0][1]=(1.)*Eamp;

         Exci1uolda[0][1]=Exc1[0][1];
         Exci1uolda[1][1]=Exc1[1][1];
         Eyci1uolda[0][1]=Eyc1[0][1];
         Eyci1uolda[1][1]=Eyc1[1][1];

         Exci2uolda[0][1]=Exc2[0][1];
         Exci2uolda[1][1]=Exc2[1][1];
         Eyci2uolda[0][1]=Eyc2[0][1];
         Eyci2uolda[1][1]=Eyc2[1][1];
      }
      if(tra2flag==1) {
         Eyc3[1][1]=(-1.)*Eamp;
         Eyc3[0][1]=(-1.)*(-1.)*Eamp;
         Exc3[1][1]=(-1.)*Eamp;
         Exc3[0][1]=(-1.)*(-1.)*Eamp;

         Eyc4[1