[fdtd]fdtd源代码一个.txt

用FDTD算法解决电磁场问题

发信人: BlueDanube (blauen Donau), 信区: Electromagnetics
标  题: [FDTD]FDTD源代码一个
发信站: 飘渺水云间 (Sun Oct  3 19:40:23 2004), 转信

#include "math.h"
#include "stdio.h"
#include "stdlib.h"

#define Ke 202
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,dt,ddx,pi,epsc,sigma;
    int n, k, kc, kstart,epsilon, NSTEPS;
    float t,freq_in;
    float to,spread;
    double pulse;

    FILE *fp;

    kstart=100;
    epsilon=4;
    ddx=0.01;
    dt=ddx/(2*3e8);
    pi=3.1415926;
    epsc=8.85419e-12;
    sigma=0.04;
    kc=5;           /* center of the problem space*/
    to=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;
    fp=fopen("all.txt","w");


    for (k=0;k<=Ke;k++)
    {
        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)        /* NSTEPS is the number of the main loops executed
*/
    {
        printf ("\nNSTEPS-->");
        scanf ("%d",&NSTEPS);
        printf("%d\n",NSTEPS);
        n=0;

        for (n=1;n<=NSTEPS;n++)
        {
            t++;            /* T keeps the track of the total number */


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

            pulse=exp(-.5*(pow((to-t)/spread,2.0)));
            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];
            }

            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];

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

            for (k=1;k<Ke-1;k++)
            {
                fprintf (fp,"%3d %6.2f %6.2f\n",k,ex[k],dx[k]);
                printf ("%3d %6.2f %6.2f\n",k,ex[k],dx[k]);
            }
            printf ("T=%5.0f\n",t);
        }

    }
    fclose (fp);
}
--

Hi, Hsin, You know what you have to do now.

※ 来源:·飘渺水云间 freecity.cn·[FROM: BlueDanube]