811阅读
2回复

[求助]求助！Sullivan书中程序FD1D_2.2.c有个地方没懂

UID ：100451

注册：2012-10-15
登录：2012-11-30
发帖：3
等级：旁观者

0楼发表于: 2012-10-15 17:50:19

程序中在对Ex进行傅里叶变换之后那段对入射脉冲作傅里叶变换的if语句有什么意义?为什么乘的是ex[10]而不是加入脉冲的位置的ex[5]?T<100在这里又怎么解释？求各位指教，小弟初学，先谢过了。

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

# define KE 200                                  /* KE is the number of cells to be used */
void main ()
{
      float dx[KE], ex[KE], hy[KE], ix[KE];
      float ga[KE],gb[KE];

      int n,m,k,kc,ke,kstart,NSTEPS;
      float T,ddx,dt,epsz,epsilon,sigma;
      float t0,spread,pi,pulse;

      FILE *fp;

      float ex_low_m1,ex_low_m2,ex_high_m1,ex_high_m2;
      float real_pt[5][KE],imag_pt[5][KE];
      float freq[5],arg[5],ampn[5][KE],phasen[5][KE];
      float real_in[5],imag_in[5],amp_in[5],phase_in[5];
     // float mag[KE];

      pi = 3.14159;
      kc = KE/2;                                    /* Center of the problem space*/
      ddx = .01;                                 /* Cell size */
      dt = ddx/6e8;                              /* Time steps */
      epsz = 8.8e-12;

      printf (" %6.4f  %10.5e \n",ddx,dt);

      /*Initialize to free space  */

      for (k = 0; k <KE; k++)
      {
          ga[k] = 1.;
          gb[k] = 0.;
          ex[k] = 0;
          dx[k] = 0;
          hy[k] = 0;
          ix[k] = 0;
         // mag[k] = 0;

          for (m=0; m<=2; m++)
          {
              real_pt[m][k] = 0;                 /* Real and imaginary parts of Fourier Transform    */
              imag_pt[m][k] = 0;
              ampn[m][k] = 0;                    /* Amplitude and phase of the Fourier Transforms */
              phasen[m][k] =0;
          }
      }

      for ( m=0; m<=2; m++)
      {
          real_in[m] = 0;                        /* Fourier Trans of input pulse */
          imag_in[m] = 0;
      }

      ex_low_m1=0;
      ex_low_m2=0;
      ex_high_m1=0;
      ex_high_m2=0;

      /* Parameters for the Fourier Transform */

      freq[0] = 100.e6;
      freq[1] = 200.e6;
      freq[2] = 500.e6;

      for (n=0; n<=2; n++)
      {
          arg[n] = 2*pi*freq[n]*dt;
          printf ("%2d %6.2f %7.5f \n",n,freq[n]*1e-6,arg[n]);
      }

      printf ("Dielectric starts at -->");
      scanf ("%d",&kstart);
      printf ("Epsilon -->");
      scanf ("%f",&epsilon);
      printf ("Conductivity -->");
      scanf ("%f",&sigma);
      printf ("%d %6.2f %6.2f \n",kstart,epsilon,sigma);

      /* Initialize the parameters */

      for (k=kstart; k<KE; k++)
      {
          ga[k] = 1./(epsilon+sigma*dt/epsz);
          gb[k] = sigma*dt/epsz;
      }

      for (k=1; k<KE; k++)
      {
          printf ("%2d  ga[k]=%4.2f gb[k]=%4.2f\n",k,ga[k],gb[k]);
      }

      /* These parameters specify the input pulse */

      t0=50.0;                                   /* Center of the incident pulse */
      spread=10.0;                                 /* Width of the incident pulse */

      T=0;
      NSTEPS=1;
            while ( NSTEPS>0)
      {
      /* Main part of the program */

      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 */
                                                 /* of times the main loop is executed */
          /* Main FDTD Loop */

          /* Calculate the Dx field */

          for (k=1; k<KE; k++)
          {
              dx[k] = dx[k] +.5*(hy[k-1]-hy[k]);
          }

          /* Initialize with a Guassian pulse */

          pulse = exp(-.5*(pow( (t0-T)/spread,2.0) ));
          dx[5] = dx[5] + pulse;
          printf ("%5.1f %6.2f  %6.2f\n",T,pulse,dx[5]);

          /* Calculate Ex from Dx */

          for ( k=1; k<KE;k++)
          {
              ex[k] = ga[k]*(dx[k] - ix[k]);
              ix[k] = ix[k] + gb[k]*ex[k];
          }

           /* Calculate the Fourier transform of Ex. */
          for ( k=0;k<KE;k++)
          {
              for (m=0;m<=2;m++)
              {
                  real_pt[m][k] = real_pt[m][k] + cos(arg[m]*T)*ex[k];
                  imag_pt[m][k] = imag_pt[m][k] - sin(arg[m]*T)*ex[k];
              }
          }

          /* Fourier Transform of the input pulse */

          if (T<100)
          {
              for (m=0; m<=2; m++)
              {
                   real_in[m] = real_in[m] + cos (arg[m]*T)*ex[10];
                   imag_in[m] = imag_in[m] - sin (arg[m]*T)*ex[10];
              }
          }


          /*Absorbing Boundary Conditions*/

          ex[0]      =       ex_low_m2;
          ex_low_m2  =       ex_low_m1;
          ex_low_m1  =       ex[1];

          ex[KE-1]   =       ex_high_m2;
          ex_high_m2 =       ex_high_m1;
          ex_high_m1 =       ex[KE-2];

          /* Calculate the Hy field */

          for ( k=0; k<KE-1; k++)
          {
              hy[k] = hy[k] + .5*(ex[k]-ex[k+1]);
          }
      }

      /* End of the Main FDTD Loop */

      /* At the end of the calculation,print out
      the Ex and Hy fields */

      for (k=0; k<KE; k++)
      {
          printf ("%2d  %6.2f  %6.2f \n",k,dx[k],ex[k]);
      }

      /* Write the E field out ot a file "Ex" */
      fp=fopen( "Ex","w");
      for (k=0;k<KE;k++)
      {
          fprintf (fp," %6.3f \n",ex[k]);
      }
      fclose(fp);

      /* Calculate the amplitude and phase of each frequency */

      /* Amplitude and phase of the input pulse */
      for ( m=0; m<=2; m++)
      {
           amp_in[m] = sqrt (pow (imag_in[m],2) + pow(real_in[m],2));
           phase_in[m] = atan2 (imag_in[m],real_in[m]);

           printf ( "%d Input Pulse : %8.4f %8.4f %8.4f %7.2f \n",
           m,real_in[m],imag_in[m],amp_in[m],(180.0/pi)*phase_in[m]);

           for ( k=0; k<KE; k++)
           {
               ampn[m][k] = (1./amp_in[m])*sqrt( pow(real_pt[m][k],2.)
                        ..

未注册仅能浏览部分内容，查看全部内容及附件请先登录或注册

离线lawrencelsq

UID ：100451

注册：2012-10-15
登录：2012-11-30
发帖：3
等级：旁观者

1楼发表于: 2012-10-15 18:16:25

没人么求大神们指教啊

离线lx84

呵呵

UID ：30566

注册：2009-04-22
登录：2020-09-16
发帖：866
等级：积极交流六级

2楼发表于: 2012-10-15 18:55:46

在自由空间内在那个观察点傅里叶变换对结果都没有影响；关键事在你做傅里叶变换的时间段内必须有完整的时域波形；这个T<100；也就是说100步内在x=10观察点形成了完整波形，这不影响最后结果；当然你也可以用T<NSTEPS；只不过100步以后都是0，做福利也变换也没什么意义了。

继续！

发帖回复