1543阅读
8回复

一段一维FDTD程序的问题

UID ：43237

注册：2009-10-07
登录：2014-01-15
发帖：55
等级：仿真一级

0楼发表于: 2009-11-14 22:09:41

这是一个叫DENIS的人写的一本书上面的例子，但是我在我的机子运行的结果就是不对，不知道哪里出了错，我是新手，望多多指教。

/*1D FDTD simulation in free space*/
#include<math.h>
#include<stdlib.h>
#include<stdio.h>

#define KE 200

main()
{
    float ex[KE],  hy[KE];
    int n, k, kc, ke, NSTEPS;
    float T;
    float t0, spread, pulse;
    FILE *fp; //, *fopen();

    /*Initialize*/
    for (k=1; k<KE; k++)
    {
        ex[k] = 0;
        hy[k] = 0;
    }

    kc = KE/2;           //Center of the problem space
    t0 = 40.0;              //Center of the incident pulse
    spread = 12;          //Width of the incident pulse
    T = 0;
    NSTEPS = 1;

    while(NSTEPS > 0)
    {
        printf("NSTEPS -->");            //NSTEPS is the number of times the main loop has executed
        scanf("%d", &NSTEPS);
        printf("%d\n", NSTEPS);
        n = 0;
        for (n=1; n<=NSTEPS; n++)
        {
            T = T+1;
            //Main FDTD Loop
            /*Calculate the Ex field*/
            for ( k=1; k < KE; k++)
            {
                ex[k] = ex[k]+ 0.5*( hy[k-1] - hy[k] );
            }
            //Put a Gaussian pulse in the middle
            pulse = exp(-0.5*(pow((t0-T)/spread, 2.0 )));
            ex[kc] = pulse;
            printf( " %5.1f  %6.2f \n", t0 - T, ex[kc] );
            //Calculate the Hy field
            for (k=0; k<KE-1; k++)
            {
                hy[k] = hy[k] +0.5* (ex[k] -  ex[k+1] );
            }
        }
        //End of the Main loop
        //At the end of the calculation ,print out the Ex and Hy fields
        for (k=1; k<=KE; k++)
        {
            printf("%3d   %6f   %6.2f\n", k, ex[k], hy[k]);
        }

        //Write the E field out to a file "Ex"
        fp = fopen("Ex.dat","w");
        for (k=1; k<=KE; k++)
        {
            fprintf(fp, "%d %6.2f \n",k, ex[k]);
        }
        fclose(fp);

        //Write the H field out to a file "Hy"
     &nb ..

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

离线learnerlcy

UID ：43237

注册：2009-10-07
登录：2014-01-15
发帖：55
等级：仿真一级

1楼发表于: 2009-11-14 22:12:12

这是计算结果画的图

图片:untitled.jpg

离线learnerlcy

UID ：43237

注册：2009-10-07
登录：2014-01-15
发帖：55
等级：仿真一级

2楼发表于: 2009-11-15 00:04:48

这个结果似乎是发散的啊~~不知道怎么弄才对啊？

离线learnerlcy

UID ：43237

注册：2009-10-07
登录：2014-01-15
发帖：55
等级：仿真一级

3楼发表于: 2009-11-15 12:08:58

怎么没人回答我，急呀~这次我把C语言的代码文件给附上~希望哪位高手指教一二啊！

附件:

fdtd_1.1a.rar (1 K) 下载次数:10

离线kerving16

UID ：48824

注册：2009-12-14
登录：2009-12-14
发帖：3
等级：旁观者

4楼发表于: 2009-12-14 21:07:46

研究研究~~~~~~~~~~~

离线gwzhao

方恨少

UID ：17098

注册：2008-08-24
登录：2019-01-09
发帖：1374
等级：荣誉管理员

5楼发表于: 2009-12-14 22:16:01

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

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

逆流而上

离线gwzhao

方恨少

UID ：17098

注册：2008-08-24
登录：2019-01-09
发帖：1374
等级：荣誉管理员

6楼发表于: 2009-12-14 22:16:26

你改一下，试试看行不行。

逆流而上

离线alang

UID ：41410

注册：2009-09-11
登录：2010-03-02
发帖：77
等级：仿真一级

7楼发表于: 2010-02-07 10:49:17

xiexielouzhufenxianga

离线learnerlcy

UID ：43237

注册：2009-10-07
登录：2014-01-15
发帖：55
等级：仿真一级

8楼发表于: 2010-03-17 23:22:43

回 5楼(gwzhao) 的帖子

谢谢啊，最近在忙别的事情，好久没来了~这个程序我已经解决掉了！谢谢你！
/*1D FDTD simulation in free space*/
#include<math.h>
#include<stdlib.h>
#include<stdio.h>

#define KE 200

main()
{
    double ex[KE],  hy[KE];
    int n, k, kc, ke, NSTEPS;
    double T;
    double t0, spread, pulse;
    double ex_low_m1, ex_low_m2, ex_high_m1, ex_high_m2;
    FILE *fp;

    //Initialize
    for (k=0; k<KE; k++)
    {
        ex[k] = 0;
        hy[k] = 0;
    }

    kc = KE/2;           //Center of the problem space
    t0 = 40.0;              //Center of the incident pulse
    spread = 25;          //Width of the incident pulse
    T = 0;
    NSTEPS = 50;
    //Initialize boundary conditions' variable
    ex_low_m1 = 0;
    ex_low_m2 = 0;
    ex_high_m1 = 0;
    ex_high_m2 = 0;
    //Initialization finished

    while(NSTEPS > 0)
    {
        printf("NSTEPS -->");            //NSTEPS is the number of times the main loop has executed
        scanf("%d", &NSTEPS);
        printf("%d\n", NSTEPS);
        n = 0;
        for (n=1; n<=NSTEPS; n++)
        {
            T = T+1;
            //Main FDTD Loop
            /*Calculate the Ex field*/
            for ( k=1; k < KE; k++)
            {
                ex[k] = ex[k]+ 0.5*( hy[k-1] - hy[k] );
            }
            //Put a Gaussian pulse in the middle
            pulse = exp(-3.0*(pow((t0-T)/spread, 2.0 ))); //exp()函数前面的系数决定了振幅的大小
            ex[kc] = pulse;
            printf( " %5.1f  %6.2f \n", t0 - T, ex[kc] );

            //Absorbing Boundary Conditions
            //第一次改变：加入边界条件。
            //这里使用的是Toflove吸收边界条件
            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];
            //Absorbing Boundary Conditions is finished

            //Calculate the Hy field
            for (k=0; k<KE-1; k++)
            {
                hy[k] = hy[k] + 0.5*( ex[k] - ex[k+1] );
            }
        }
        //End of the Main loop
        //At the end of the calculation ,print out the Ex and Hy fields
        for (k=0; k<KE; k++)
        {
            printf("%3d   %6f   %6.2f\n", k, ex[k], hy[k]);
        }

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

        //Write the H field out to a file "Hy"
        fp = fopen("Hy.dat","w");
        for (k=0; k<KE; k++)
        {
            fprintf(fp, "%d %6.2f \n", k, hy[k]);
        }
        fclose(fp);
        printf("T= %5.0f\n",T);
    }
}

发帖回复