1008阅读
3回复

[求助]一维FDTD源程序

I do

UID ：91162

注册：2012-03-29
登录：2013-03-15
发帖：19
等级：仿真新人

0楼发表于: 2012-04-22 15:24:49

关键词： FDTD

     这个程序可以运行，但是结果不准确，画出的图形走势不对麻烦大神提个意见

#include < math.h >
#include < stdio.h >
#include < stdlib.h >
#define KE 200
main()
{
    double dx[KE],ex[KE],hy[KE],ix[KE],ga[KE],gb[KE];
    double ex_low_m1,ex_low_m2,ex_high_m1,ex_high_m2;
    int n, k, kc, kstart,epsilon, NSTEPS;
    float t0,spread,dt,ddx,pi,epsc,sigma;
float T,freq_in;
    double pulse;
    FILE *fp;
    kstart=100;
    epsilon=4;
    ddx=0.01;           /* Cell size */
    dt=ddx/(2*6e8);     /* time steps */
    pi=3.1415926;
    epsc=8.8e-12;
    sigma=0.04;
    kc=KE/2;           /* center of the problem space*/
    t0=50.0;            /* center of the incident pulse */
    spread=12;          /* width of the incident pulse */
    T=0;
    NSTEPS=1;
    ex_low_m1=0;
    ex_low_m2=0;
    ex_high_m1=0;
    ex_high_m2=0;


    for (k=0;k<=KE;k++)      /*Initialize to free space */
    {
        ga[k]=1.0;
        gb[k]=0;
        dx[k]=0;
        ex[k]=0;
        hy[k]=0;
        ix[k]=0;
    }
    for (k=kstart;k<=KE;k++)
    {
        ga[k]=1/(epsilon+sigma*dt/epsc);
        gb[k]=sigma*dt/epsc;
    }

    while (NSTEPS>0)        /* main part of the program */
    {
        printf ("\nNSTEPS-->");
        scanf ("%d",&NSTEPS);
        printf("%d\n",NSTEPS);
        n=0;
        for (n=1;n<=NSTEPS;n++)
        {
            T++;

            for (k=1;k<=KE; k++)
            {
                dx[k]=dx[k]+0.5*(hy[k-1]-hy[k]);
            }
            pulse=exp(-.5*(pow((t0-T)/spread,2.0)));    /* put a Gassian pulse at the low end */
            dx[kc]=dx[kc]+pulse;
            printf("%5.1f %6.2f %6.2f \n",T,pulse,ex[5]);

            for (k=0;k<=KE; k++)
            {
                ex[k]=ga[k]*(dx[k]-ix[k]);
                ix[k]=ix[k]+gb[k]*ex[k];
            }
         /* 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; k++)
            {
                hy[k]=hy[k]+.5*(ex[k]-ex[k+1]);
            }
            hy[0]=0;
            hy[KE]=0;
        /* write the E field out to file " EX " */
   fp = fopen("EX", "w" ) ;
            for (k=1;k<KE-1;k++)
            {