1252阅读
4回复

[求助]这是网上下载的一维FDTD的C语言模拟程序，哪位高手能帮我看下这几行代码

UID ：65748

注册：2010-09-03
登录：2019-04-05
发帖：122
等级：仿真二级

0楼发表于: 2010-09-12 22:20:10

我在VS2008中调试，结果应该存放在fp文件中吧，但是我找不到这个文件，请帮我分析一下代码是否正确。
另，我刚刚学习用C语言写FDTD，请给点思路吧。非常感谢！

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define NUM_of_ZAXIS 400
float gauss_pulse(float T,float t0,float spread);
int file_save(float* data,char* filename);
void main()
{
    float ex[NUM_of_ZAXIS],hy[NUM_of_ZAXIS];
    float obj_parameters[NUM_of_ZAXIS][4]; /*模型参数设置 */
    float CA[NUM_of_ZAXIS],CB[NUM_of_ZAXIS],CP[NUM_of_ZAXIS],CQ[NUM_of_ZAXIS];
    float ca,cb,cp,cq;        /*真空时的参数*/
    float var_ca,var_cb,var_cp,var_cq;        /* */
    float var_border,ex_low_m1,ex_low_m2,ex_high_m1,ex_high_m2,    ex_low_s1,ex_low_s2,ex_high_s1,ex_high_s2;    /*边界吸收参数*/
    float Epsilon,Mu,Pi,C;          /*介电系数 Epsilon 0，磁导系数 Mu 0*/
    float rel_epsz,rel_mu;   /*相对介电系数，磁导系数*/
    float e_sigma,h_sigma;   /*电导率，磁导率*/

    float dt,ddz;
    float source,T;
    int pos_driv_source;
    int i,n,Nsteps;
    int e_low,e_high;

    FILE* fp;

    /*初始化各个变量 */
    Pi=3.14159;
    Epsilon=8.85e-12;
    Mu=4*Pi*1e-7;

    C=pow( (float)(1/(Epsilon*Mu)),(float)(0.5));

    ddz=0.01;
    dt=ddz/(2*C);
    T=0;

                    // 采用葛德彪书中一维FDTD公式中的参数
    e_sigma=0;
    h_sigma=0;
    var_ca=0.5*e_sigma*dt/Epsilon;
    var_cb=dt/Epsilon;
    var_cp=0.5*h_sigma*dt/Mu;
    var_cq=dt/Mu;

    ca=(1-var_ca)/(1+var_ca);
    cb=var_cb/(1+var_ca);
    cp=(1-var_cp)/(1+var_cp);
    cq=var_cq/(1+var_cp);

    var_border=(C*dt-ddz)/(C*dt+ddz);


                   // 边界处用来存数据的临时变量
    ex_low_m1=0;
    ex_low_m2=0;
    ex_high_m1=0;
    ex_high_m2=0;
    ex_low_s1=0;
    ex_low_s2=0;
    ex_high_s1=0;
    ex_high_s2=0;

                    // 上，下边界位置
    e_low=20;
    e_high=380;

    pos_driv_source=200;


    for(i=0;i<NUM_of_ZAXIS;i++)
    {
        ex=0;
        hy=0;
        obj_parameters[0]=1;
        obj_parameters[1]=1;
        obj_parameters[2]=0;
        obj_parameters[3]=0;
        CA=ca;
        CB=cb;
        CP=cp;
        CQ=cq;
    }


    // 开始主循环

    Nsteps=1;

    while(Nsteps>0)
    {
        printf("Nsteps-->");
        scanf("%d",&Nsteps);

        for(n=0;n<Nsteps;n++)
        {
            T+=1;

            // 计算 Ex
            for(i=e_low+1;i<e_high;i++)
            {
                                                                 // 采用葛德彪书中的一维FDTD 公式
                ex=CA*ex-CB*(hy-hy[i-1])/ddz;
            }

            if(T<200){

            ex[pos_driv_source]+=gauss_pulse(T,40,12);

            }


            //Mur 边界吸收

            ex_low_m2=ex_low_m1;
            ex_low_m1=ex[e_low];

            ex_low_s2=ex_low_s1;
            ex_low_s1=ex[e_low+1];

            ex_high_m2=ex_high_m1;
            ex_high_m1=ex[e_high];

            ex_high_s2=ex_high_s1;
            ex_high_s1=ex[e_high-1];


                             ex[e_low]=ex_low_s2+var_border*(ex[e_low+1]-ex_low_m2);
                             ex[e_high]=ex_high_s2+var_border*(ex[e_high-1]-ex_high_m2);

                                          //                行波时延法
        //    ex[e_high]=0.5*(ex_high_m2+ex_high_s2);
        //    ex[e_low]=0.5*(ex_low_m2+ex_low_s2);

            // 计算 Hy

            for(i=e_low;i<e_high;i++)
            {
                hy=CP*hy-CQ*(ex[i+1]-ex)/ddz;
            }

        }

        // 存储数据

        fp=fopen("Ex","w");
        for(i=0;i<NUM_of_ZAXIS;i++) <szD"p|K

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

离线wjxylazypig

UID ：65748

注册：2010-09-03
登录：2019-04-05
发帖：122
等级：仿真二级

1楼发表于: 2010-09-13 09:06:52

版主帮我看看吧。

离线lx84

呵呵

UID ：30566

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

2楼发表于: 2010-10-05 00:16:38

楼主的C有待加强啊。

继续！

离线wjxylazypig

UID ：65748

注册：2010-09-03
登录：2019-04-05
发帖：122
等级：仿真二级

3楼发表于: 2010-10-05 11:19:06

回 2楼(lx84) 的帖子

期待能帮助解决问题的回复。感谢。

离线lx84

呵呵

UID ：30566

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

4楼发表于: 2010-10-08 18:21:02

数据，保存在同文件夹，Ex.Hy文件中，由于没加后缀名称，看不到正常图标，打开时选记事本就行了。

继续！

发帖回复