main.c

这是一个用FDTD计算远场的c程序

            for(k=kb+1;k<KE;k++){
                kzh=k-kb-1;
                 curl_h=(hy[i][j][k]-hy[i-1][j][k]-hx[i][j][k]+hx[i][j-1][k]);
                  idzh[i][j][kzh]=idzh[i][j][kzh]+curl_h;
             dz[i][j][k] = gi3[i]*gj3[j]*dz[i][j][k]+gi2[i]*gj2[j]*0.5*(curl_h+
             gk1[k]*idzh[i][j][kzh]);
            }
          }
      }


    /* Calculate the E from D field */
      for(k=1;k<KE-1;k++){
          for(j=1;j<JE-1;j++){
            for(i=1;i<IE-1;i++){
               ez[i][j][k]=gaz[i][j][k]*dz[i][j][k];
               ex[i][j][k]=gax[i][j][k]*dx[i][j][k];
               ey[i][j][k]=gay[i][j][k]*dy[i][j][k];
            }
          }
      }


    fprintf(fptime,"%7.3f",T*0.017);
    fprintf(fpt1,"%8.4f ",ez[ic][jmetal+2][kc]);
    fprintf(fpt2,"%8.4f ",ez[ic][jmetal+10][kc]);
    fprintf(fpt3,"%8.4f ",ez[ic][jmetal+15][kc]);
    fprintf(fpt4,"%8.4f ",ez[ic][jmetal+40][kc]);
    fprintf(fpt5,"%8.4f ",ez[ic][jmetal+50][kc]);
    fprintf(fpt6,"%8.4f ",ez[ic][jmetal+60][kc]);
    fprintf(fpt7,"%8.4f ",ez[ic][jmetal+70][kc]);





  /* Calculate the incident field */

  for(j=0;j<JE-1;j++){
      hx_inc[j]=hx_inc[j]+0.5*(ez_inc[j]-ez_inc[j+1]);
  }

    /* Remember : part of the PML is E=0 at the edge */

     for(j=0;j<JE-1;j++){
         for(k=0;k<KE-1;k++){
             ez[0][j][k]=0.0;
             ez[IE-1][j][k]=0.0;
         }
     }
     for(i=0;i<IE-1;i++){
         for(k=0;k<KE;k++){
             ez[i][0][k]=0.0;
             ez[i][JE-1][k]=0.0;
         }
     }
     for(i=0;i<IE-1;i++){
         for(j=0;j<JE-1;j++){
             ez[i][j][0]=0.0;
             ez[i][j][KE-1]=0.0;
         }
     }
      /* Calculate the Hx field */
      for(i=0;i<ia;i++){
          for(j=0;j<JE-1;j++){
            for(k=0;k<KE-1;k++){
              curl_e=(ey[i][j][k+1]-ey[i][j][k]-ez[i][j+1][k]+ez[i][j][k]);
                ihxl[i][j][k]=ihxl[i][j][k]+curl_e;
               hx[i][j][k] = hx[i][j][k]+fk2[k]*0.5*(
                             curl_e+fi1[i]*ihxl[i][j][k]);
            }
          }
      }
      for(i=ia;i<=ib;i++){
          for(j=0;j<JE-1;j++){
            for(k=0;k<KE-1;k++){
            curl_e=(ey[i][j][k+1]-ey[i][j][k]-ez[i][j+1][k]+ez[i][j][k]);
                hx[i][j][k] = fj3[j]*fk3[k]*hx[i][j][k]+fj2[j]*fk2[k]*0.5*
                             curl_e;
            }
          }
      }
       for(i=ib+1;i<IE;i++){
          ixh=i-ib-1;
          for(j=0;j<JE-1;j++){
            for(k=0;k<KE-1;k++){
                curl_e=(ey[i][j][k+1]-ey[i][j][k]-ez[i][j+1][k]+ez[i][j][k]);
                 ihxh[ixh][j][k]=ihxh[ixh][j][k]+curl_e;
               hx[i][j][k] = fj3[j]*fk3[k]*hx[i][j][k]+fj2[j]*fk2[k]*0.5*(
                             curl_e+fi1[i]*ihxh[ixh][j][k]);
            }
          }
       }

        /*  Incident Hx values  */
     for(i=ia;i<=ib;i++){
         for(k=ka;k<=kb;k++){
         hx[i][ja-1][k]=hx[i][ja-1][k]+0.5*ez_inc[ja];
     }
     }

        /* Calculate the Hy field */
    for(i=0;i<IE-1;i++){
          for(j=0;j<ja;j++){
            for(k=0;k<KE-1;k++){
                curl_e=(ez[i+1][j][k]-ez[i][j][k]-ex[i][j][k+1]+ex[i][j][k]);
        ihyl[i][j][k]=ihyl[i][j][k]+curl_e;
               hy[i][j][k] = fi3[i]*fk3[k]*hy[i][j][k]+fi2[i]*fk3[k]*0.5*(
                             curl_e+fj1[j]*ihyl[i][j][k]);
            }
          }
    }
     */
    for(i=0;i<IE-1;i++){
          for(j=jmetal+1;j<=jb;j++){
            for(k=0;k<KE-1;k++){
                 curl_e=(ez[i+1][j][k]-ez[i][j][k]-ex[i][j][k+1]+ex[i][j][k]);
              hy[i][j][k] = fi3[i]*fk3[k]*hy[i][j][k]+fi2[i]*fk3[k]*0.5*
                             curl_e;
            }
          }
    }

     /* Incident Hy values */
       for(j=ja;j<=jb;j++){
            for(k=ka;k<=kb;k++){

           hy[ib][j][k]  =hy[ib][j][k]+0.5*ez_inc[j];
       }
    }
        /* Calculate the Hz field */
       for(i=0;i<IE-1;i++){
          for(j=jmetal+1;j<JE-1;j++){
            for(k=0;k<ka;k++){
                curl_e=(ex[i][j+1][k]-ex[i][j][k]-ey[i+1][j][k]+ey[i][j][k]);
                ihzl[i][j][k]=ihzl[i][j][k]+curl_e;
               hz[i][j][k] = fi3[i]*fj3[j]*hz[i][j][k]+fi2[i]*fj2[j]*0.5*(
                             curl_e+fk1[k]*ihzl[i][j][k]);

            }
          }
       }



if(T==10){
 fp=fopen("Ez10.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ez[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }
if(T==20){
 fp=fopen("Ez20.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ez[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }
if(T==30){
 fp=fopen("Ez30.dat","w");
  for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ez[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }
  if(T==40){
 fp=fopen("Ez40.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ez[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }

 if(T==50){
 fp=fopen("Ez50.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ez[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }


 if(T==60){
 fp=fopen("Ey60.dat","w");
  for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ey[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }
 if(T==70){
 fp=fopen("Ex70.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ey[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }
 if(T==80){
 fp=fopen("Ey80.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ey[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }

 if(T==81){
 fp=fopen("Ez81.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ez[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }

 if(T==82){
 fp=fopen("Ex82.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ex[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }

 if(T==83){
 fp=fopen("Hz83.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",hz[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }

 if(T==90){
 fp=fopen("Ey90.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ey[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }
 if(T==91){
 fp=fopen("Ex91.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ex[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }
 if(T==92){
 fp=fopen("Hz92.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",hz[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }



 if(T==121){
 fp=fopen("Ez121.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ez[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }

  if(T==150){
 fp=fopen("Ez150.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ez[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }
  if(T==151){
 fp=fopen("Ex151.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ex[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }

  if(T==160){
 fp=fopen("Ez160.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",ez[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }
  if(T==161){
 fp=fopen("Hz161.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%21.10f",hz[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
 }


       /*   End of the Main FDTD Loop  */



}