e_iter.c

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

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

        The Iteration Formulae of E Component

*/

# include "E_Iter.h"

void E_Inner_Iter(MODULE *pModule, double ****Ex,double ****Ey,double ****Ez,
                  double ****Hx,double ****Hy,double ****Hz,int time)
{
  int i,j,k,p;
  
  int NExx,NExy,NExz;
  int NEyx,NEyy,NEyz;
  int NEzx,NEzy,NEzz;

  double tmp1,tmp2,er;

  FILE *fp;

  fp=fopen("Ez.dat","a+");


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

  /* Iteration Formulae of Ex Inner Points  */

   for(p=0;p<pModule->nLayerNum;p++)
     {

       er=pModule->pPermit[p];

       for(k=pModule->pLayer[p].z_start+1;k<pModule->pLayer[p].z_end; k++)
        {
          for(i=0;i<NExx;i++)
           {
             for(j=1;j<NExy-1;j++)
              {
                tmp1=(Hz[1][i][j][k]-Hz[1][i][j-1][k])/dhy;
                tmp2=(Hy[1][i][j][k]-Hy[1][i][j][k-1])/dhz;

                Ex[1][i][j][k]=Ex[0][i][j][k]+dt*C/er*(tmp1-tmp2);
              
               }
            }
         }
     }


  /* Iteration Formulae of Ey Inner Points  */

   for(p=0;p<pModule->nLayerNum;p++)
     {

       er=pModule->pPermit[p];

       for(k=pModule->pLayer[p].z_start+1;k<pModule->pLayer[p].z_end; k++)
        {
          for(i=1;i<NEyx-1;i++)
           {
             for(j=0;j<NEyy;j++)
              {
                tmp1=(Hx[1][i][j][k]-Hx[1][i][j][k-1])/dhz;
                tmp2=(Hz[1][i][j][k]-Hz[1][i-1][j][k])/dhx;

                Ey[1][i][j][k]=Ey[0][i][j][k]+dt*C/er*(tmp1-tmp2);
              
               }
            }
         }
     }

                                       
    /* Iteration Formulae of Ez Inner Points */
      
   for(p=0;p<pModule->nLayerNum;p++)
     {
       er=pModule->pPermit[p];
 
       for(k=pModule->pLayer[p].z_start; k<pModule->pLayer[p].z_end;k++)
        {
          for(i=1;i<NEzx-1;i++)
           {
            for(j=1;j<NEzy-1;j++)
             {
                  tmp1=(Hy[1][i][j][k]-Hy[1][i-1][j][k])/dhx;
                  tmp2=(Hx[1][i][j][k]-Hx[1][i][j-1][k])/dhy;

                  Ez[1][i][j][k]=Ez[0][i][j][k]+dt*C/er*(tmp1-tmp2);

              }
           }
        }
     }

      fprintf(fp,"%d  %lf \n",time,Ez[1][38][60][3]);

      fclose(fp);  

}

       
void E_Surface_Iter(MODULE *pModule, double ****Ex, double ****Ey,
                    double ****Hx, double ****Hy, double ****Hz)
{
  int i,j,l,m,tmp_k;
  int NExx,NExy,NExz,NEyx,NEyy,NEyz;  
  double tmp1,tmp2,et;
  BOOL Flagc;
  POINT *p;

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

  p=(POINT*) malloc(sizeof(POINT));

  /* Iteration Formulae of Ex Surface Points */

  for(l=0;l<pModule->nLayerNum-1;l++)
   {

   tmp_k=pModule->pLayer[l+1].z_start; 
   et=0.5*(pModule->pPermit[l]+pModule->pPermit[l+1]); 

   if(pModule->pLayer[l].nCondPolyNum!=0)  /* Has conductor */  
    {  
     for(i=0;i<NExx;i++)
     { 
        for(j=1;j<NExy-1;j++)
        { 
              p->x=0.5+1.0*i;
              p->y=1.0*j;
        for(m=0;m<pModule->pLayer[l].nCondPolyNum;m++) /* scan all Conductors */
          {     
               Flagc=OnPolygon(&pModule->pLayer[l].pCondPoly[m],p);   
               if(Flagc==TRUE ) break;
            }
        if(Flagc==TRUE) /* It means that the point p is a conductor point */
           {
              
            Ex[1][i][j][tmp_k]=0.0; 

             
           }
         else
           {
                tmp1=(Hz[1][i][j][tmp_k]-Hz[1][i][j-1][tmp_k])/dhy;
                tmp2=(Hy[1][i][j][tmp_k]-Hy[1][i][j][tmp_k-1])/dhz;

                Ex[1][i][j][tmp_k]=Ex[0][i][j][tmp_k]+dt*C/et*(tmp1-tmp2);            
           }
        }
          
       }
       
      }
    else  /* No Conductors */
     {

       for(i=0;i<NExx;i++)
         { 
          for(j=1;j<NExy-1;j++)
          {
                tmp1=(Hz[1][i][j][tmp_k]-Hz[1][i][j-1][tmp_k])/dhy;
                tmp2=(Hy[1][i][j][tmp_k]-Hy[1][i][j][tmp_k-1])/dhz;

                Ex[1][i][j][tmp_k]=Ex[0][i][j][tmp_k]+dt*C/et*(tmp1-tmp2); 
            
 
           }
          }
       }
   }



  /* Iteration Formulae of Ey Surface Points */

  for(l=0;l<pModule->nLayerNum-1;l++)
   {

   tmp_k=pModule->pLayer[l+1].z_start; 
   et=0.5*(pModule->pPermit[l]+pModule->pPermit[l+1]); 

   if(pModule->pLayer[l].nCondPolyNum!=0)  /* Has conductor */  
    {  
     for(i=1;i<NEyx-1;i++)
     { 
        for(j=0;j<NEyy;j++)
        { 
              p->x=1.0*i;
              p->y=1.0*j+0.5;
        for(m=0;m<pModule->pLayer[l].nCondPolyNum;m++) /* scan all Conductors */
          {     
               Flagc=OnPolygon(&pModule->pLayer[l].pCondPoly[m],p);   
               if(Flagc==TRUE ) break;
            }
        if(Flagc==TRUE) /* It means that the point p is a conductor point */
           {
              
            Ey[1][i][j][tmp_k]=0.0; 

             
           }
         else
           {
                tmp1=(Hx[1][i][j][tmp_k]-Hx[1][i][j][tmp_k-1])/dhz;
                tmp2=(Hz[1][i][j][tmp_k]-Hz[1][i-1][j][tmp_k])/dhx;

                Ey[1][i][j][tmp_k]=Ey[0][i][j][tmp_k]+dt*C/et*(tmp1-tmp2);            
           }
        }
          
       }
       
      }
    else  /* No Conductors */
     {
       for(i=1;i<NEyx-1;i++)
         { 
          for(j=0;j<NEyy;j++)
          {
                tmp1=(Hx[1][i][j][tmp_k]-Hx[1][i][j][tmp_k-1])/dhz;
                tmp2=(Hz[1][i][j][tmp_k]-Hz[1][i-1][j][tmp_k])/dhx;

                Ey[1][i][j][tmp_k]=Ey[0][i][j][tmp_k]+dt*C/et*(tmp1-tmp2);
            
 
           }
          }
       }
   }

    free(p); 

}

void E_2_Surface_Iter(MODULE *pModule, double ****Ex, double ****Ey,
                    double ****Hx, double ****Hy, double ****Hz)
{
  int i,j,l,m,tmp_k;
  int NExx,NExy,NExz,NEyx,NEyy,NEyz;  
  double tmp1,tmp2,et;
  BOOL Flagc;
  POINT *p;
  POLYGON *pCondPoly;

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

  p=(POINT*) malloc(sizeof(POINT));  
 
  pCondPoly=POLYGON_malloc(1);  /* !! Note: Modify the Parameters  Future */

  pCondPoly->nVertex=4;   /*  The Microstrip Parameters  */  

  pCondPoly->pVertex=POINT_malloc(pCondPoly->nVertex);

  pCondPoly->pVertex[0].x=1.0*pModule->Excite.nFw_start;
  pCondPoly->pVertex[0].y=0.0;

  pCondPoly->pVertex[1].x=1.0*pModule->Excite.nFw_start;
  pCondPoly->pVertex[1].y=1.0*Ny;

  pCondPoly->pVertex[2].x=1.0*pModule->Excite.nFw_end;
  pCondPoly->pVertex[2].y=1.0*Ny;

  pCondPoly->pVertex[3].x=1.0*pModule->Excite.nFw_end;
  pCondPoly->pVertex[3].y=0.0;

  /* Iteration Formulae of Ex Surface Points */

   tmp_k=pModule->pLayer[1].z_start; 
   et=0.5*(pModule->pPermit[0]+pModule->pPermit[1]); 

     for(i=0;i<NExx;i++)
     { 
        for(j=1;j<NExy-1;j++)
        { 
              p->x=0.5+1.0*i;
              p->y=1.0*j;

         Flagc=OnPolygon(pCondPoly,p);   

        if(Flagc==TRUE) /* It means that the point p is a conductor point */
           {
              
            Ex[1][i][j][tmp_k]=0.0; 

             
           }
         else
           {
                tmp1=(Hz[1][i][j][tmp_k]-Hz[1][i][j-1][tmp_k])/dhy;
                tmp2=(Hy[1][i][j][tmp_k]-Hy[1][i][j][tmp_k-1])/dhz;

                Ex[1][i][j][tmp_k]=Ex[0][i][j][tmp_k]+dt*C/et*(tmp1-tmp2);            
           }
        }
          
     }
       

  /* Iteration Formulae of Ey Surface Points */


   tmp_k=pModule->pLayer[1].z_start; 
   et=0.5*(pModule->pPermit[0]+pModule->pPermit[1]); 

     for(i=1;i<NEyx-1;i++)
     { 
        for(j=0;j<NEyy;j++)
        { 
              p->x=1.0*i;
              p->y=1.0*j+0.5;

               Flagc=OnPolygon(pCondPoly,p);   

        if(Flagc==TRUE) /* It means that the point p is a conductor point */
           {
              
            Ey[1][i][j][tmp_k]=0.0; 

             
           }
         else
           {
                tmp1=(Hx[1][i][j][tmp_k]-Hx[1][i][j][tmp_k-1])/dhz;
                tmp2=(Hz[1][i][j][tmp_k]-Hz[1][i-1][j][tmp_k])/dhx;

                Ey[1][i][j][tmp_k]=Ey[0][i][j][tmp_k]+dt*C/et*(tmp1-tmp2);            
           }
        }
          
     }

    if (p!=NULL) {
                   free(p); 
                 }


    if( pCondPoly->pVertex !=NULL)
       {
         free(pCondPoly->pVertex); 
       } 
     

    if( pCondPoly !=NULL);
      { 
        free(pCondPoly);
      }
     
  
}