e_abc_iter.c

FDTD 模拟 天线 包含 近场－远场 外推程序（时谐场）

/*-------------------------------------------------
      
       The Iteration Formulae of E Component

       at the absorbing boundary condition        

*/


# include "E_ABc_Iter.h"

void E_Down_Iter(MODULE *pModule, double ****Ex, double ****Ey)
{

 int i,j;
 int NExx,NExy,NExz;
 int NEyx,NEyy,NEyz;

 double tmp,er,vel;

 NExx=Nx;    NExy=Ny+1;  NExz=Nz+1;
 NEyx=Nx+1;  NEyy=Ny;    NEyz=Nz+1;

     er=pModule->pPermit[0];
     vel=C/sqrt(er);  

     tmp=(vel*dt-dhz)/(vel*dt+dhz);

   /*  Iteration Formulae of Ex Points  */

    switch(pModule->nSurf[1])
     {
       case 0:  /* It means that Down Plane is electric wall */

         for(i=0;i<NExx;i++)
          {
            for(j=1;j<NExy-1;j++)
              {
                Ex[1][i][j][0]=0.0;
              }
           }
           
            break;

          default:  /* G.Mur. First Order */

            for(i=0; i<NExx; i++)
              {
                 for(j=1;j<NExy-1; j++)
                  {
                        
                                  
           Ex[1][i][j][0]=Ex[0][i][j][1]+tmp*(Ex[1][i][j][1]-Ex[0][i][j][0]); 

                   }
               }
             
             break;

         }   

   /* --------------------------------------------  */

        /*  Iteration Formulae of Ey Points  */

    switch(pModule->nSurf[1])
     {
       case 0:  /* It means that Down Plane is electric wall */

         for(i=1;i<NEyx-1;i++)
          {
            for(j=0;j<NEyy;j++)
              {
                Ey[1][i][j][0]=0.0;
              }
           }
           
            break;

          default:  /* G.Mur. First Order */

            for(i=1; i<NEyx-1; i++)
              {
                 for(j=0;j<NEyy; j++)
                  {
                        
                                  
           Ey[1][i][j][0]=Ey[0][i][j][1]+tmp*(Ey[1][i][j][1]-Ey[0][i][j][0]); 

                   }
               }
             
             break;

         }   

}   
      

void E_Up_Iter(MODULE *pModule, double ****Ex, double ****Ey)
{

    int i,j,N_k_End;
    int NExx,NExy,NExz;
    int NEyx,NEyy,NEyz;

    double tmp,er,vel;

    NExx=Nx;    NExy=Ny+1;  NExz=Nz+1;
    NEyx=Nx+1;  NEyy=Ny;    NEyz=Nz+1;

     er=pModule->pPermit[pModule->nLayerNum-1];
     vel=C/sqrt(er);  

     tmp=(vel*dt-dhz)/(vel*dt+dhz);
     N_k_End=pModule->pLayer[pModule->nLayerNum-1].z_end;

   /*  Iteration Formulae of Ex Points  */

    switch(pModule->nSurf[0])
     {
       case 0:  /* It means that Down Plane is electric wall */

         for(i=0;i<NExx;i++)
          {
            for(j=1;j<NExy-1;j++)
              {
                Ex[1][i][j][N_k_End]=0.0;
              }
           }
           
            break;

          default:  /* G.Mur. First Order */

            for(i=0; i<NExx; i++)
              {
                 for(j=1;j<NExy-1; j++)
                  {
                        
                                  
        Ex[1][i][j][N_k_End]=Ex[0][i][j][N_k_End-1]+
                   tmp*(Ex[1][i][j][N_k_End-1]-Ex[0][i][j][N_k_End]); 

                   }
               }
             
             break;

         }   

   /* --------------------------------------------  */

        /*  Iteration Formulae of Ey Points  */

    switch(pModule->nSurf[0])
     {
       case 0:  /* It means that Down Plane is electric wall */

         for(i=1;i<NEyx-1;i++)
          {
            for(j=0;j<NEyy;j++)
              {
                Ey[1][i][j][N_k_End]=0.0;
              }
           }
           
            break;

          default:  /* G.Mur. First Order */

            for(i=1; i<NEyx-1; i++)
              {
                 for(j=0;j<NEyy; j++)
                  {
                        
                                  
          Ey[1][i][j][N_k_End]=Ey[0][i][j][N_k_End-1]+
                        tmp*(Ey[1][i][j][N_k_End-1]-Ey[0][i][j][N_k_End]); 

                   }
               }
             
             break;

         }   

} 


void E_Front_Iter(MODULE *pModule, double ****Ey, double ****Ez)
{

  int l,k,j,tmp_k;
  int NEyx,NEyy,NEyz;
  int NEzx,NEzy,NEzz;

  double tmp,er,et,vel;

  NEyx=Nx+1; NEyy=Ny;   NEyz=Nz+1;
  NEzx=Nx+1; NEzy=Ny+1; NEzz=Nz;

  /* Iteration Formulae of Ey Points */

   for(l=0;l<pModule->nLayerNum; l++)
    {
      er=pModule->pPermit[l];
      vel=C/sqrt(er);  
      tmp=(vel*dt-dhx)/(vel*dt+dhx);     

      for(k=pModule->pLayer[l].z_start+1; k<pModule->pLayer[l].z_end; k++)
        {
          for(j=0;j<NEyy;j++)
           {
               
         Ey[1][NEyx-1][j][k]=Ey[0][NEyx-2][j][k]+
                        tmp*(Ey[1][NEyx-2][j][k]-Ey[0][NEyx-1][j][k]);

            }
        }
    }               

/* -------------------------------------------- */
    
    for(l=0;l<pModule->nLayerNum-1; l++)
    {
      et=0.5*(pModule->pPermit[l]+pModule->pPermit[l+1]);
      vel=C/sqrt(et);  
      tmp=(vel*dt-dhx)/(vel*dt+dhx);     

      tmp_k=pModule->pLayer[l+1].z_start;

      for(j=0;j<NEyy;j++)
        {
               
         Ey[1][NEyx-1][j][tmp_k]=Ey[0][NEyx-2][j][tmp_k]+
                        tmp*(Ey[1][NEyx-2][j][tmp_k]-Ey[0][NEyx-1][j][tmp_k]);

         }
        
    }  

 /*  Iteration Formulae of Ez Points */

 for(l=0; l<pModule->nLayerNum; l++)
   {
     er=pModule->pPermit[l];
     vel=C/sqrt(er);  
     tmp=(vel*dt-dhx)/(vel*dt+dhx); 

       for(k=pModule->pLayer[l].z_start; k<pModule->pLayer[l].z_end;k++)
        {
          for(j=1;j<NEzy-1;j++)
           {

          Ez[1][NEzx-1][j][k]=Ez[0][NEzx-2][j][k]+
                        tmp*(Ez[1][NEzx-2][j][k]-Ez[0][NEzx-1][j][k]);   

            }
         }
     }

         
}             
 

void E_Back_Iter(MODULE *pModule, double ****Ey, double ****Ez)
{
  int j,k,l,tmp_k;

  int NEyx,NEyy,NEyz;
  int NEzx,NEzy,NEzz;

  double tmp,er,et,vel;

  NEyx=Nx+1; NEyy=Ny;   NEyz=Nz+1;
  NEzx=Nx+1; NEzy=Ny+1; NEzz=Nz;

  /* Iteration Formulae of Ey Points */

   for(l=0;l<pModule->nLayerNum; l++)
    {
      er=pModule->pPermit[l];
      vel=C/sqrt(er);  
      tmp=(vel*dt-dhx)/(vel*dt+dhx);     

      for(k=pModule->pLayer[l].z_start+1; k<pModule->pLayer[l].z_end; k++)
        {
          for(j=0;j<NEyy;j++)
           {
               
         Ey[1][0][j][k]=Ey[0][1][j][k]+
                        tmp*(Ey[1][1][j][k]-Ey[0][0][j][k]);

            }
        }
    }               

/* -------------------------------------------- */
    
    for(l=0;l<pModule->nLayerNum-1; l++)
    {
      et=0.5*(pModule->pPermit[l]+pModule->pPermit[l+1]);
      vel=C/sqrt(et);  
      tmp=(vel*dt-dhx)/(vel*dt+dhx);     

      tmp_k=pModule->pLayer[l+1].z_start;

      for(j=0;j<NEyy;j++)
        {
               
         Ey[1][0][j][tmp_k]=Ey[0][1][j][tmp_k]+
                        tmp*(Ey[1][1][j][tmp_k]-Ey[0][0][j][tmp_k]);

         }
        
    }  

 /*  Iteration Formulae of Ez Points */

 for(l=0; l<pModule->nLayerNum; l++)
   {
     er=pModule->pPermit[l];
     vel=C/sqrt(er);  
     tmp=(vel*dt-dhx)/(vel*dt+dhx); 

       for(k=pModule->pLayer[l].z_start; k<pModule->pLayer[l].z_end;k++)
        {
          for(j=1;j<NEzy-1;j++)
           {

          Ez[1][0][j][k]=Ez[0][1][j][k]+
                        tmp*(Ez[1][1][j][k]-Ez[0][0][j][k]);   

            }
         }
     }

}              

void EMur_Left_Iter(MODULE *pModule,double ****Ex,double ****Ez)
{

   int i,k,l,tmp_k;

   int NExx,NExy,NExz;
   int NEzx,NEzy,NEzz;

   double tmp,vel,er,et;

   NExx=Nx;     NExy=Ny+1;    NExz=Nz+1;
   NEzx=Nx+1;   NEzy=Ny+1;    NEzz=Nz;

    er=0.5*(pModule->pPermit[0]+pModule->pPermit[1]);

    er=1.95; 

    vel=C/sqrt(er);  
    tmp=(vel*dt-dhy)/(vel*dt+dhy);


   /* Iteration Formulae of Ex Points  */

   for(l=0;l<pModule->nLayerNum; l++)
    {
      for(k=pModule->pLayer[l].z_start+1; k<pModule->pLayer[l].z_end; k++)
        {
          for(i=0;i<NExx;i++)
           {
               
         Ex[1][i][0][k]=Ex[0][i][1][k]+
                        tmp*(Ex[1][i][1][k]-Ex[0][i][0][k]);

            }
        }
    }               

/* -------------------------------------------- */
    
    for(l=0;l<pModule->nLayerNum-1; l++)
    {
      tmp_k=pModule->pLayer[l+1].z_start;

      for(i=0;i<NExx;i++)
        {
               
         Ex[1][i][0][tmp_k]=Ex[0][i][1][tmp_k]+
                        tmp*(Ex[1][i][1][tmp_k]-Ex[0][i][0][tmp_k]);

         }
        
    }  


   /*  Iteration Formulae of Ez Points  */
               
 for(l=0; l<pModule->nLayerNum; l++)
   {
       for(k=pModule->pLayer[l].z_start; k<pModule->pLayer[l].z_end;k++)
        {
          for(i=1;i<NEzx-1;i++)
           {

          Ez[1][i][0][k]=Ez[0][i][1][k]+
                        tmp*(Ez[1][i][1][k]-Ez[0][i][0][k]);   

            }
         }
     }

}

void EMur_Right_Iter(MODULE *pModule,double ****Ex,double ****Ez)
{

   int i,k,l,tmp_k;

   int NExx,NExy,NExz;
   int NEzx,NEzy,NEzz;

   double tmp,vel,er,et;

   NExx=Nx;     NExy=Ny+1;    NExz=Nz+1;
   NEzx=Nx+1;   NEzy=Ny+1;    NEzz=Nz;


   /* Iteration Formulae of Ex Points  */

      er=0.5*(pModule->pPermit[0]+pModule->pPermit[1]);
      er=1.95;

      vel=C/sqrt(er);  
      tmp=(vel*dt-dhy)/(vel*dt+dhy);


   for(l=0;l<pModule->nLayerNum; l++)
    {
      for(k=pModule->pLayer[l].z_start+1; k<pModule->pLayer[l].z_end; k++)
        {
          for(i=0;i<NExx;i++)
           {
               
         Ex[1][i][NExy-1][k]=Ex[0][i][NExy-2][k]+
                        tmp*(Ex[1][i][NExy-2][k]-Ex[0][i][NExy-1][k]);

            }
        }
    }               

/* -------------------------------------------- */
    
    for(l=0;l<pModule->nLayerNum-1; l++)
    {
      tmp_k=pModule->pLayer[l+1].z_start;

      for(i=0;i<NExx;i++)
        {
               
         Ex[1][i][NExy-1][tmp_k]=Ex[0][i][NExy-2][tmp_k]+
                        tmp*(Ex[1][i][NExy-2][tmp_k]-Ex[0][i][NExy-1][tmp_k]);

         }
        
    }  




   /*  Iteration Formulae of Ez Points  */
               
 for(l=0; l<pModule->nLayerNum; l++)
   {
       for(k=pModule->pLayer[l].z_start; k<pModule->pLayer[l].z_end;k++)
        {
          for(i=1;i<NEzx-1;i++)
           {

          Ez[1][i][NEzy-1][k]=Ez[0][i][NEzy-2][k]+
                        tmp*(Ez[1][i][NEzy-2][k]-Ez[0][i][NEzy-1][k]);   

            }
         }
     }

}





/*
   
2th Order Dispersive Absorbing Boundary Condition

*/

void E_Left_Iter(MODULE *pModule,double ****Ex,double ****Ez,double ***pLE_X, double ***pLE_Z)
{

   int i,k;

   int NExx,NExy,NExz;
   int NEzx,NEzy,NEzz;

   double gama1,gama2,v1,v2;

   v1=C/sqrt(erff1);  v2=C/sqrt(erff2);

   gama1=(dhy-v1*dt)/(v1*dt+dhy);  gama2=(dhy-v2*dt)/(v2*dt+dhy);


   NExx=Nx;     NExy=Ny+1;    NExz=Nz+1;
   NEzx=Nx+1;   NEzy=Ny+1;    NEzz=Nz;


   /* Iteration Formulae of Ex Points  */

    for(i=0; i<NExx; i++)
     {
       for(k=1;k<pModule->pLayer[pModule->nLayerNum-1].z_end; k++)
        {
          
      Ex[1][i][0][k]=2.0*Ex[0][i][1][k]-pLE_X[2][i][k]+(gama1+gama2)*
                     (Ex[0][i][0][k]-Ex[1][i][1][k]-pLE_X[1][i][k]
                      +Ex[0][i][2][k])-gama1*gama2*(pLE_X[0][i][k]-2.0*
                       Ex[0][i][1][k]+Ex[1][i][2][k]);
         }
      } 


   /*  Iteration Formulae of Ez Points  */
               
    for(i=1; i<NEzx-1; i++)
     {
       for(k=0;k<pModule->pLayer[pModule->nLayerNum-1].z_end; k++)
        {
          
      Ez[1][i][0][k]=2.0*Ez[0][i][1][k]-pLE_Z[2][i][k]+(gama1+gama2)*
                     (Ez[0][i][0][k]-Ez[1][i][1][k]-pLE_Z[1][i][k]
                      +Ez[0][i][2][k])-gama1*gama2*(pLE_Z[0][i][k]-2.0*
                       Ez[0][i][1][k]+Ez[1][i][2][k]);
         }
      } 

}


     
void E_Right_Iter(MODULE *pModule,double ****Ex,double ****Ez,double ***pRE_X, double ***pRE_Z)
{

   int i,k;

   int NExx,NExy,NExz;
   int NEzx,NEzy,NEzz;

   double v1,v2,gama1,gama2;

   v1=C/sqrt(erff1);  v2=C/sqrt(erff2);

   gama1=(dhy-v1*dt)/(v1*dt+dhy);  gama2=(dhy-v2*dt)/(v2*dt+dhy);


   NExx=Nx;     NExy=Ny+1;    NExz=Nz+1;
   NEzx=Nx+1;   NEzy=Ny+1;    NEzz=Nz;


   /* Iteration Formulae of Ex Points  */

    for(i=0; i<NExx; i++)
     {
       for(k=1;k<pModule->pLayer[pModule->nLayerNum-1].z_end; k++)
        {
          
      Ex[1][i][NExy-1][k]=2.0*Ex[0][i][NExy-2][k]-pRE_X[2][i][k]+(gama1+gama2)*
                     (Ex[0][i][NExy-1][k]-Ex[1][i][NExy-2][k]-pRE_X[1][i][k]
                      +Ex[0][i][NExy-3][k])-gama1*gama2*(pRE_X[0][i][k]-2.0*
                       Ex[0][i][NExy-2][k]+Ex[1][i][NExy-3][k]);
         }
      } 


   /*  Iteration Formulae of Ez Points  */
               
    for(i=1; i<NEzx-1; i++)
     {
       for(k=0;k<pModule->pLayer[pModule->nLayerNum-1].z_end; k++)
        {

       Ez[1][i][NEzy-1][k]=2.0*Ez[0][i][NEzy-2][k]-pRE_Z[2][i][k]+(gama1+gama2)*
                     (Ez[0][i][NEzy-1][k]-Ez[1][i][NEzy-2][k]-pRE_Z[1][i][k]
                      +Ez[0][i][NEzy-3][k])-gama1*gama2*(pRE_Z[0][i][k]-2.0*
                       Ez[0][i][NEzy-2][k]+Ez[1][i][NEzy-3][k]);    
      
         }
      } 


}   

void Save_B_Ex(MODULE *pModule,double ****Ex, double ***pLE_X,double ***pRE_X)
{
  int i,k,l;
  int NExx,NExy,NExz;

  NExx=Nx;   NExy=Ny+1;  NExz=Nz+1;



   for(l=0;l<3; l++)
    {
      for(i=0; i<NExx; i++)
       {
         for(k=1; k<NExz-1; k++)
          {
            pLE_X[l][i][k]=Ex[1][i][l][k];
            pRE_X[l][i][k]=Ex[1][i][NExy-1-l][k];

           }
        }
     }            
  
}


void Save_B_Ez(MODULE *pModule,double ****Ez, double ***pLE_Z,double ***pRE_Z)
{
  int i,k,l;
  int NEzx,NEzy,NEzz;

  NEzx=Nx+1;   NEzy=Ny+1;  NEzz=Nz;

   for(l=0;l<3; l++)
    {
      for(i=1; i<NEzx-1; i++)
       {
         for(k=0; k<NEzz; k++)
          {
            pLE_Z[l][i][k]=Ez[1][i][l][k];
            pRE_Z[l][i][k]=Ez[1][i][NEzy-1-l][k];

           }
        }
     }            
  
}


void Save_B_Hy(MODULE *pModule,double ****Hy, double ***pLH_Y,double ***pRH_Y)
{
  int i,k,l;
  int NHyx,NHyy,NHyz;

  NHyx=Nx;   NHyy=Ny+1;  NHyz=Nz;

   for(l=0;l<3; l++)
    {
      for(i=0; i<NHyx; i++)
       {
         for(k=0; k<NHyz; k++)
          {
            pLH_Y[l][i][k]=Hy[1][i][l][k];
            pRH_Y[l][i][k]=Hy[1][i][NHyy-1-l][k];

           }
        }
     }            
  
}