fdtd_1.1_change.c

来自「一维FDTD」· C语言 代码 · 共 78 行

/* FDTD_1.1_change.c FDTD simulation in free space */

#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#define KE 200 /* KE is the number of cells to be used */

main(int argc, char *argv[])
{
	float ex[KE], hy[KE];
	int n, k, kc, ke, NSTEPS;
	int T, t0;
	float spread, pulse;
	FILE *fp;
	
	/* Initialize */
	kc = KE / 2;	/* Center of the problem space */
	t0 = 40; 	/* Center of the incident pulse */
	spread = 12;	/* Width of the incident pulse */
	if (!(argc == 3 && fopen("Ex", "r") > 0 && fopen("Hy", "r") > 0)) {
		T = 0;
		NSTEPS = 1;
		
		for (k = 0; k < KE; k++) {
			ex[k] = 0;
			hy[k] = 0;
		}
	} else {
		sscanf(argv[1], "%d", &T);
		sscanf(argv[2], "%d", &NSTEPS);
		
		fp = fopen("Ex", "r");
		for (k = 1; k < KE; k++)
			fscanf(fp, "%f",&ex[k]);
		fclose(fp);
		ex[0] = 0;
		
		fp = fopen("Hy", "r");
		for (k = 0; k < KE - 1; k++)
			fscanf(fp, "%f", &hy[k]);
		fclose(fp);
		hy[KE - 1] = 0;
	}
	printf("OK\tT = %d\tNSTEPS = %d\n", T, NSTEPS);
	for (n = 1; n <= NSTEPS; n++) {
		T = T + 1;	/* T keeps track of the total number of times the main loop is executed */
		
		/* Main FDTD Loop */
		for (k = 1; k < KE; k++)
			ex[k] = ex[k] + 0.5*(hy[k-1] - hy[k]);
		
		/* Put a Gaussian pulse in the middle */
		pulse = exp(-0.5*(pow((t0 - T) / spread, 2.0)));
		ex[kc] = pulse;
		printf("%d	%f\n", t0 - T, ex[kc]);
		
		/* Calculate the Hy field */
		for (k = 0; k < KE - 1; k++)
			hy[k] = hy[k] + 0.5*(ex[k] - ex[k+1]);
	}
	
	/* End of the Main FDTD Loop */
	
	/* Write the E field out to a file "Ex" */
	fp = fopen("Ex", "w");
	for (k = 1; k < KE; k++)
		fprintf(fp, "%f\n", ex[k]);
	fclose(fp);
	
	fp = fopen("Hy", "w");
	for (k = 0; k < KE - 1; k++)
		fprintf(fp, "%f\n", hy[k]);
	fclose(fp);
	
	printf("T = %d\tNSTEPS = %d\n", T, NSTEPS);
}