电磁波时域有限差分方法（FDTD）-C语言-FDTD和激励源和MUR边界

本文主要是介绍电磁波时域有限差分方法（FDTD）-C语言-FDTD和激励源和MUR边界，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

参考：

葛德彪, 闫玉波. 电磁波时域有限差分方法[M]. 西安电子科技大学出版社, 2011.

Schneider J B. Understanding the Finite-Difference Time-Domain Method[J]. 2013.

最近要用FDTD做仿真，所以找了以上两本书，FDTD确实是有点难的。

想了想流程：

1：FDTD

2：激励源

3：吸收边界

4：散射

现在完成了FDTD和激励源，的TM部分，虽然都是抄书的，但对于一个不懂电磁波电磁场，数学极差的我来说，光是看懂代码就已经竭尽全力了。

何况还遇到了C语言宏定义的坑。

然后做了一部分的MUR吸收，效果不咋的，不知道就是这样的还是我代码有问题。

上代码，基本就是抄书的，但是简洁了一点。。

就不先做注释了，等到最后全部写完再做注释。

// Boundary.cpp: 定义控制台应用程序的入口点。
//#include "stdafx.h"
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define _CRT_SECURE_NO_WARNINGS#define ALLOC_2D(PNTR, NUMX, NUMY, TYPE)					\PNTR = (TYPE *)calloc((NUMX) * (NUMY), sizeof(TYPE));	\if (!PNTR) {												\perror("ALLOC_2D");									\fprintf(stderr, "Allocation failed for" #PNTR ". Terminating.\n");	\exit(-1);											\}
#define SizeX g->sizeX
#define SizeY g->sizeY
#define CDTDS g->cdtds
#define Time g->time
#define MaxTime g->maxTime#define Ez(M, N) g->ez[M * SizeY + N]
#define CA(M, N) g->ceze[M * SizeY + N]
#define CB(M, N) g->cezh[M * SizeY + N]#define Hx(M, N) g->hx[M * (SizeY) + N]
#define CP(M, N) g->chxh[M * (SizeY) + N]
#define CQ(M, N) g->chxe[M * (SizeY) + N]#define Hy(M, N) g->hy[M * SizeY + N]
//#define Chyh(M, N) g->chyh[M * SizeY + N]
//#define Chye(M, N) g->chye[M * SizeY + N]struct Grid {double *hx, *chxh, *chxe;double *hy, *chyh, *chye;double *ez, *ceze, *cezh;int sizeX, sizeY;int time, maxTime;int type;double cdtds;	//库朗数
};
typedef struct Grid Grid;void gridInit(Grid *g);
double driveSource(Grid *g, double time, double locX, double locY);
void writeFile(Grid *g);
void updateE(Grid *g);
void updateH(Grid *g);
void updateBoundary(Grid *g);
void EHInit(Grid *g);float c = 3e8;
float PI = 3.1415926535;
float MU0 = PI * 4e-7;
float EPS0 = (1 / (36 * PI)) * 1e-9;
float SIGMA0 = 0;
float SIGMAm0 = 0;
float DELTA = 5e-3;
float DELTAT = DELTA / (sqrt(2.0) * c);int main()
{int m, n;double source;Grid *g;ALLOC_2D(g, 100, 100, Grid);EHInit(g);gridInit(g);for (Time = 0; Time < 150; Time++){updateH(g);updateBoundary(g);updateE(g);source = driveSource(g, Time, 0, 0);Ez(SizeX / 2, SizeY / 2) = source;}writeFile(g);getchar();return 0;
}
void EHInit(Grid *g)
{int n, m;for (m = 0; m < SizeX; m++){for (n = 0; n < SizeY; n++){Hy(m, n) = 0;Hx(m, n) = 0;Ez(m, n) = 0;}}
}void gridInit(Grid *g)
{int m, n;SizeX = 100;SizeY = 100;MaxTime = 300;CDTDS = 1.0 / sqrt(2.0);ALLOC_2D(g->hx, SizeX, SizeY, double);ALLOC_2D(g->chxh, SizeX, SizeY, double);ALLOC_2D(g->chxe, SizeX, SizeY, double);ALLOC_2D(g->hy, SizeX, SizeY, double);ALLOC_2D(g->chyh, SizeX, SizeY, double);ALLOC_2D(g->chye, SizeX, SizeY, double);ALLOC_2D(g->ez, SizeX, SizeY, double);ALLOC_2D(g->ceze, SizeX, SizeY, double);ALLOC_2D(g->cezh, SizeX, SizeY, double);for (m = 0; m < SizeX; m++){for (n = 0; n < SizeY; n++){//Ez(m, n) = 0;CA(m, n) = (1 - 0.5 * SIGMA0 * DELTAT / MU0) / (1 + 0.5 * SIGMA0 * DELTAT / MU0);CB(m, n) = (DELTAT / EPS0) / (1 + (SIGMA0 * DELTAT / (2 * EPS0)));}}for (m = 0; m < SizeX; m++){for (n = 0; n < SizeY; n++){//Hx(m, n) = 0;CP(m, n) = (1 - (SIGMAm0 * DELTAT) / (2 * MU0)) / (1 + (SIGMAm0 * DELTAT) / (2 * MU0));CQ(m, n) = (DELTAT / MU0) / (1 + (SIGMAm0 * DELTAT / (2 * MU0)));}}printf("The CA is :%f\n", CA(50, 50));printf("The CB is :%f\n", CB(50, 50));printf("The CP is :%f\n", CP(50, 50));printf("The CQ is :%f\n", CQ(50, 50));
}double driveSource(Grid *g, double time, double locX, double locY)
{double result;double ricker_result;double ppw = 20;double cdtds = CDTDS;double w0 = 1000000000.0; //天线的中心频率double arg;result = time * time * sin(2 * 3.1415 * w0 * time) * exp(-0.93 * w0 * time);arg = 3.14159 * ((cdtds * time - 0.0) / ppw - 1.0);arg = arg * arg;ricker_result = (1.0 - 2.0 * arg) * exp(-arg);return ricker_result;
}void updateE(Grid *g)
{int m, n;for (m = 1; m < (SizeX - 1); m++){for (n = 1; n < (SizeY - 1); n++){Ez(m, n) = CA(m, n) * Ez(m, n) + CB(m, n) * ((Hy(m, n) - Hy((m - 1), n)) / DELTA - (Hx(m, n) - Hx(m, (n - 1))) / DELTA);}}
}void updateH(Grid *g)
{int m, n;for (m = 0; m < (SizeX); m++){for (n = 0; n < (SizeY - 1); n++){Hx(m, n) = CP(m, n) * Hx(m, n) - CQ(m, n) * (Ez(m, (n + 1)) - Ez(m, n)) / DELTA;}}for (m = 0; m < (SizeX - 1); m++){for (n = 0; n < (SizeY); n++){Hy(m, n) = CP(m, n) * Hy(m, n) + CQ(m, n) * (Ez((m + 1), n) - Ez(m, n)) / DELTA;}}
}void updateBoundary(Grid *g)
{int m, n;int firstX = 0, firstY = 0, lastX = 99, lastY = 99;Grid *g1;ALLOC_2D(g1, 100, 100, Grid);memcpy(g1, g, sizeof(Grid));gridInit(g1);//updateH(g1);updateE(g1);//修正Hy左边界m = firstX;for (n = 1; n < lastY; n++){Ez(m, n) = Ez((m + 1), n) + ((c * DELTAT - DELTA) / (c * DELTAT + DELTA)) * (g1->ez[(m+1) * SizeY + n] - Ez(m, n))//-((c * c * DELTAT * MU0) / (2 * (c * DELTAT + DELTA))) * (g1->hx[m*SizeY+(n+1)] - g1->hx[m*SizeY + (n - 1)] + g1->hx[(m+1)*SizeY + (n + 1)] - g1->hx[(m + 1)*SizeY + (n - 1)]);-((c * c * DELTAT * MU0) / (2 * (c * DELTAT + DELTA))) * (Hx(m, (n + 1)) - Hx(m, (n - 1)) + Hx((m + 1), (n + 1)) - Hx((m + 1), (n - 1)));}//修正Hy右边界m = lastX;for (n = 1; n < lastY; n++){Ez(m, n) = Ez((m - 1), n) + ((c * DELTAT - DELTA) / (c * DELTAT + DELTA)) * (g1->ez[(m - 1) * SizeY + n] - Ez(m, n))- ((c * c * DELTAT * MU0) / (2 * (c * DELTAT + DELTA))) * (Hx(m, (n + 1)) - Hx(m, (n - 1)) + Hx((m - 1), (n + 1)) - Hx((m - 1), (n - 1)));}/*//修正Hx下边界//修正Hx上边界*/
}void writeFile(Grid *g)
{int m, n;FILE *f = fopen("C:\\Users\\Administrator\\Desktop\\data.txt", "a");if (f == NULL) printf("open file error");for (n = SizeY - 1; n >= 0; n--){for (m = 0; m < SizeX; m++){fprintf(f, "%f,", Ez(m, n));}fprintf(f, "\n");}fclose(f);
}

然后用matlab对生成的东西成像：