📄 fdtd.txt
字号:
#include <stdio.h>
//#include "tokenizer.h"
#include <math.h>
#define PCepsilonNot 8.85e-12
#define PCpi 3.14159265359
#define PCmuNot 12.57e-7
#define PCdeltaT 2.5e-9
#define PCdeltaX 0.02
#define PCdeltaY 0.02
#define PCdeltaZ 0.02
#define PCdeltaS 0.02
#define PCtimeSteps 20
#define PCfreq 400e6
#define PCsize 100
//int saveDataX(int x,char *outputFileName);
//int saveDataY(int y,char *outputFileName);
//int saveDataZ(int z,char *outputFileName);
int initialise();
//int setCoefficients();
inline int solve (int n, float array[][n]);
int fdtd();
float Exsub1[PCsize][PCsize][PCsize];
float Exsub2[PCsize][PCsize][PCsize];
float Eysub1[PCsize][PCsize][PCsize];
float Eysub2[PCsize][PCsize][PCsize];
float Ezsub1[PCsize][PCsize][PCsize];
float Ezsub2[PCsize][PCsize][PCsize];
float Ix[PCsize][PCsize][PCsize];
float Hx[PCsize][PCsize][PCsize];
float Hy[PCsize][PCsize][PCsize];
float Hz[PCsize][PCsize][PCsize];
float triDiag[4][PCsize-2];
float epsilon[PCsize][PCsize][PCsize];
float sigma[PCsize][PCsize][PCsize];
float coef1xsub1[PCsize][PCsize][PCsize];
float coef2xsub1[PCsize][PCsize][PCsize];
float coef3xsub1[PCsize][PCsize][PCsize];
float coef4xsub1[PCsize][PCsize][PCsize];
float coef5xsub1[PCsize][PCsize][PCsize];
float coef1ysub1[PCsize][PCsize][PCsize];
float coef2ysub1[PCsize][PCsize][PCsize];
float coef3ysub1[PCsize][PCsize][PCsize];
float coef4ysub1[PCsize][PCsize][PCsize];
float coef1zsub1[PCsize][PCsize][PCsize];
float coef2zsub1[PCsize][PCsize][PCsize];
float coef3zsub1[PCsize][PCsize][PCsize];
float coef4zsub1[PCsize][PCsize][PCsize];
float Hcoef;
int main()
{
printf("Begin\n");
initialise();
// setCoefficients();
fdtd();
// crop();
printf("Success\n");
}
int fdtd()
{
printf("FDTD\n");
int x,y,z,t;
for (t=0;t<=PCtimeSteps-1;t++)
{
///////////////////////////////////////Sub-iteration 1//////////////////////////////////////////////////////////
for (z=1;z<=PCsize-2;z++)
for (x=1;x<=PCsize-2;x++)
{
for (y=2;y<=PCsize-3;y++) // Set values within the triDiag matrix
{
triDiag[1][y-1] = coef1xsub1[x][y][z];
triDiag[0][y-1] = -coef2xsub1[x][y][z];
triDiag[2][y-1] = -coef2xsub1[x][y][z];
triDiag[3][y-1] = coef3xsub1[x][y][z]*Exsub1[x][y][z]
+ coef4xsub1[x][y][z]*(Hz[x][y+1][z]-Hz[x][y][z]-Hy[x][y][z+1]+Hy[x][y][z])
- coef2xsub1[x][y][z]*(Eysub1[x][y+1][z]-Eysub1[x-1][y+1][z]-Eysub1[x][y][z]+Eysub1[x-1][y][z])
;//Ix -> - coef5xsub1[x][y][z]*Ix[x][y][z];
}
triDiag[1][0] = coef1xsub1[x][1][z]; // Set values on the ends of the tridiagonal
triDiag[2][0] = -coef2xsub1[x][1][z];
triDiag[3][0] = coef3xsub1[x][1][z]*Exsub1[x][1][z]
+ coef4xsub1[x][1][z]*(Hz[x][2][z]-Hz[x][1][z]-Hy[x][1][z+1]+Hy[x][1][z])
- coef2xsub1[x][1][z]*(Eysub1[x][2][z]-Eysub1[x-1][2][z]-Eysub1[x][1][z]+Eysub1[x-1][1][z])
;//Ix -> - coef5xsub1[x][1][z]*Ix[x][1][z];
triDiag[1][PCsize-3] = coef1xsub1[x][PCsize-2][z];
triDiag[0][PCsize-3] = -coef2xsub1[x][PCsize-2][z];
triDiag[3][PCsize-3] = coef3xsub1[x][PCsize-2][z]*Exsub1[x][PCsize-2][z]
+ coef4xsub1[x][PCsize-2][z]*(Hz[x][PCsize-1][z]-Hz[x][PCsize-2][z]-Hy[x][PCsize-2][z+1]+Hy[x][PCsize-2][z])
- coef2xsub1[x][PCsize-2][z]*(Eysub1[x][PCsize-1][z]-Eysub1[x-1][PCsize-1][z]-Eysub1[x][PCsize-2][z]+Eysub1[x-1][PCsize-2][z])
;//Ix -> - coef5xsub1[x][PCsize-2][z]*Ix[x][PCsize-2][z];
solve(PCsize-2,triDiag); // Solve Tridiagonal Matrix
for (y=1;y<=PCsize-2;y++)
{
Exsub2[x][y][z]=triDiag[3][y-1]; // Set values in alternate array
}
}
for (y=1;y<=PCsize-2;y++)
for (x=1;x<=PCsize-2;x++)
{
for (z=2;z<=PCsize-3;z++) // Set values within the triDiag matrix
{
triDiag[1][z-1] = coef1ysub1[x][y][z];
triDiag[0][z-1] = -coef2ysub1[x][y][z];
triDiag[2][z-1] = -coef2ysub1[x][y][z];
triDiag[3][z-1] = coef3ysub1[x][y][z]*Eysub1[x][y][z]
+ coef4ysub1[x][y][z]*(Hx[x][y][z+1]-Hx[x][y][z]-Hz[x+1][y][z]+Hz[x][y][z])
- coef2ysub1[x][y][z]*(Ezsub1[x][y][z+1]-Ezsub1[x][y-1][z+1]-Ezsub1[x][y][z]+Ezsub1[x][y-1][z]);
}
triDiag[1][0] = coef1ysub1[x][y][1];
triDiag[2][0] = -coef2ysub1[x][y][1];
triDiag[3][0] = coef3ysub1[x][y][1]*Eysub1[x][y][1]
+ coef4ysub1[x][y][1]*(Hx[x][y][2]-Hx[x][y][z]-Hz[2][y][z]+Hz[x][y][1])
- coef2ysub1[x][y][1]*(Ezsub1[x][y][2]-Ezsub1[x][y-1][2]-Ezsub1[x][y][1]+Ezsub1[x][y-1][1]);
triDiag[1][PCsize-3] = coef1ysub1[x][y][PCsize-2];
triDiag[0][PCsize-3] = -coef2ysub1[x][y][PCsize-2];
triDiag[3][PCsize-3] = coef3ysub1[x][y][PCsize-2]*Eysub1[x][y][PCsize-2]
+ coef4ysub1[x][y][PCsize-2]*(Hx[x][y][PCsize-1]-Hx[x][y][PCsize-2]-Hz[x+1][y][PCsize-2]+Hz[x][y][PCsize-2])
- coef2ysub1[x][y][PCsize-2]*(Ezsub1[x][y][PCsize-1]-Ezsub1[x][y-1][PCsize-1]-Ezsub1[x][y][PCsize-2]+Ezsub1[x][y-1][PCsize-2]);
solve(PCsize-2,triDiag); // Solve Tridiagonal Matrix
for (z=1;z<=PCsize-2;z++)
{
Eysub2[x][y][z]=triDiag[3][z-1]; // Set values in alternate array
}
}
for (z=1;z<=PCsize-2;z++)
for (y=1;y<=PCsize-2;y++)
{
for (x=2;x<=PCsize-3;x++) // Set values within the triDiag matrix
{
triDiag[1][x-1] = coef1zsub1[x][y][z];
triDiag[0][x-1] = -coef2zsub1[x][y][z];
triDiag[2][x-1] = -coef2zsub1[x][y][z];
triDiag[3][x-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]*(Exsub1[x+1][y][z]-Exsub1[x+1][y][z-1]-Exsub1[x][y][z]+Exsub1[x][y][z-1]);
}
triDiag[1][0] = coef1zsub1[1][y][z];
triDiag[2][0] = -coef2zsub1[1][y][z];
triDiag[3][0] = coef3zsub1[1][y][z]*Ezsub1[1][y][z]
+ coef4zsub1[1][y][z]*(Hy[2][y][z]-Hy[1][y][z]-Hx[1][y+1][z]+Hx[1][y][z])
- coef2zsub1[1][y][z]*(Exsub1[2][y][z]-Exsub1[2][y][z-1]-Exsub1[1][y][z]+Exsub1[1][y][z-1]);
triDiag[1][PCsize-3] = coef1zsub1[PCsize-2][y][z];
triDiag[0][PCsize-3] = -coef2zsub1[PCsize-2][y][z];
triDiag[3][PCsize-3] = coef3zsub1[PCsize-2][y][z]*Ezsub1[PCsize-2][y][z]
+ coef4zsub1[PCsize-2][y][z]*(Hy[PCsize-1][y][z]-Hy[PCsize-2][y][z]-Hx[PCsize-2][y+1][z]+Hx[PCsize-2][y][z])
- coef2zsub1[PCsize-2][y][z]*(Exsub1[PCsize-1][y][z]-Exsub1[PCsize-1][y][z-1]-Exsub1[PCsize-2][y][z]+Exsub1[PCsize-2][y][z-1]);
solve(PCsize-2,triDiag); // Solve Tridiagonal Matrix
for (y=1;y<=PCsize-2;y++)
{
Ezsub2[x][y][z]=triDiag[3][y-1]; // Set values in alternate array
}
}
// H field FDTD equations are identical to explicit scheme
Exsub2[50][50][50] = 1; //Forcing Function For Detection Purposes
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*(Eysub2[x][y][z]-Eysub2[x][y][z-1]-Ezsub2[x][y][z]+Ezsub2[x][y-1][z]);
Hy[x][y][z] += Hcoef*(Ezsub2[x][y][z]-Ezsub2[x-1][y][z]-Exsub2[x][y][z]+Exsub2[x][y][z-1]);
Hz[x][y][z] += Hcoef*(Exsub2[x][y][z]-Exsub2[x][y-1][z]-Eysub2[x][y][z]+Eysub2[x-1][y][z]);
}
printf("Hy 50,50,50 - %f \n",Hy[50][50][50]);
printf("Hz 50,50,50 - %f \n",Hz[50][50][50]);
exit(0);
//////////////////////////End Of Sub-iteration 1///////////////////////////////////////////////////////
///////////////////////////////Sub-iteration 2//////////////////////////////////////////////////////////
for (y=1;y<=PCsize-2;y++)
for (x=1;x<=PCsize-2;x++)
{
for (z=2;z<=PCsize-3;z++) // Set values within the triDiag matrix
{
triDiag[1][z-1] = coef1xsub1[x][y][z];
triDiag[0][z-1] = -coef2xsub1[x][y][z];
triDiag[2][z-1] = -coef2xsub1[x][y][z];
triDiag[3][z-1] = coef3xsub1[x][y][z]*Exsub1[x][y][z]
+ coef4xsub1[x][y][z]*(Hz[x][y+1][z]-Hz[x][y][z]-Hy[x][y][z+1]+Hy[x][y][z])
- coef2xsub1[x][y][z]*(Ezsub1[x][y][z+1]-Ezsub1[x-1][y][z+1]-Ezsub1[x][y][z]+Ezsub1[x-1][y][z])
;//Ix -> - coef5xsub1[x][y][z]*Ix[x][y][z];
}
triDiag[1][0] = coef1xsub1[x][y][1];
triDiag[2][0] = -coef2xsub1[x][y][1];
triDiag[3][0] = coef3xsub1[x][y][1]*Exsub1[x][y][1]
+ coef4xsub1[x][y][1]*(Hz[x][y+1][1]-Hz[x][y][1]-Hy[x][y][2]+Hy[x][y][1])
- coef2xsub1[x][y][1]*(Ezsub1[x][y][2]-Ezsub1[x-1][y][2]-Ezsub1[x][y][1]+Ezsub1[x-1][y][1])
;//Ix-> - coef5xsub1[x][y][1]*Ix[x][y][1];
triDiag[1][PCsize-3] = coef1xsub1[x][y][PCsize-2];
triDiag[0][PCsize-3] = -coef2xsub1[x][y][PCsize-2];
triDiag[3][PCsize-3] = coef3xsub1[x][y][PCsize-2]*Exsub1[x][y][PCsize-2]
+ coef4xsub1[x][y][PCsize-2]*(Hz[x][y+1][PCsize-2]-Hz[x][y][PCsize-2]-Hy[x][y][PCsize-1]+Hy[x][y][PCsize-2])
- coef2xsub1[x][y][PCsize-2]*(Ezsub1[x][y][PCsize-1]-Ezsub1[x-1][y][PCsize-1]-Ezsub1[x][y][PCsize-2]+Ezsub1[x-1][y][PCsize-2])
;//Ix -> - coef5xsub1[x][y][PCsize-2]*Ix[x][y][PCsize-2];
solve(PCsize-2,triDiag); // Solve Tridiagonal Matrix
for (z=1;z<=PCsize-2;z++)
{
Exsub1[x][y][z]=triDiag[3][z-1]; // Set values in alternate array
}
}
for (z=1;z<=PCsize-2;z++)
for (y=1;y<=PCsize-2;y++)
{
for (x=2;x<=PCsize-3;x++) // Set values within the triDiag matrix
{
triDiag[1][x-1] = coef1ysub1[x][y][z];
triDiag[0][x-1] = -coef2ysub1[x][y][z];
triDiag[2][x-1] = -coef2ysub1[x][y][z];
triDiag[3][x-1] = coef3ysub1[x][y][z]*Eysub1[x][y][z]
+ coef4ysub1[x][y][z]*(Hx[x][y][z+1]-Hx[x][y][z]-Hz[x+1][y][z]+Hz[x][y][z])
- coef2ysub1[x][y][z]*(Exsub1[x+1][y][z]-Exsub1[x+1][y-1][z]-Exsub1[x][y][z]+Exsub1[x][y-1][z]);
}
triDiag[1][0] = coef1ysub1[1][y][z];
triDiag[2][0] = -coef2ysub1[1][y][z];
triDiag[3][0] = coef3ysub1[1][y][z]*Eysub1[1][y][z]
+ coef4ysub1[1][y][z]*(Hx[1][y][z+1]-Hx[1][y][z]-Hz[2][y][z]+Hz[1][y][z])
- coef2ysub1[1][y][z]*(Exsub1[2][y][z]-Exsub1[2][y-1][z]-Exsub1[1][y][z]+Exsub1[1][y-1][z]);
triDiag[1][PCsize-3] = coef1ysub1[PCsize-2][y][z];
triDiag[0][PCsize-3] = -coef2ysub1[PCsize-2][y][z];
triDiag[3][PCsize-3] = coef3ysub1[PCsize-2][y][z]*Eysub1[PCsize-2][y][z]
+ coef4ysub1[PCsize-2][y][z]*(Hx[PCsize-2][y][z+1]-Hx[PCsize-2][y][z]-Hz[PCsize-1][y][z]+Hz[PCsize-2][y][z])
- coef2ysub1[PCsize-2][y][z]*(Exsub1[PCsize-1][y][z]-Exsub1[PCsize-1][y-1][z]-Exsub1[PCsize-2][y][z]+Exsub1[PCsize-2][y-1][z]);
solve(PCsize-2,triDiag); // Solve Tridiagonal Matrix
for (x=1;x<=PCsize-2;x++)
{
Eysub1[x][y][z]=triDiag[3][x-1]; // Set values in alternate array
}
}
for (z=1;z<=PCsize-2;z++)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -