二维FDTD正演C代码

#include“stdio.h“
#include“math.h“
#include“malloc.h“
#define c 3e8
#define dt 0.1
/****************  TM wave   ****************/
/******* CA=0.1    CB=0.1   CP=0.1   CQ=0.1***********/
int main（）
{
	printf（“please input compute times:/n“）;
	const float X_0 = -100 X_1 = 100 Y_0 = -100 Y_1 = 100;
	const float dx = 5 dy = 5;
	int i jnk;
	float M;
	M = （c*dt - sqrt（dx*dx + dy * dy）） / （c*dt + sqrt（dx*dx + dy * dy））;
	scanf（“%d“ &n）;
	float *N_Hx = （float *）malloc（sizeof（float）*40*40）;
	float *N_Hy = （float *）malloc（sizeof（float） * 40*40）;
	float *N_Ez = （float *）malloc（sizeof（float） * 41*41）;
	float *O_Hx = （float *）malloc（sizeof（float） * 40 * 40）;
	float *O_Hy = （float *）malloc（sizeof（float） * 40 * 40）;
	float *O_Ez = （float *）malloc（sizeof（float） * 41 * 41）;
	float *n_1Ez = （float *）malloc（sizeof（float） * 39 * 8）;
	for （i = 0; i < 40 * 40; i++）
	{
		*（O_Hx + i） = 0;
		*（O_Hy + i） = 0;
		*（O_Ez + i） = 0;
	}
	for （i = 0; i < 41+40; i++）
	{
		*（O_Ez + i + 40 * 40） = 0;
	}
	for （i = 0; i < 39 * 8; i++）
	{
		*（n_1Ez + i） = 0;
	}
	/*************** （00） as sin（0.001pit） *****************/
	for （k = 1; k <= n; k++）
	{
		for （i = 1; i < 40; i++）
			for （j = 1; j < 40; j++）
			{
				*（N_Ez + i * 41 + j） = *（O_Ez + i * 41 + j）*0.1 + 0.1*（（*（O_Hy + i * 40 + j ） - *（O_Hy + i * 40 + j-1）） / dx - （*（O_Hx + （i） * 40 + j-1） - *（O_Hy + （i-1） * 40 + j-1）） / dy）;
			}
		*（N_Ez） = *（O_Ez+42）+M * （*（N_Ez+42）-*（O_Ez））;
		*（N_Ez + 40） = *（O_Ez + 41+39） + M * （*（N_Ez + 41+39） - *（O_Ez +40））;
		*（N_Ez + 41*40） = *（O_Ez + 41*39+1） + M * （*（N_Ez + 41*39+1） - *（O_Ez + 41*40））;
		*（N_Ez + 41*41-1） = *（O_Ez + 41*40-2） + M * （*（N_Ez + 41*40-2） - *（O_Ez + 41*41-1））;
		for （i = 0; i < 39; i++）
		{
			*（N_Ez + i + 1） = -*（n_1Ez + i + 39） + （c*dt - dx） / （c*dt + dx）*（*（N_Ez + i + 41 + 1） + *（n_1Ez + i）） + 2 * dx / （c*dt + dx）*（*（O_Ez + i + 1） + *（O_Ez + i + 1 + 41）） + dx * c*dt*c*dt / （2 * dy*dy*（c*dt + dx））*（（*（O_Ez + i） - 2 * （*（O_Ez + i + 1）） + *（O_Ez + i + 2）） + *（O_Ez + i + 41） - 2 * （*（O_Ez + i + 1 + 41）） + *（O_Ez + i + 2 + 41））;
			*（N_Ez + i + 1+41*40） = -*（n_1Ez + i + 39*3）