fdtd_1.3.c

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

/* FDTD_1.3.c	1D FDTD simulation in free space		*/
/* Absorbing Boundary Condition added					*/
/* method :												*/
/* First, run "FDTD_1D.exe" to initialize the variable	*/
/* 		  and make a file named "output" 				*/
/* Second, change the variable in the "output" file		*/
/* 		   and run "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, kstart,NSTEPS;
	int T, t0, KE;
	float spread, pulse, epsilon, *cb;
	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);
				fscanf(fp, "%10*c%d", &kstart);
				fscanf(fp, "%11*c%f", &epsilon);
				fscanf(fp, "%2*c");
				
				cb = (float *)malloc(KE*sizeof(float));
				for (k = 0; k < KE; k++)
					fscanf(fp, "%f", &cb[k]);
				
				Ex = (float *)malloc(KE*sizeof(float));
				Hy = (float *)malloc(KE*sizeof(float));
				
				for (k = 0; k < KE + 1; k++) {
					fscanf(fp, "%f", &Ex[k]);
				}
				for (k = 0; k < KE; k++) {
					fscanf(fp, "%f", &Hy[k]);
				}
			}
			fclose(fp);
			break;
		}
	}
	if (argc == 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 this time					*/
		Ex_low_m1 = 0;		/* set Absorbing Boundary Conditions	*/
		Ex_low_m2 = 0;
		Ex_high_m1 = 0;
		Ex_high_m2 = 0;
		
		/* Initialize to free space */
		kstart = kc;
		epsilon = 4;
		cb = (float *)malloc(KE * sizeof(float));
		for (k = 1; k< KE; k++)
			cb[k] = 0.5;
		cb[0] = 0.5;
		for ( k = kstart; k < KE; k++)
			cb[k] = 0.5/epsilon;
		
		Ex = (float *)malloc((KE + 1) * sizeof(float));
		Hy = (float *)malloc(KE * sizeof(float));
		for (k = 0; k < KE + 1; k++)
			Ex[k] = 0.0;	/* Initial Ex field */
		for (k = 0; k < KE; k++)
			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] + cb[k] * (Hy[k-1] - Hy[k]);
		
		/* Put a Gaussian pulse in the middle */
		pulse = exp(-0.5 * (pow((t0 - T) / spread, 2.0)));
		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 */
	
	/* Write the Ex,Hy field and settings out to a file "output" */
	/*fp = fopen("output", "w");
	fprintf(fp, "########################################################################\n");
	fprintf(fp, "#KE =\n\t%d\n", KE);
	fprintf(fp, "#kc =\n\t%d\n", kc);
	fprintf(fp, "#t0 =\n\t%d\n", t0);
	fprintf(fp, "#spread =\n\t%f\n", spread);
	fprintf(fp, "#T =\n\t%d\n", T);
	fprintf(fp, "#NSTEPS =\n\t%d\n", NSTEPS);
	fprintf(fp, "#Ex_low_m1 =\n\t%f\n", Ex_low_m1);
	fprintf(fp, "#Ex_low_m2 =\n\t%f\n", Ex_low_m2);
	fprintf(fp, "#Ex_high_m1 =\n\t%f\n", Ex_high_m1);
	fprintf(fp, "#Ex_high_m2 =\n\t%f\n", Ex_high_m2);
	fprintf(fp, "\n");
	fprintf(fp, "########################################################################\n");
	fprintf(fp, "#Ex\n");
	for (k = 0; k < KE + 1; k++)
		fprintf(fp, "%f\t", Ex[k]);
	fprintf(fp, "\n");
	fprintf(fp, "#Hy\n");
	for (k = 0; k < KE; k++)
		fprintf(fp, "%f\t", Hy[k]);
	fprintf(fp, "\n");
	fclose(fp);*/
	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, "kstart = %d;\n", kstart);
	fprintf(fp, "epsilon = %f;\n", epsilon);
	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);
	
	printf("KE=%d\tkc=%d\tt0=%d\tspread=%f\tT=%d\tNSTEPS=%d\n", KE, kc, t0, spread, T, NSTEPS);
	free(Ex);
	free(Hy);
}