fdtd.txt

fdtd的一些一维到三维的程序,通常为pml边界,用txt文本编辑,无需解压密码

			for (x=1;x<=PCsize-2;x++)
			{
				for (y=2;y<=PCsize-3;y++)	// Set values within the triDiag matrix
				{
					triDiag[1][y-1] = coef1zsub1[x][y][z];
					triDiag[0][y-1] = -coef2zsub1[x][y][z];
					triDiag[2][y-1] = -coef2zsub1[x][y][z];
					triDiag[3][y-1] = coef3zsub1[x][y][z]*Ezsub1[x][y][z]
						+ coef4zsub1[x][y][z]*(Hy[x+1][y][z]-Hy[x][y][z]-Hx[x][y+1][z]+Hx[x][y][z])
						- coef2zsub1[x][y][z]*(Eysub1[x][y+1][z]-Eysub1[x][y+1][z-1]-Eysub1[x][y][z]+Eysub1[x][y][z-1]);
				}
				triDiag[1][0] = coef1zsub1[x][1][z];
				triDiag[2][0] = -coef2zsub1[x][1][z];
				triDiag[3][0] = coef3zsub1[x][1][z]*Ezsub1[x][1][z]
					+ coef4zsub1[x][1][z]*(Hy[x+1][1][z]-Hy[x][1][z]-Hx[x][2][z]+Hx[x][1][z])
					- coef2zsub1[x][1][z]*(Eysub1[x][2][z]-Eysub1[x][2][z-1]-Eysub1[x][1][z]+Eysub1[x][1][z-1]);
				triDiag[1][PCsize-3] = coef1zsub1[x][PCsize-2][z];
				triDiag[0][PCsize-3] = -coef2zsub1[x][PCsize-2][z];
				triDiag[3][PCsize-3] = coef3zsub1[x][PCsize-2][z]*Ezsub1[x][PCsize-2][z]
					+ coef4zsub1[x][PCsize-2][z]*(Hy[x+1][PCsize-2][z]-Hy[x][PCsize-2][z]-Hx[x][PCsize-1][z]+Hx[x][PCsize-2][z])
					- coef2zsub1[x][PCsize-2][z]*(Eysub1[x][PCsize-1][z]-Eysub1[x][PCsize-1][z-1]-Eysub1[x][PCsize-2][z]+Eysub1[x][PCsize-2][z-1]);
				solve(PCsize-2,triDiag);			// Solve Tridiagonal Matrix
				for (y=1;y<=PCsize-2;y++)
				{
					Ezsub1[x][y][z]=triDiag[3][y-1];	// Set values in alternate array
				}
			}
			// H field FDTD equations are identical to explicit scheme
		for (z=1;z<=PCsize-2;z++)
			for (y=1;y<=PCsize-2;y++)
				for (x=1;x<=PCsize-2;x++)
				{
					Hx[x][y][z] +=  Hcoef*(Eysub1[x][y][z]-Eysub1[x][y][z-1]-Ezsub1[x][y][z]+Ezsub1[x][y-1][z]);
					Hy[x][y][z] +=  Hcoef*(Ezsub1[x][y][z]-Ezsub1[x-1][y][z]-Exsub1[x][y][z]+Exsub1[x][y][z-1]);
					Hz[x][y][z] +=  Hcoef*(Exsub1[x][y][z]-Exsub1[x][y-1][z]-Eysub1[x][y][z]+Eysub1[x-1][y][z]);
				}
		//////////////////////////End Of Sub-iteration 2///////////////////////////////////////////////////////
	}

			//	for(t=0;t<PCtimeSteps;t++)
//	{
//		for (z=10;z<PCzWidth-10;z++)	//H portion of FDTD - DFT computed within//
//			for (yIndex=0;yIndex<PCloBoundsFile;yIndex++)//index value of the bounds matrix//
//			{
//				y = bounds[yIndex].y; //actual y value is stored in bounds//
//				for (x=bounds[yIndex].xmin+1;x<bounds[yIndex].xmax;x++)
//				{
//					Hx[PMxyz]=Hx[PMxyz]-HCoef2*((Ez[PMxyz+PCyOffset]-Ez[PMxyz])/PCdeltaY-(Ey[PMxyz+PCzOffset]-Ey[PMxyz])/PCdeltaZ);
//					HxFDR[PMxyz]=HxFDR[PMxyz]+Hx[PMxyz]*realDFT;//DFT//
//					HxFDI[PMxyz]=HxFDI[PMxyz]+Hx[PMxyz]*imagDFT;//DFT//
//					Hy[PMxyz]=Hy[PMxyz]-HCoef2*((Ex[PMxyz+PCzOffset]-Ex[PMxyz])/PCdeltaZ-(Ez[PMxyz+PCxOffset]-Ez[PMxyz])/PCdeltaX);
//					HyFDR[PMxyz]=HyFDR[PMxyz]+Hy[PMxyz]*realDFT;//DFT//
//					HyFDI[PMxyz]=HyFDI[PMxyz]+Hy[PMxyz]*imagDFT;//DFT//
//					Hz[PMxyz]=Hz[PMxyz]-HCoef2*((Ey[PMxyz+PCxOffset]-Ey[PMxyz])/PCdeltaX-(Ex[PMxyz+PCyOffset]-Ex[PMxyz])/PCdeltaY);
//				}
//			}
//	}
}
int initialise()	//Initialise Data//
{
	printf("Initialising\n");
	int x,y,z;
	for (z=0;z<=PCsize-1;z++)
		for (y=0;y<=PCsize-1;y++)
			for (x=0;x<=PCsize-1;x++)
			{
				epsilon[x][y][z]=PCepsilonNot;
  // Setinitial Values
				sigma[x][y][z]=1e-30;;
	  // Use 1e-30 as zero;
				Exsub1[x][y][z]=1e-30;
				Exsub2[x][y][z]=1e-30;
				Eysub1[x][y][z]=1e-30;
				Eysub2[x][y][z]=1e-30;
				Ezsub1[x][y][z]=1e-30;
				Ezsub2[x][y][z]=1e-30;
				Hx[x][y][z]=1e-30;
				Hy[x][y][z]=1e-30;
				Hz[x][y][z]=1e-30;
				Ix[x][y][z]=1e-30;
			}
	for (z=1;z<=PCsize-2;z++)
		for (y=1;y<=PCsize-2;y++)
			for (x=1;x<=PCsize-2;x++)
			{
				coef1xsub1[x][y][z] = 1+((PCdeltaT*PCdeltaT/(PCmuNot*(epsilon[x][y][z]+epsilon[x-1][y][z])*PCdeltaS*PCdeltaS))/(1+((sigma[x][y][z]+sigma[x-1][y][z])*PCdeltaT)/(2*(epsilon[x][y][z]+epsilon[x-1][y][z]))));
				coef2xsub1[x][y][z] = (coef1xsub1[x][y][z] - 1) * 0.5;
				coef3xsub1[x][y][z] = (1-((sigma[x][y][z]+sigma[x-1][y][z])*PCdeltaT)/(2*(epsilon[x][y][z]+epsilon[x-1][y][z])))/(1+((sigma[x][y][z]+sigma[x-1][y][z])*PCdeltaT)/(2*(epsilon[x][y][z]+epsilon[x-1][y][z])));
				coef4xsub1[z][y][z] = ((PCdeltaT)/((epsilon[x][y][z]+epsilon[x-1][y][z])*PCdeltaS))/(1+((sigma[x][y][z]+sigma[x-1][y][z])*PCdeltaT)/(2*(epsilon[x][y][z]+epsilon[x-1][y][z])));
				coef5xsub1[x][y][z] = (PCdeltaT/(epsilon[x][y][z]+epsilon[x-1][y][z]))/(1+((sigma[x][y][z]+sigma[x-1][y][z])*PCdeltaT)/(2*(epsilon[x][y][z]+epsilon[x-1][y][z])));
					// Ey Coeffecients
				coef1ysub1[x][y][z] = 1+((PCdeltaT*PCdeltaT/(PCmuNot*(epsilon[x][y][z]+epsilon[x][y-1][z])*PCdeltaS*PCdeltaS))/(1+((sigma[x][y][z]+sigma[x][y-1][z])*PCdeltaT)/(2*(epsilon[x][y][z]+epsilon[x][y-1][z]))));
				coef2ysub1[x][y][z] = (coef1xsub1[x][y][z] - 1) * 0.5;
				coef3ysub1[x][y][z] = (1-((sigma[x][y][z]+sigma[x][y-1][z])*PCdeltaT)/(2*(epsilon[x][y][z]+epsilon[x][y-1][z])))/(1+((sigma[x][y][z]+sigma[x][y-1][z])*PCdeltaT)/(2*(epsilon[x][y][z]+epsilon[x][y-1][z])));
				coef4ysub1[z][y][z] = ((PCdeltaT)/((epsilon[x][y][z]+epsilon[x][y-1][z])*PCdeltaS))/(1+((sigma[x][y][z]+sigma[x][y-1][z])*PCdeltaT)/(2*(epsilon[x][y][z]+epsilon[x][y-1][z])));
					// Ez Coefficients
				coef1zsub1[x][y][z] = 1+((PCdeltaT*PCdeltaT/(PCmuNot*(epsilon[x][y][z]+epsilon[x][y][z-1])*PCdeltaS*PCdeltaS))/(1+((sigma[x][y][z]+sigma[x][y][z-1])*PCdeltaT)/(2*(epsilon[x][y][z]+epsilon[x][y][z-1]))));
				coef2zsub1[x][y][z] = (coef1xsub1[x][y][z] - 1) * 0.5;
				coef3zsub1[x][y][z] = (1-((sigma[x][y][z]+sigma[x][y][z-1])*PCdeltaT)/(2*(epsilon[x][y][z]+epsilon[x][y][z-1])))/(1+((sigma[x][y][z]+sigma[x][y][z-1])*PCdeltaT)/(2*(epsilon[x][y][z]+epsilon[x][y][z-1])));
				coef4zsub1[z][y][z] = ((PCdeltaT)/((epsilon[x][y][z]+epsilon[x][y][z-1])*PCdeltaS))/(1+((sigma[x][y][z]+sigma[x][y][z-1])*PCdeltaT)/(2*(epsilon[x][y][z]+epsilon[x][y][z-1])));
			}
	Hcoef= PCdeltaT / (2*PCmuNot*PCdeltaS);
}

/*int saveDataX(int x,char *outputFileName)//Saves data in x plane//
{
	FILE *output;

	int y,z;
	float xAbs, yAbs, tot;

	output = fopen(outputFileName,"w");
	for(y=0;y<=PCyWidth-1;y++)
	{
		for (z=0;z<=PCzWidth-1;z++)
		{
			xAbs = sqrt(HxFDR[PMxyz]*HxFDR[PMxyz]+HxFDI[PMxyz]*HxFDI[PMxyz]);//absolute Hx//
			yAbs = sqrt(HyFDR[PMxyz]*HyFDR[PMxyz]+HyFDI[PMxyz]*HyFDI[PMxyz]);//absolute Hy//
			tot = sqrt(xAbs*xAbs+yAbs*yAbs);//B1 total//
			fprintf(output,"%g ",tot);
		}
		fprintf(output,"\n");
	}
	 fclose(output);
	 return 0;
}
int saveDataY(int y,char *outputFileName)//Saves data in y plane//
{
	FILE *output;

	int x,z;
	float xAbs, yAbs, tot;

	output = fopen(outputFileName,"w");
	for(x=0;x<=PCxWidth-1;x++)
	{
		for (z=0;z<=PCzWidth-1;z++)
		{
			xAbs = sqrt(HxFDR[PMxyz]*HxFDR[PMxyz]+HxFDI[PMxyz]*HxFDI[PMxyz]);//absolute Hx//
			yAbs = sqrt(HyFDR[PMxyz]*HyFDR[PMxyz]+HyFDI[PMxyz]*HyFDI[PMxyz]);//absolute Hy//
			tot = sqrt(xAbs*xAbs+yAbs*yAbs);//B1 total//
			fprintf(output,"%g ",tot);
		}
		fprintf(output,"\n");
	}
	 fclose(output);
	 return 0;
}
int saveDataZ(int z,char *outputFileName)//Saves data in z plane//
{
	FILE *output;

	int x,y;
	float xAbs, yAbs, tot;

	output = fopen(outputFileName,"w");
	for(x=0;x<=PCxWidth-1;x++)
	{
		for (y=0;y<=PCyWidth-1;y++)
		{
			xAbs = sqrt(HxFDR[PMxyz]*HxFDR[PMxyz]+HxFDI[PMxyz]*HxFDI[PMxyz]);//absolute Hx//
			yAbs = sqrt(HyFDR[PMxyz]*HyFDR[PMxyz]+HyFDI[PMxyz]*HyFDI[PMxyz]);//absolute Hy//
			tot = sqrt(xAbs*xAbs+yAbs*yAbs);//B1 total//
			fprintf(output,"%g ",tot);
		}
		fprintf(output,"\n");
	}
	 fclose(output);
	 return 0;
}
*/

//int loadData(char *inputFileName)//Loads material data from a given file//
//{
//	FILE *inputFile;	//input file object//
//
//	long size;			//size of the text file//
//	int x,y,z;			//cartesian values//
//	float mEpsilon, mSigma;	//material data//
//
//	inputFile = fopen(inputFileName,"r");	//opens input file//
//	if (inputFile == NULL) return 1;		//error condition//
//	fseek(inputFile,0,SEEK_END);			//seek to the end//
//	size = ftell(inputFile);				//get file size//
//	fseek(inputFile,0,SEEK_SET);			//go back to the begining//
//	do{
//		x = getnexttoken(inputFile)+PCxShift;	//get cartesian coordinates//
//		y = getnexttoken(inputFile)+PCyShift;	//shift down to account for numbering//
//		z = getnexttoken(inputFile)+PCzShift;
//		mEpsilon = getnexttokenf(inputFile);	//get material data//
//		mSigma = getnexttokenf(inputFile);
//		epsilon[PMxyz]=mEpsilon*PCepsilonNot;
//		sigma[PMxyz]=mSigma;
//		unitArray[PMxyz]=1;
//	}while(size-ftell(inputFile)>10);	//loop until near eof//
//	fclose(inputFile);
//	return 0;
//}

//int crop()
//{
//	int x,y,z,yIndex;
//	for (z=10;z<PCzWidth-10;z++)	//H portion of FDTD - DFT computed within//
//		for (yIndex=0;yIndex<PCloBoundsFile;yIndex++)//index value of the bounds matrix//
//		{
//			y = bounds[yIndex].y; //actual y value is stored in bounds//
//			for (x=bounds[yIndex].xmin+1;x<bounds[yIndex].xmax;x++)
//			{
//				HxFDR[PMxyz]=HxFDR[PMxyz]*((double)unitArray[PMxyz]);
//				HyFDR[PMxyz]=HxFDR[PMxyz]*((double)unitArray[PMxyz]);
//				HxFDI[PMxyz]=HxFDR[PMxyz]*((double)unitArray[PMxyz]);
//				HyFDI[PMxyz]=HxFDR[PMxyz]*((double)unitArray[PMxyz]);
//			}
//		}
//}

inline int solve (int n, float array[][n])  //Gaussian Elimination (No zero checking)
{
	int j;  //array index
	float c;  //array manipulation operation value
	for (j=0;j<n-1;j++)  //main loop, step thru j's
	{
		c = array[0][j+1]/array[1][j];
		array[1][j+1] -= array[2][j] * c;
		array[3][j+1] -= array[3][j] * c;
	}
	for (j=n-1;j>=0;j--)  //main loop, step thru j's
	{
		c = array[2][j]/array[1][j+1];//
		array[3][j] -= array[3][j+1] * c;
		array[3][j] *= 1/array[1][j];
		array[1][j] *= 1/array[1][j];
	}
	array[3][n-1] *= 1/array[1][n-1];
}