yanshe2.c

这是一个用FDTD方法计算二维TM波通过金属屏的衍射程序。很少见。不可多得。

/* fd2d_3.3c 2D TM program with plane source */


#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define IE 70          /*IE is the number of cells to the X lable*/
#define JE 70          /*IE is the number of cells to the Y lable*/
#define pi 3.1415926
#define epsz 8.85419e-12

main()
{
    float ga[IE][IE],dz[IE][IE],ez[IE][JE],hy[IE][JE],hx[IE][JE];
    float ez_inc_low_m1=0.,ez_inc_low_m2=0.,ez_inc_high_m1=0.,
          ez_inc_high_m2=0.;
    float ez_inc[JE],hx_inc[JE];
    int i,j,n,jc,ic,NSTEPS,npml;
    int ia,ib,ja,jb;
    float T,ddx,dt,epsilon,omega;
    float xn,xxn,xnum,xd,curl_e;
    float t0,spread,pulse;
    float gi2[IE],gi3[IE];
    float gj2[JE],gj3[JE];
    float fi1[IE],fi2[IE],fi3[IE];
    float fj1[JE],fj2[JE],fj3[JE];
    float ihx[IE][JE],ihy[IE][JE];
    FILE *fp,*fopen(),*fpt;

    jc=JE/2;
    ic=IE/2;
      ia=14;   /* Total/scattered field boundaries  */
      ib=IE-ia-1;
      ja=14;
      jb=JE-ja-10;
    ddx=.01;             /* Set the cell size to 1cm */
    dt=ddx/(2*3e8);      /* Calculate the time stem */
    /*Initialize to free space */
    for(j=0;j<JE;j++){
        ez_inc[j]=0.;
        hx_inc[j]=0.;
    }
    for(j=0;j<JE;j++){
        printf("%2d ",j);
        for(i=0;i<IE;i++){
          dz[i][j]=0.;
          ez[i][j]=0.;
          hx[i][j]=0.;
          hy[i][j]=1.;
          ihx[i][j]=0.;
          ihy[i][j]=0.;
          ga[i][j]=1.0;
          printf("%5.2f ",ga[i][j]);
        }
        printf( "\n");
    }

     /* Calculate the PML parameters */

     for(i=0;i<IE;i++){
          gi2[i]=1.0;
          gi3[i]=1.0;
          fi1[i]=0.0;
          fi2[i]=1.0;
          fi3[i]=1.0;
      }
       for(j=0;j<JE;j++){
          gj2[j]=1.0;
          gj3[j]=1.0;
          fj1[j]=0.0;
          fj2[j]=1.0;
          fj3[j]=1.0;
       }
       printf("Number of PML cells--->");
       scanf("%d",&npml);

       for(j=0;j<=npml;j++){
           xnum =npml-i;
           xd   =npml;
           xxn  =xnum/xd;
           xn   =0.33*pow(xxn,3.0);
           printf(" %d  %7.4f  %7.4f  \n",i,xxn,xn);
              gi2[i]      =1.0/(1.0+xn);
              gi2[IE-1-j] =1.0/(1.0+xn);
              gi3[i]      =(1.0-xn)/(1.0+xn);
              gi3[IE-j-1] =(1.0-xn)/(1.0+xn);
            xxn=(xnum-0.5)/xd;
            xn =0.33*pow(xxn,3.0);
              fi1[i]     =xn;
              fi1[IE-2-i]=xn;
              fi2[i]     =1.0/(1.0+xn);
              fi2[IE-2-j]=1.0/(1.0+xn);
              fi3[i]     =(1.0-xn)/(1.0+xn);
              fi3[IE-2-j] =(1.0-xn)/(1.0+xn);
       }

        for(j=0;j<=npml;j++){
           xnum =npml-j;
           xd   =npml;
           xxn  =xnum/xd;
           xn   =0.333*pow(xxn,3.0);
           printf(" %d  %7.4f  %7.4f  \n",i,xxn,xn);
              gj2[j]      =1.0/(1.0+xn);
              gj2[JE-1-i] =1.0/(1.0+xn);
              gj3[j]      =(1.0-xn)/(1.0+xn);
              gj3[JE-j-1] =(1.0-xn)/(1.0+xn);
            xxn=(xnum-0.5)/xd;
            xn =0.33*pow(xxn,3.0);
              fj1[j]     =xn;
              fj1[JE-2-i]=xn;
              fj2[i]     =1.0/(1.0+xn);
              fj2[JE-2-i]=1.0/(1.0+xn);
              fj3[j]     =(1.0-xn)/(1.0+xn);
              fj3[JE-2-i] =(1.0-xn)/(1.0+xn);
       }

       printf("gi +fi \n");
       for(i=0;i<IE;i++){
           printf("%2d   %5.2f  %5.2f \n",i,gi2[i],gi3[i]),
           printf("%5.2f   %5.2f  %5.2f \n",fi1[i],fi2[i],fi3[i]);
       }

       printf("gi +fj \n");
       for(j=0;j<JE;j++){
           printf("%2d   %5.2f  %5.2f \n",j,gj2[j],gj3[j]),
           printf("%5.2f   %5.2f  %5.2f \n",i,fj1[j],fj2[j],fj3[j]);
       }
    for(i=0;i<IE;i++){
        j=40;
        ga[i][j]=0.0;
    }
 for(i=30;i<=35;i++){
     j=40;
      ga[i][j]=1.0;
    }
     t0=20;              /*Center of the incident pulse*/
     spread=8.0;          /*Width of the incident pulse*/
     T=0;
     NSTEPS=1;
   fpt=fopen("Ez_time1(yuanchang20,120).dat","w");

     while (NSTEPS>0) {
       printf("NSTEPS-->");  /*NSTEPS is the number of times the*/
       scanf("%d",&NSTEPS);  /* main loop has executed*/
       printf("%d\n",NSTEPS);

     for (n=1; n<=NSTEPS; n++)
    {

      T=T+1;           /*T keeps track of the total number*/
      /*    main FDTD Loop    */
       for(j=1;j<JE;j++){
            ez_inc[j]=ez_inc[j]+0.5*(hx_inc[j-1]-hx_inc[j]);
       }
        ez_inc[40]=0.0;
       /* ABC for the incident buffer*/
       ez_inc[1]       =ez_inc_low_m2;
       ez_inc_low_m2   =ez_inc_low_m1;
       ez_inc_low_m1   =ez_inc[1];

       ez_inc[JE-1]    =ez_inc_high_m2;
       ez_inc_high_m2  =ez_inc_high_m1;
       ez_inc_high_m1  =ez_inc[IE-2];


     /*   Calculate the Dz field  */
        for(j=1;j<JE;j++){
             for(i=1;i<IE;i++){
         dz[i][j]=gi3[i]*gj3[j]*dz[i][j]+gi2[i]*gj2[j]*0.5*(hy[i][j]
                  -hy[i-1][j]-hx[i][j]+ex[i][j-1]);
            }
        }

       /*  Source  */

   /*   pulse = sin(2*pi*1500*1e6*dt*T);*/
     pulse = exp(-0.5*(pow((t0-T)/spread,2.0)));
      ez_inc[3]= ez_inc[3]+pulse;

      /* Incident Dz values  */
      for(i=ia;i<=ib;i++){
          dz[i][ja]=dz[i][ja]+0.5*hx_inc[ja-1];
          dz[i][jb]=dz[i][jb]-0.5*hx_inc[jb];
      }

      /* Calculate the Ez field */
       for(j=0;j<JE;j++){
             for(i=0;i<IE;i++){
           ez[i][j]=ga[i][j]*dz[i][j];
             }
       }


      fprintf(fpt,"%8.4f ",ez[9+0][120]);
      /* Set the Ez edges to 0,as part of the PML */
      for(j=0;j<JE-1;j++){
        ez[0][j]=0.0;
        ez[jE-1][j]=0.0;
      }
       for(i=0;i<IE-1;i++){
        ez[i][0]=0.0;
        ez[i][JE-1]=0.0;
      }
     for(j=0;j<JE;j++){
         hx_inc[j]=hx_inc[j]+0.5*(ez_inc[j]-ez_inc[j+1]);
     }

        /* Calculate the Hx field */
         for(j=0;j<JE-1;j++){
             for(i=0;i<IE;i++){
               curl_e =ez[i][j] -ez[i][j+1];
               ihx[i][j]=ihx[i][j] + fi1[i]*curl_e;
               hx[i][j]=fj3[j]*hx[i][j]+fj2[j]*0.5*(curl_e + ihx[i][j]);
             }
         }


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

    /* Calculate the Hy field */
         for(j=0;j<IE-1;j++){
             for(i=0;i<IE-1;i++){
           curl_e =ez[i+1][j] -ez[i][j];
           ihy[i][j]=ihy[i][j]+fj1[j]*curl_e;
           hy[i][j]= fi3[j]*hy[i][j]+fi2[j]*0.5*(curl_e + ihy[i][j]);
             }
       }

       /* Incident Hy values */
       for(j=ja;j<=jb;j++){
           hy[ia-1][j] =hy[ia-1][j]-0.5*ez_inc[j];
           hy[ib][j]  =hy[ib][j]+0.5*ez_inc[j];
       }

       /* write the E field out to a field "Ez"  */
    if(T==10){   fp=fopen("Ez10.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 10.5f",ez[i][j]);
        }
           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,"% 10.5f",ez[i][j]);
        }
           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,"% 10.5f",ez[i][j]);
        }
           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,"% 10.5f",ez[i][j]);
        }
           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,"% 10.5f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }

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

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

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

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

    if(T==100){   fp=fopen("Ez100.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 10.5f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
    if(T==110){   fp=fopen("Ez110.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 10.5f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }

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


    if(T==130){   fp=fopen("Ez130.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 10.5f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
     if(T==140){   fp=fopen("Ez140.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 10.5f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
    if(T==115){   fp=fopen("Ez115.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 18.7f",ez[i][j]);
        }
           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,"% 18.7f",ez[i][j]);
        }
           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,"% 18.7f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
     if(T==170){   fp=fopen("Ez170.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 18.7f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
      if(T==180){   fp=fopen("Ez180.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 18.7f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
      if(T==190){   fp=fopen("Ez190.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 18.7f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
      if(T==200){   fp=fopen("Ez200.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 18.7f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
       if(T==250){   fp=fopen("Ez250.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 18.7f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
    if(T==270){   fp=fopen("Ez270.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 18.7f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
    if(T==290){   fp=fopen("Ez290.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 18.7f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
    if(T==310){   fp=fopen("Ez310.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 18.7f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
    if(T==330){   fp=fopen("Ez330.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 18.7f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
    if(T==350){   fp=fopen("Ez350.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 18.7f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
    if(T==370){   fp=fopen("Ez370.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 18.7f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
    if(T==400){   fp=fopen("Ez400.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"% 18.7f",ez[i][j]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
    }
    }

       /*   End of the Main FDTD Loop  */



     }
}