fdtd_1d.c

本程序是使用有限时域差分法模拟电磁波在光学薄膜中的传输

/* FDTD_1.2.c	1D FDTD simulation in free space */
/* Absorbing Boundary Condition added */
/* method: FDTD_1D.exe [output]*/

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

int main(int argc, char *argv[])
{
	float *Ex, *Hy;
	float Ex_low_m1, Ex_low_m2, Ex_high_m1, Ex_high_m2;
	int n, k, kc, NSTEPS;
	int T, t0, KE;
	float spread, pulse;
	FILE *fp;
	
	/* Initialize */
	for (k = 1; k < argc; k++) {
		if (!(strcmp(argv[k], "output"))) {
			fp = fopen("output", "r");
			if (fp <= 0) {
				printf("error: cannot open the output file\n");
				exit(0);
			} else {
				fscanf(fp, "%5*c%d", &KE);
				fscanf(fp, "%6*c%d", &kc);
				fscanf(fp, "%6*c%d", &t0);
				fscanf(fp, "%10*c%f", &spread);
				fscanf(fp, "%5*c%d", &T);
				fscanf(fp, "%10*c%d", &NSTEPS);
				fscanf(fp, "%13*c%f", &Ex_low_m1);
				fscanf(fp, "%13*c%f", &Ex_low_m2);
				fscanf(fp, "%14*c%f", &Ex_high_m1);
				fscanf(fp, "%14*c%f", &Ex_high_m2);
				
				Ex = (float *)malloc(KE*sizeof(float));
				Hy = (float *)malloc(KE*sizeof(float));
				fscanf(fp, "%2*c");
				for (k = 1; k < KE; k++) {
					fscanf(fp, "%f", &Ex[k]);
				}
				Ex[0] = 0.0;
				for (k = 0; k < KE - 1; k++) {
					fscanf(fp, "%f", &Hy[k]);
				}
				Hy[KE] = 0.0;
			}
			fclose(fp);
			break;
		}
	}
	
	if (k == 1 || k == argc) {
		KE = 200;		/* KE is the number of cells to be used */
		kc = KE / 2;	/* Center of the problem space */
		t0 = 40;		/* Center of the incident pulse */
		spread = 12;	/* Width of the incident pulse */
		T = 0;			/* Initial step */
		NSTEPS = 1;		/* run step */
        Ex_low_m1 = 0; Ex_low_m2 = 0; Ex_high_m1 = 0; Ex_high_m2 = 0;
		
		Ex = (float *)malloc(KE*sizeof(float));
		Hy = (float *)malloc(KE*sizeof(float));
		for (k = 0; k < KE; k++) {
			Ex[k] = 0.0;	/* Initial Ex field */
			Hy[k] = 0.0;	/* Initial Hy field */
		}
	}
	

	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 */
		/* Calculate the Ex field */
		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] = Ex[kc] + pulse;
		
		/* Absorbing 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 - 1; k++)
			Hy[k] = Hy[k] + 0.5 * (Ex[k] - Ex[k + 1]);
	}
	
	/* End of the Main FDTD Loop */
	
	/* Write the Ex,Hy field and settings out to a file "output" */
	fp = fopen("output", "w");
	fprintf(fp, "KE = %d;\n", KE);
	fprintf(fp, "kc = %d;\n", kc);
	fprintf(fp, "t0 = %d;\n", t0);
	fprintf(fp, "spread = %f;\n", spread);
	fprintf(fp, "T = %d;\n", T);
	fprintf(fp, "NSTEPS = %d;\n", NSTEPS);
	fprintf(fp, "Ex_low_m1 = %f;\n", Ex_low_m1);
	fprintf(fp, "Ex_low_m2 = %f;\n", Ex_low_m2);
	fprintf(fp, "Ex_high_m1 = %f;\n", Ex_high_m1);
	fprintf(fp, "Ex_high_m2 = %f;\n", Ex_high_m2);
	fprintf(fp, "\n");
	for (k = 1; k < KE; k++)
		fprintf(fp, "%f\t", Ex[k]);
	fprintf(fp, "\n");
	for (k = 0; k < KE - 1; k++)
		fprintf(fp, "%f\t", Hy[k]);
	fprintf(fp, "\n");
	fclose(fp);
	
	printf("KE=%d\tkc=%d\tt0=%d\tspread=%f\tT=%d\tNSTEPS=%d\n", KE, kc, t0, spread, T, NSTEPS);
	free(Ex);
	free(Hy);
}