fdtd3d.c

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

   /* V=Vo*exp(-(tlength/2.)*(tlength/2.)/(timvar*timvar)); */
   /* Using (V/Vo)=0.001 to pick length of Gaussian voltage impulse */
   timvar=pow((tlength/2.)*(tlength/2.)/2.3,0.5);
   /* printf("denom of Eamp exponential: %f\n",timvar); */
   sigpr=1.0/(10.*fd1.dy); /* for underwater probe */

   /* ####################################################################### */
   hmul=fd1.dt*0.00079577471546/(fd1.dz);
   zcent=(fd1.zinter-fd1.anthei)/fd1.dz;
   zcere=zcent;
   xcent=fd1.xcen/fd1.dx;
   ycent=fd1.ycen/fd1.dy;
   xcere=fd1.xrcen/fd1.dx;
   if((fd1.anttype==2)||(fd1.anttype==3)||(fd1.anttype==0)) ycere=fd1.yrcen/fd1.dx;
   else ycere=ycent;
   if(fd1.anttype!=0) {
      if((xcere!=xas+midr1x)&&(fd1.anttype==3)) {
         fprintf(fparms,"xcere from parameter file = %d\n",xcere);
         xcere=xas+midr1x;
         fprintf(fparms,"adjusted xcere = %d\n",xcere);
      }
   }
   if(fd1.anttype>=25) {
      fprintf(fparms,"xcere from parameter file = %d\n",xcere);
      xcere=xas+midr1x;
      fprintf(fparms," xcere now= %d\n",xcere);
      xc=xcent;
      xc1=xc;
      yc1=ycent-1;
      xc2=xc;
      yc2=ycent+1;
      xcr1=xcere;
      ycr1=ycere-1;
      xcr2=xcere;
      ycr2=ycere+1;
      if(fd1.anttype==25) {
         tra1flag=1;
         rec1flag=1;
         tra2flag=0;
         rec2flag=0;
      }
      if(fd1.anttype==12345) {
         tra1flag=1;
         rec1flag=1;
         tra2flag=0;
         rec2flag=0;
         xcr1=xcere;
         ycr1=ycere-1;
         xcr2=xcere;
         ycr2=ycere+1;
      }
      if(fd1.anttype==12346) {
         tra1flag=1;
         rec1flag=1;
         tra2flag=0;
         rec2flag=0;
         xcr1=xcere-1;
         ycr1=ycere;
         xcr2=xcere+1;
         ycr2=ycere;
      }
   }
   if(fd1.anttype==2) {
      xc1=xcent;
      yc1=ycent-1;
      xc2=xcent;
      yc2=ycent+1;
      tra1flag=1;
      rec1flag=0;
      tra2flag=0;
      rec2flag=0;
   }
   if(fd1.anttype==3) {
      xc1=xcent;
      yc1=ycent-1;
      xc2=xcent;
      yc2=ycent+1;
      xcr1=xcere;
      ycr1=ycere-1;
      xcr2=xcere;
      ycr2=ycere+1;
      tra1flag=1;
      rec1flag=1;
      tra2flag=0;
      rec2flag=0;
   }
   if(fd1.anttype==0) {
      xcere=fd1.xrcen/fd1.dx;
      ycere=fd1.yrcen/fd1.dx;
      xc1=xcent;
      yc1=ycent-1;
      xc2=xcent;
      yc2=ycent+1;
      xcr2=xcent;
      ycr2=ycent+1;
      xcr1=xcent;
      ycr1=ycent+1;
      tra1flag=0;
      rec1flag=0;
      tra2flag=0;
      rec2flag=0;
   }
   if((xcent>xm-5)||(ycent>ym-5)||(xcent<5)||(ycent<5)) {
      fprintf(fparms,"XCENT OR YCENT ARE OUT OF BOUNDS - EXIT FLAG SET\n");
      exitflag=1;
   }
   if((xcere>xm-5)||(ycere>ym-5)||(xcere<5)||(ycere<5)) {
      fprintf(fparms,"XCERE OR YCERE ARE OUT OF BOUNDS - EXIT FLAG SET\n");
      exitflag=1;
   }
   if((interfac>zm-2)||(interfac<2)) {
      fprintf(fparms,"interfac OUT OF BOUNDS - EXIT FLAG SET\n");
      exitflag=1;
   }
   if((zcent>zm-2)||(zcent<2)||(zcere>zm-2)||(zcere<2)) {
      fprintf(fparms,"zcent or zcere OUT OF BOUNDS - EXIT FLAG SET\n");
      exitflag=1;
   }
   fprintf(fparms,"tra1flag,tra2flag,rec1flag,rec2flag = %d,%d,%d,%d\n",tra1flag,tra2flag,rec1flag,rec2flag);
   fprintf(fparms,"xc1,yc1 = %d,%d\n",xc1,yc1);
   fprintf(fparms,"xc2,yc2 = %d,%d\n",xc2,yc2);
   fprintf(fparms,"xc3,yc3 = %d,%d\n",xc3,yc3);
   fprintf(fparms,"xc4,yc4 = %d,%d\n",xc4,yc4);
   fprintf(fparms,"xcr1,ycr1 = %d,%d\n",xcr1,ycr1);
   fprintf(fparms,"xcr2,ycr2 = %d,%d\n",xcr2,ycr2);
   len=fd1.xlen/fd1.dx;
 
   zas=zcent-osudimz; /* FOR OSUANT */
   zcfin=zcent-2;
   zc=zcfin; 
   xc=xcent;
   yc=ycent-1;
   zcm=90;
   fprintf(fparms,"zcm = %d\n",zcm);
   if(zcm>CODIMABO-2) {
      printf("zcm > CODIMABO-2\n");
      zcm=CODIMABO-10;
      fprintf(fparms,"zcm now %d\n",zcm);
   }
   /* w=angular frequency for cw measurements */
   w=2.0*pi*0.20;
   fprintf(fparms,"xcent= %d\n",xcent); 
   fprintf(fparms,"ycent= %d\n",ycent);
   fprintf(fparms,"zcent= %d\n",zcent);
   fprintf(fparms,"interfac= %d\n",interfac);
   fprintf(fparms,"xcere= %d\n",xcere);
   fprintf(fparms,"ycere= %d\n",ycere);
   fprintf(fparms,"zcere= %d\n",zcere);
   fprintf(fparms,"ground conductivity: %f\n",fd1.zcond);
   fprintf(fparms,"ground permittivity: %f\n",fd1.zperm);
   prpt=interfac+fd1.prdepth/fd1.dz; 
   fprintf(fparms,"Z-LOCATION IN OF PROBE IN GRID (prpt): %d\n",prpt);
   if(prpt>zm-2) {
      fprintf(fparms,"prpt > zm-2\n");
      prpt=zm-10;
      fprintf(fparms,"prpt now %d\n",prpt);
   }

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

   CALCULATE THIN SHEET PARAMETERS FOR USE DURING EXECUTION

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

   if(fd1.sheet==1) {
      printf("THIN SHEET OPTION ACTIVATED\n");
      if(fd1.anttype==2) {
         shtxst=(fd1.xcen-fd1.shtxdim/2)/fd1.dx;
         shtyst=(fd1.ycen-fd1.shtydim/2)/fd1.dy;
      }
      else if(fd1.anttype>=25) {
         shtxst=(fd1.xcen)/fd1.dx-1;
         shtyst=(fd1.ycen-fd1.shtydim/2)/fd1.dy;
      }
      else if(fd1.anttype==3) {
         shtxst=((fd1.xcen+fd1.xrcen)/2.-fd1.shtxdim/2)/fd1.dx;
         shtyst=((fd1.ycen+fd1.yrcen)/2.-fd1.shtydim/2)/fd1.dy;
      }
      shtxlen=fd1.shtxdim/fd1.dx;
      shtylen=fd1.shtydim/fd1.dy;
      eave=(1.0-fd1.shtthk/fd1.dz)*1.0+fd1.shtthk*fd1.shtperm/fd1.dz;
      save=(fd1.shtthk*fd1.shtcond)/fd1.dz;
      eoldmulsht=(1.0-(fd1.shtcond*fd1.dt/(2.0*fd1.shtperm*eonan)))/(1.0+(fd1.shtcond*fd1.dt/(2.0*fd1.shtperm*eonan)));
      eoldmshave=(1.0-(save*fd1.dt/(2.0*eave*eonan)))/(1.0+(save*fd1.dt/(2.0*eave*eonan)));
      enmulsht  =fd1.dt/(fd1.shtperm*eonan*(1.0+fd1.shtcond*fd1.dt/(2.*fd1.shtperm*eonan)))*dxi;
      enmulshave=fd1.dt/(eave*eonan*(1.0+save*fd1.dt/(2.*eave*eonan)))*dxi;  
      zsh=interfac-(fd1.shthei/fd1.dz);
      if(zsh>zm-2) {
        fprintf(fparms,"zsh>zm-2 - EXIT FLAG SET \n");
        exitflag=1;
      }
      if((shtxst<4)||(shtxst+shtxlen>xm-4)) {
        fprintf(fparms,"PROBLEMS WITH shtxst - EXIT FLAG SET\n");
        exitflag=1;
      }
      if((shtyst<4)||(shtyst+shtylen>ym-4)) {
        fprintf(fparms,"PROBLEMS WITH shyxst - EXIT FLAG SET\n");
        exitflag=1;
      }
      if((zsh<4)||(zsh>zm-4)) {
        fprintf(fparms,"PROBLEMS WITH zsh - EXIT FLAG SET\n");
        exitflag=1;
      }
      fprintf(fparms,"SHEET STARTING X-GRID VALUE: %d\n",shtxst);
      fprintf(fparms,"SHEET X LENGTH (GRID PTS): %d\n",shtxlen);
      fprintf(fparms,"SHEET Y LENGTH (GRID PTS): %d\n",shtylen );
      fprintf(fparms,"SHEET STARTING Y-GRID VALUE: %d\n",shtyst);
      fprintf(fparms,"POSITION OF SHEET IN GRID (GRID PT.) %d\n",zsh); 
   }


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

   INITIALIZE PERMITTIVITY AND CONDUCTIVITY VALUES IN AIR/SUBSURFACE PORTIONS 
   OF GRID

   Ex, Ey, and Hz ARE TANGENT TO HORIZONTAL INTERFACES  

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

    if((fd1.zperm<0.01)||(fd1.zur<0.01)) {
       fprintf(fparms,"PERMITTIVITY/PERMEABILITY OF SUBSURFACE (fd1.zperm or fd1.zur) TOO LOW - EXITING\n");
       fclose(fparms);
       exit();
    }
    if(fd1.dt<0.00001) {
       fprintf(fparms,"TIMESTEP (fd1.dt) TOO SMALL - EXITING\n");
       fclose(fparms);
       exit();
    }
    for (z=0;z<DIMSIZZ;z++) {
       if (z<interfac) {
          perme[z]=1.0;
          sige[z]=0.0;
          permez[z]=1.0;
          sigez[z]=0.0;
	  ur[z]=1.0;
          hmulx[z]=fd1.dt*0.00079577471546/fd1.dx;
          hmuly[z]=fd1.dt*0.00079577471546/fd1.dy;
          hmulz[z]=fd1.dt*0.00079577471546/fd1.dz;
       }
       else if(z==interfac) {
         sige[z]=(0.0+fd1.zcond)/2.0; /* Ex and Ey continuous on interface */
         perme[z]=(1.0+fd1.zperm)/2.0; /* THEREFORE TAKE AVERAGE SIGMA AND EPSILON */
         sigez[z]=0.0; /* Ez does not lie on the interface */
         permez[z]=1.0;/* Ez does not lie on the interface */
         ur[z]=fd1.zur;/* Hz lies on the interface */
         hmulx[z]=fd1.dt*0.00079577471546/fd1.dx;
         hmuly[z]=fd1.dt*0.00079577471546/fd1.dy;
         hmulz[z]=fd1.dt*0.00079577471546/(fd1.dz*fd1.zur);
      }
      else {
        if(fd1.obj==10) {
              zo=interfac+(fd1.zpos)/fd1.dz;          
              if(z<zo) {
                   perme[z]=fd1.zperm;
                    sige[z]=fd1.zcond;
                    sigez[z]=fd1.zcond;
                    permez[z]=fd1.zperm;
                    ur[z]=fd1.zur;/* Hz lies on the interface */
                    hmulx[z]=fd1.dt*0.00079577471546/(fd1.dx*fd1.zur);
                    hmuly[z]=fd1.dt*0.00079577471546/(fd1.dy*fd1.zur);
                    hmulz[z]=fd1.dt*0.00079577471546/(fd1.dz*fd1.zur);
              }
              else if(z==zo) {
                   perme[z]=(fd1.targperm+fd1.zperm)/2.;
                    sige[z]=(fd1.targcond+fd1.zcond)/2.;
                    sigez[z]=fd1.zcond;
                    permez[z]=fd1.zperm;
                    ur[z]=fd1.zur;/* Hz lies on the interface */
                    hmulx[z]=fd1.dt*0.00079577471546/(fd1.dx*fd1.zur);
                    hmuly[z]=fd1.dt*0.00079577471546/(fd1.dy*fd1.zur);
                    hmulz[z]=fd1.dt*0.00079577471546/(fd1.dz*fd1.zur);
              }
              else if(z>zo) {
                   perme[z]=fd1.targperm;
                    sige[z]=fd1.targcond;
                    sigez[z]=fd1.targcond;
                    permez[z]=fd1.targperm;
                    ur[z]=fd1.zur;/* Hz lies on the interface */
                    hmulx[z]=fd1.dt*0.00079577471546/(fd1.dx*fd1.zur);
                    hmuly[z]=fd1.dt*0.00079577471546/(fd1.dy*fd1.zur);
                    hmulz[z]=fd1.dt*0.00079577471546/(fd1.dz*fd1.zur);
             }
         }
         else {                  
                   perme[z]=fd1.zperm;
                    sige[z]=fd1.zcond;
                    sigez[z]=fd1.zcond;
                    permez[z]=fd1.zperm;
                    ur[z]=fd1.zur;/* Hz lies on the interface */