fdtd_init.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_init.exe fileName"							*/

#include <stdlib.h>
#include <stdio.h>
#define epsz 8.854187817e-12

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, epsilon, sigma, eaf;
	int n, k, kc, T, NSTEPS, KE, kstart;
	FILE *fp;
	
	T = 0;				/* Initial step							*/
	NSTEPS = 10;			/* run step this time					*/
	KE = 200;			/* KE is the number of cells to be used */
	ddx = 0.01;			/* Set the cell size to 1cm				*/
	dt = ddx / (2*3e8);	/* Calculate the time step				*/
	freq_in = 700*1e6;	/* frequency (Hz) 						*/
	Ex_low_m1 = 0;		/* set Absorbing Boundary Conditions	*/
	Ex_low_m2 = 0;
	Ex_high_m1 = 0;
	Ex_high_m2 = 0;
	
	/* Initialize */
	kstart = KE / 2;
	epsilon = 4;
	sigma = 0.04;
	ca = (float *)malloc(KE * sizeof(float));
	cb = (float *)malloc(KE * sizeof(float));
	for (k = 0; k < KE; k++) {
		ca[k] = 1.0;
		cb[k] = 0.5;
	}
	eaf = dt * sigma / (2 * epsz * epsilon);
	for (k = kstart; k < KE; k++) {
		ca[k] = (1.0 - eaf) / (1 + eaf);
		cb[k] = 0.5/(epsilon * (1 + eaf));
	}
	
	Ex = (float *)malloc((KE + 1) * sizeof(float));
	for (k = 0; k < KE + 1; k++)
		Ex[k] = 0.0;	/* Initialize Ex field	*/
	
	Hy = (float *)malloc(KE * sizeof(float));
	for (k = 0; k < KE; k++)
		Hy[k] = 0.0; 	/* Initialize Hy field	*/
	
	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);
	
	return 0;
}