fdtd_1.5.c

本程序是使用加入边界条件的有限时域差分法计算模拟电磁波在光学薄膜中的传输

/* FDTD_1.5.c FDTD simulation									*/
/* Simulation of a sinusoidal wave hitting a dielectric medium	*/
/* 1D FDTD simulation of a lossy dielectric medium				*/
/* Absorbing Boundary Condition added							*/
/* method : run "FDTD_1D.exe fileName"							*/

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#define pi 3.14159265358979323846

int main(int argc, char *argv[])
{
	float *Ex, *Hy, *ca, *cb;
	float Ex_low_m1, Ex_low_m2, Ex_high_m1, Ex_high_m2;
	float ddx, dt, freq_in, pulse;
	int n, k, T, NSTEPS, KE;
	FILE *fp;
	
	/* Initialize */
	if (argc == 1) {
		printf("error: no file to open!");
		exit(0);
	} else if ((fp = fopen(argv[1], "r")) <= 0) {
		printf("error: cannot open the file named \"%s\"", argv[1]);
		exit(0);
	} else {
		fscanf(fp, "%4*c%d", &T);
		fscanf(fp, "%10*c%d", &NSTEPS);
		fscanf(fp, "%6*c%d", &KE);
		fscanf(fp, "%7*c%f", &ddx);
		fscanf(fp, "%6*c%e", &dt);
		fscanf(fp, "%11*c%f", &freq_in);
		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);
		fscanf(fp, "%2*c");
		
		ca = (float *)malloc(KE * sizeof(float));
		for (k = 0; k < KE; k++)
			fscanf(fp, "%f", &ca[k]);
		
		cb = (float *)malloc(KE * sizeof(float));
		for (k = 0; k < KE; k++)
			fscanf(fp, "%f", &cb[k]);
			
		Ex = (float *)malloc((KE + 1) * sizeof(float));
		for (k = 0; k < KE + 1; k++)
			fscanf(fp, "%f", &Ex[k]);
		
		Hy = (float *)malloc(KE * sizeof(float));
		for (k = 0; k < KE; k++)
			fscanf(fp, "%f", &Hy[k]);
	}
	fclose(fp);
	
	/* Main FDTD Loop */
	for (n = 1; n <= NSTEPS; n++) {
		T = T + 1; /* T keeps track of the total number of times the main loop is executed */
		
		/* Calculate the Ex field */
		for (k = 1; k < KE; k++)
			Ex[k] = ca[k] * Ex[k] + cb[k] * (Hy[k-1] - Hy[k]);
		
		/* Put a Sine wave */
		pulse = sin(2 * pi * freq_in * dt * T);
		Ex[5] = Ex[5] + pulse;
		
		/* Absorbing Boundary Conditions */
		Ex[0] = Ex_low_m2;
		Ex_low_m2 = Ex_low_m1;
		Ex_low_m1 = Ex[1];
		Ex[KE] = Ex_high_m2;
		Ex_high_m2 = Ex_high_m1;
		Ex_high_m1 = Ex[KE-1];
		
		/* Calculate the Hy field */
		for (k = 0; k < KE; k++)
			Hy[k] = Hy[k] + 0.5 * (Ex[k] - Ex[k + 1]);
	}
	
	/* End of the Main FDTD Loop */
	fp = fopen(argv[1], "w");
	fprintf(fp, "T = %d;\n", T);
	fprintf(fp, "NSTEPS = %d;\n", NSTEPS);
	fprintf(fp, "KE = %d;\n", KE);
	fprintf(fp, "ddx = %f;\n", ddx);
	fprintf(fp, "dt = %e;\n", dt);
	fprintf(fp, "freq_in = %f;\n", freq_in);
	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 = 0; k < KE; k++)
		fprintf(fp, "%f\t", ca[k]);
	fprintf(fp, "\n");
	for (k = 0; k < KE; k++)
		fprintf(fp, "%f\t", cb[k]);
	fprintf(fp, "\n");
	for (k = 0; k < KE + 1; k++)
		fprintf(fp, "%f\t", Ex[k]);
	fprintf(fp, "\n");
	for (k = 0; k < KE; k++)
		fprintf(fp, "%f\t", Hy[k]);
	fprintf(fp, "\n");
	fclose(fp);
	free(ca);free(cb);free(Ex);free(Hy);
	
	printf("T=%d NSTEPS=%d KE=%d\n", T, NSTEPS, KE);
	
	return 0;
}