fdtd.cpp

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

#include "fdtd.h"
#define epsz 8.854187817e-12
#define pi 3.14159265358979323846
#define LightSpeed 299792458

Materials::Materials()
{
	KE = 200;
	ddx = 0.01;
	dt = ddx / LightSpeed;
	epsilon = 4;
	sigma = 0.04;
}

double *Materials::setGa()
{
	double *ga;int k, kstart = KE/2;
	ga = new double[KE - 1];
	for (k = 0; k < KE - 1; k++)
		ga[k] = 1.0;
	for (k = kstart; k < KE - 1; k++)
		ga[k] = 1.0 / (epsilon + sigma * dt / epsz);
	return ga;
}

double *Materials::setGb()
{
	double *gb;int k, kstart = KE/2;
	gb = new double[KE - 1];
	for (k = 0; k < KE - 1; k++)
		gb[k] = 0.0;
	for (k = kstart; k < KE - 1; k++)
		gb[k] = sigma * dt / epsz;
	return gb;
}

int Materials::setKE()
{
	return KE;
}

void FDTD::initialize()
{
	int k;
	/* Initialize */
	material = new Materials;
	ga = material->setGa();
	gb = material->setGb();
	
	T = 0;
	NSTEPS = 5;
	KE = material->setKE();
	
	Ex = new double[KE + 1];
	for (k = 0; k < KE + 1; k++)
		Ex[k] = 0.0;
	Hy = new double[KE];
	for (k = 0; k < KE; k++)
		Hy[k] = 0.0;
		
	//ga = new double[KE - 1];
	//gb = new double[KE - 1];
	Ix = new double[KE - 1];
	Dx = new double[KE - 1];
	for (k = 0; k < KE - 1; k++) {
		//ga[k] = 1.0;
		//gb[k] = 0.0;
		Ix[k] = 0.0;
		Dx[k] = 0.0;
	}
	
	/* set Absorbing Boundary Conditions	*/
	boundaryType = 1;	/* Boundary type equals one				*/
	boundary = new Boundary;
	
	/* Set Source */
	sourceType = 2;		/* Source type ? =1:Sine,=2:Gaussian	*/
	source = new Source(sourceType);
}

FDTD::~FDTD()
{
	delete []ga;
	delete []gb;
	delete []Ix;
	delete []Dx;
	delete []Ex;
	delete []Hy;
	delete material;
	delete boundary;
	delete source;
}

void FDTD::readFile(QTextStream &stream)
{
	QString line;
	QStringList list;
	int n = 0;
	do {
		line = stream.readLine();
		if (line.isEmpty() || line[0] == '#')
			continue;
		n = n + 1;
		switch (n) {
		case 1:
			T = line.toInt();
			break;
		case 2:
			NSTEPS = line.toInt();
			break;
		case 3:
			KE = line.toInt();
			Ex = new double[KE + 1];
			Hy = new double[KE];
			ga = new double[KE - 1];
			gb = new double[KE - 1];
			Ix = new double[KE - 1];
			Dx = new double[KE - 1];
			break;
		case 4:
			boundaryType = line.toInt();
			boundary = new Boundary;
			break;
		case 5:
			sourceType = line.toInt();
			source = new Source(sourceType);
			break;
		case 6:
			list = line.split('\t');
			for (int k = 0; k < KE + 1; k++)
				Ex[k] = list[k].toDouble();
			break;
		case 7:
			list = line.split('\t');
			for (int k = 0; k < KE; k++)
				Hy[k] = list[k].toDouble();
			break;
		case 8:
			list = line.split('\t');
			for (int k = 0; k < KE - 1; k++)
				ga[k] = list[k].toDouble();
			break;
		case 9:
			list = line.split('\t');
			for (int k = 0; k < KE - 1; k++)
				gb[k] = list[k].toDouble();
			break;
		case 10:
			list = line.split('\t');
			for (int k = 0; k < KE - 1; k++)
				Ix[k] = list[k].toDouble();
			break;
		case 11:
			list = line.split('\t');
			for (int k = 0; k < KE - 1; k++)
				Dx[k] = list[k].toDouble();
			break;
		}
	} while (n != 11);
	boundary->readFile(stream);
	source->readFile(stream);
}

void FDTD::writeFile(QTextStream &stream)
{
	int k;
	
	stream << "########################################################################" << endl;
	stream << "####################Electromagnetic simulation##########################" << endl;
	stream << endl;
	stream << "########################################################################" << endl;
	stream << "#The parameters of FDTD" << endl;
	stream << "#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::" << endl;
	stream << "#The number of total time steps up to now" << endl;
	stream << "\t" << T << endl;
	stream << "#The number of steps for running this time" << endl;
	stream << "\t" << NSTEPS << endl;
	stream << "#The number of cells to be used" << endl;
	stream << "\t" << KE << endl;
	stream << "#Boundary Type? =1:simple" << endl;
	stream << "\t" << boundaryType << endl;
	stream << "#Source Type? =1:Sine source; =2:Gaussian source" << endl;
	stream << "\t" << sourceType << endl;
	stream << "#The Ex field" << endl;
	for (k = 0; k < KE + 1; k++)
		stream << Ex[k] << "\t";
	stream << endl << "#The Hy field" << endl;
	for (k = 0; k < KE; k++)
		stream << Hy[k] << "\t";
	stream << endl << "#The ga value" << endl;
	for (k = 0; k < KE - 1; k++)
		stream << ga[k] << "\t";
	stream << endl << "#The gb value" << endl;
	for (k = 0; k < KE - 1; k++)
		stream << gb[k] << "\t";
	stream << endl << "#The Ix field" << endl;
	for (k = 0; k < KE - 1; k++)
		stream << Ix[k] << "\t";
	stream << endl << "#The Dx field" << endl;
	for (k = 0; k < KE - 1; k++)
		stream << Dx[k] << "\t";
	stream << endl;
	stream << "########################################################################\n" << endl;
	boundary->writeFile(stream);
	source->writeFile(stream);
}

void FDTD::simulate()
{
	int n,k;
	
	/* Main part of the FDTD simulation */
	for (n = 1; n <= NSTEPS; n++) {
		T = T + 1;	/* T keeps track of the total number of times the main loop is executed */
		
		qDebug() << "T=" << T << "\tNSTEPS=" << NSTEPS << endl;
		
		/* Calculate the Dx field */
		for (k = 0; k < KE -1; k++)
			Dx[k] = Dx[k] + 0.5 * (Hy[k] - Hy[k+1]);
		
		/* Put a Source at the low end */
		source->setSource(Dx, T);
		
		/* Calculate Ex from Dx */
		for(k = 1; k < KE; k++) {
			Ex[k] = ga[k-1] * (Dx[k-1] - Ix[k-1]);
			Ix[k-1] = Ix[k-1] + gb[k-1]*Ex[k];
		}
		
		/* Boundary conditions */
		boundary->setBoundary(Ex, KE);
		
		/*Calculate the Hy field */
		for (k = 0; k < KE; k++)
			Hy[k] = Hy[k] + 0.5 * (Ex[k] - Ex[k+1]);
	}
}

Boundary::Boundary()
{
	Ex_low_m1 = 0;
	Ex_low_m2 = 0;
	Ex_high_m1 = 0;
	Ex_high_m2 = 0;
}

void Boundary::readFile(QTextStream &stream)
{
	QString line;
	int n = 0;
	do {
		line = stream.readLine();
		if (line.isEmpty() || line[0] == '#')
			continue;
		n = n + 1;
		switch (n) {
		case 1:
			Ex_low_m1 = line.toDouble();
			break;
		case 2:
			Ex_low_m2 = line.toDouble();
			break;
		case 3:
			Ex_high_m1 = line.toDouble();
			break;
		case 4:
			Ex_high_m2 = line.toDouble();
			break;
		}
	} while (n != 4);
}

void Boundary::writeFile(QTextStream &stream)
{
	stream << "########################################################################" << endl;
	stream << "#The parameters of Boundary" << endl;
	stream << "#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::" << endl;
	stream << "#Ex_low_m1:" << endl;
	stream << "\t" << Ex_low_m1 << endl;
	stream << "#Ex_low_m2:" << endl;
	stream << "\t" << Ex_low_m2 << endl;
	stream << "#Ex_high_m1:" << endl;
	stream << "\t" << Ex_high_m1 << endl;
	stream << "#Ex_high_m2:" << endl;
	stream << "\t" << Ex_high_m2 << endl;
	stream << "########################################################################\n" << endl;
}

void Boundary::setBoundary(double *Ex, int KE)
{
	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];
}

Source::Source(int sourceType)
{
	type = sourceType;
	switch (type) {
	case 1:
		SineInit();
		break;
	case 2:
		GaussianInit();
		break;
	default :
		SineInit();
	}
}

void Source::readFile(QTextStream &stream)
{
	switch (type) {
	case 1:
		SineRead(stream);
		break;
	case 2:
		GaussianRead(stream);
		break;
	default :
		SineRead(stream);
	}
}

void Source::writeFile(QTextStream &stream)
{
	switch (type) {
	case 1:
		SineWrite(stream);
		break;
	case 2:
		GaussianWrite(stream);
		break;
	default :
		SineWrite(stream);
	}
}

void Source::setSource(double *Dx, int T)
{
	int position = 5;
	switch (type) {
	case 1:
		Dx[position] = Dx[position] + SinePulse(T);
		break;
	case 2:
		Dx[position] = Dx[position] + GaussianPulse(T);
		break;
	default :
		Dx[position] = Dx[position] + SinePulse(T);
	}
}

void Source::SineInit()
{
	ddx = 0.01;;
	dt = ddx / LightSpeed;
	freq_in = 700*1e6;
}

void Source::SineRead(QTextStream &stream)
{
	QString line;
	int n = 0;
	do {
		line = stream.readLine();
		if (line.isEmpty() || line[0] == '#')
			continue;
		n = n + 1;
		switch (n) {
		case 1:
			ddx = line.toDouble();
			break;
		case 2:
			dt = line.toDouble();
			break;
		case 3:
			freq_in = line.toDouble();
			break;
		}
	} while (n != 3);
}

void Source::SineWrite(QTextStream &stream)
{
	stream << "########################################################################" << endl;
	stream << "#The parameters of Sine Source" << endl;
	stream << "#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::" << endl;
	stream << "#the cell size:" << endl;
	stream << "\t" << ddx << endl;
	stream << "#the time step size:" << endl;
	stream << "\t" << dt << endl;
	stream << "#frequency (Hz):" << endl;
	stream << "\t" << freq_in << endl;
	stream << "########################################################################\n" << endl;
}

double Source::SinePulse(int T)
{
	return sin(2 * pi * freq_in * dt * T);
}

void Source::GaussianInit()
{
	t0 = 40;
	spread = 12;
}

void Source::GaussianRead(QTextStream &stream)
{
	QString line;
	int n = 0;
	do {
		line = stream.readLine();
		if (line.isEmpty() || line[0] == '#')
			continue;
		n = n + 1;
		switch (n) {
		case 1:
			t0 = line.toDouble();
			break;
		case 2:
			spread = line.toDouble();
			break;
		}
	} while (n != 2);
}

void Source::GaussianWrite(QTextStream &stream)
{
	stream << "########################################################################" << endl;
	stream << "#The parameters of Gaussian Source" << endl;
	stream << "#:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::" << endl;
	stream << "#Center of the incident pulse	:" << endl;
	stream << "\t" << t0 << endl;
	stream << "#Width of the incident pulse	:" << endl;
	stream << "\t" << spread << endl;
	stream << "########################################################################\n" << endl;
}

double Source::GaussianPulse(int T)
{
	return exp(-0.5 * ( pow((t0 - T) / spread, 2.0)));
}