📄 fdtd.cpp
字号:
/***********************
* FDTD.cpp
* Written by LZQ
* 2006.06.08
***********************/
#include "stdafx.h"
#include "FDTD.h"
#include <memory.h>
#include <string.h>
int ReadDataEpsilon() //从e_LC,mat文件中读取epsilon值
{
MATFile *fmat;
mxArray *pdata;
double *A;
const char *file = "e_LC.mat";
fmat = matOpen(file,"r");
if (NULL==fmat)
{
return(EXIT_FAILURE);
}
pdata = matGetArray(fmat,"thita1");
A = mxGetPr(pdata);
for (i=0;i<1801;i++)
{
for (j=0;j<251;j++)
{
thita[i][j]=*(A+i*251+j);
thita[i][j]=thita[i][j]*PI/180;
}
}
mxDestroyArray(pdata);
matClose(fmat);
for (i=0;i<1800;i++)
{
epsilonzxf[i] = 3.45*(cos(thita[i][0])*sin(thita[i][0])+cos(thita[i][1])*sin(thita[i][1])+cos(thita[i+1][0])*sin(thita[i+1][0])+cos(thita[i+1][1])*sin(thita[i+1][1]));
epsilonzzf[i] = 19.0+3.45*(pow(sin(thita[i][0]),2)+pow(sin(thita[i][1]),2)+pow(sin(thita[i+1][0]),2)+pow(sin(thita[i+1][1]),2));
epsilonzxb[i] = 3.45*(cos(thita[i][1799])*sin(thita[i][1799])+cos(thita[i][1798])*sin(thita[i][1798])+cos(thita[i+1][1798])*sin(thita[i+1][0])+cos(thita[i+1][1])*sin(thita[i+1][1]));
epsilonzzb[i] = 19.0+3.45*(pow(sin(thita[i][1799]),2)+pow(sin(thita[i][1798]),2)+pow(sin(thita[i+1][1799]),2)+pow(sin(thita[i+1][1798]),2));
}
for (k=0;k<250;k++)
{
epsilonxxl[k] = 19.0+3.45*(pow(cos(thita[0][k]),2)+pow(cos(thita[1][k]),2)+pow(cos(thita[0][k+1]),2)+pow(cos(thita[1][k+1]),2));
epsilonxxr[k] = 19.0+3.45*(pow(cos(thita[249][k]),2)+pow(cos(thita[248][k]),2)+pow(cos(thita[248][k+1]),2)+pow(cos(thita[249][k+1]),2));
epsilonxzl[k] = 3.45*(cos(thita[0][k])*sin(thita[0][k])+cos(thita[1][k])*sin(thita[1][k])+cos(thita[0][k+1])*sin(thita[0][k+1])+cos(thita[1][k+1])*sin(thita[1][k+1]));
epsilonxzr[k] = 3.45*(cos(thita[249][k])*sin(thita[249][k])+cos(thita[248][k])*sin(thita[248][k])+cos(thita[249][k+1])*sin(thita[249][k+1])+cos(thita[248][k+1])*sin(thita[248][k+1]));
}
return 0;
}
int WriteMatFile()
{
/************************************************************************/
/* 将矩阵输出到data.mat */
/************************************************************************/
MATFile *fmat;
mxArray *pdata1,*pdata2,*pdata3;
const char *file = "data.mat";
int status;
printf("\nCreating file %s...\n\n", file);
fmat = matOpen(file,"w");
if (fmat == NULL)
{
printf("Error creating file %s\n", file);
return(EXIT_FAILURE);
}
pdata1 = mxCreateDoubleMatrix(1956,406,mxREAL);
memcpy((void *)(mxGetPr(pdata1)), (void *)Ex_final, sizeof(Ex_final));
status = matPutVariable(fmat, "Ex", pdata1);
if (status != 0)
{
printf("%s : Error using matPutVariable on line %d\n", __FILE__, __LINE__);
return(EXIT_FAILURE);
}
pdata2 = mxCreateDoubleMatrix(1956,406,mxREAL);
memcpy((void *)(mxGetPr(pdata2)), (void *)Ez_final, sizeof(Ez_final));
status = matPutVariable(fmat, "Ez", pdata2);
if (status != 0)
{
printf("%s : Error using matPutVariable on line %d\n", __FILE__, __LINE__);
return(EXIT_FAILURE);
}
pdata3 = mxCreateDoubleMatrix(1956,406,mxREAL);
memcpy((void *)(mxGetPr(pdata3)), (void *)Hy_final, sizeof(Hy_final));
status = matPutVariable(fmat, "Hy", pdata3);
if (status != 0)
{
printf("%s : Error using matPutVariable on line %d\n", __FILE__, __LINE__);
return(EXIT_FAILURE);
}
mxDestroyArray(pdata1);
mxDestroyArray(pdata2);
mxDestroyArray(pdata3);
if (matClose(fmat) != 0)
{
printf("Error closing file %s\n",file);
return(EXIT_FAILURE);
}
return(EXIT_SUCCESS);
}
void ShowError(char *msg)
{
MessageBox(NULL,msg,"ERROR",0);
//printf("\n%s failed:%d",GetLastError());
}
DWORD WINAPI ShowProInfo()
{
float m;
CONSOLE_SCREEN_BUFFER_INFO ConsoleScreenBufferInfo;
hStdout=GetStdHandle(STD_OUTPUT_HANDLE);
if(hStdout==INVALID_HANDLE_VALUE)
{
ShowError("GetStdHandle");
return(EXIT_FAILURE);
}
if(!GetConsoleScreenBufferInfo(hStdout,&ConsoleScreenBufferInfo))
{
ShowError("GetConsoleScreenBufferInfo");
return(EXIT_FAILURE);
}
ConsoleScreenBufferInfo.dwCursorPosition.X = 0;
SetConsoleCursorPosition(hStdout,ConsoleScreenBufferInfo.dwCursorPosition);
m = (float)n*100/steps;
printf("已完成总计算量的--------------------------------------------------%6.0f%%",m);
return(EXIT_SUCCESS);
}
int main(int argc,char **argv)
{
T = 0;
steps = 1;
for (i=0;i<NUM_EX;i++)
{
for (k=0;k<NUM_EZ;k++)
{
Ex_pre[i][k] = 0.0;
Ez_pre[i][k] = 0.0;
Hy_pre[i][k] = 0.0;
Hyx_pre[i][k] = 0.0;
Hyz_pre[i][k] = 0.0;
}
}
if (ReadDataEpsilon()==1)
{
printf("Error openning file!\n");
return(EXIT_FAILURE);
}
//PML介质参数
for (i=0;i<9;i++)
{
deltax[i] = 1.25*EPSILON0*C0/dx*pow((9.0-i)/9,4);
deltaz[i] = 1.25*EPSILON0*C0/dz*pow((9.0-i)/9,4);
deltamx[i] = 1.25*MU0*C0/dx*pow((9.0-i)/9,4);
deltamz[i] = 1.25*MU0*C0/dz*pow((9.0-i)/9,4);
}
//开始FDTD主循环
//while (steps>0)
//{
printf("\n输入计算所需的时间步数--->");
scanf("%d",&steps);
printf("%d \n",steps);
for (n=1;n<=steps;n++) //时间步长循环
{
T = T+1;
/************************************************************************/
/* PML吸收边界 (8) */
/************************************************************************/
//下边界
for (i=8;i<1948;i++)
{
Hyx_nxt[i][8] = Hyx_pre[i][8]+0.5/Z0*(Ez_pre[i+1][8]-Ez_pre[i][8]);
Hyz_nxt[i][8] = Hyz_pre[i][8]*exp(-1*deltamz[8]*dt/MU0)-(1-exp(-1*deltamz[8]*dt/MU0))/dz/deltamz[8]*(Ex_pre[i][9]-Ex_pre[i][8]);
Hy_nxt[i][8] = Hyx_nxt[i][8]+Hyz_nxt[i][8];
Ex_nxt[i][8] = Ex_pre[i][8]-0.5*Z0*(Hy_nxt[i][8]-Hy_nxt[i][7]);
Ez_nxt[i][8] = Ez_pre[i][8]+0.5*Z0*(Hy_nxt[i][8]-Hy_nxt[i-1][8]);
}
//上边界
for (i=8;i<1948;i++)
{
Hyx_nxt[i][398] = Hyx_pre[i][398]+0.5/Z0*(Ez_pre[i+1][398]-Ez_pre[i][398]);
Hyz_nxt[i][398] = Hyz_pre[i][398]*exp(-1*deltamz[8]*dt/MU0)-(1-exp(-1*deltamz[8]*dt/MU0))/dz/deltamz[8]*(Ex_pre[i][399]-Ex_pre[i][398]);
Hy_nxt[i][398] = Hyx_nxt[i][398]+Hyz_nxt[i][398];
Ex_nxt[i][398] = Ex_pre[i][398]-0.5*Z0*(Hy_nxt[i][398]-Hy_nxt[i][397]);
Ez_nxt[i][398] = Ez_pre[i][398]+0.5*Z0*(Hy_nxt[i][398]-Hy_nxt[i-1][398]);
}
//左边界
for (k=8;k<398;k++)
{
Hyx_nxt[8][k] = Hyx_pre[8][k]*exp(-1*deltamx[8]*dt/MU0)+(1-exp(-1*deltamx[8]*dt/MU0))/dx/deltamx[8]*(Ez_pre[9][k]-Ez_pre[8][k]);
Hyz_nxt[8][k] = Hyz_pre[8][k]-0.5/Z0*(Ex_pre[8][k+1]-Ex_pre[8][k]);
Hy_nxt[8][k] = Hyx_nxt[8][k]+Hyz_nxt[8][k];
Ex_nxt[8][k] = Ex_pre[8][k]-0.5*Z0*(Hy_nxt[8][k]-Hy_nxt[8][k-1]);
Ez_nxt[8][k] = Ez_pre[8][k]+0.5*Z0*(Hy_nxt[8][k]-Hy_nxt[7][k]);
}
//右边界
for (k=8;k<398;k++)
{
Hyx_nxt[1948][k] = Hyx_pre[1948][k]*exp(-1*deltamx[8]*dt/MU0)+(1-exp(-1*deltamx[8]*dt/MU0))/dx/deltamx[8]*(Ez_pre[1949][k]-Ez_pre[1948][k]);
Hyz_nxt[1948][k] = Hyz_pre[1948][k]-0.5/Z0*(Ex_pre[1948][k+1]-Ex_pre[1948][k]);
Hy_nxt[1948][k] = Hyx_nxt[1948][k]+Hyz_nxt[1948][k];
Ex_nxt[1948][k] = Ex_pre[1948][k]-0.5*Z0*(Hy_nxt[1948][k]-Hy_nxt[1948][k-1]);
Ez_nxt[1948][k] = Ez_pre[1948][k]+0.5*Z0*(Hy_nxt[1948][k]-Hy_nxt[1947][k]);
}
/************************************************************************/
/* PML内部层(1-7) */
/************************************************************************/
//下边
for (i=8;i<1498;i++)
{
for (k=1;k<8;k++)
{
Hyx_nxt[i][k] = Hyx_pre[i][k]+0.5/Z0*(Ez_pre[i+1][k]-Ez_pre[i][k]);
Hyz_nxt[i][k] = Hyz_pre[i][k]*exp(-1*deltamz[k]*dt/MU0)-(1-exp(-1*deltamz[k]*dt/MU0))/dz/deltamz[k]*(Ex_pre[i][k+1]-Ex_pre[i][k]);
Hy_nxt[i][k] = Hyz_nxt[i][k]+Hyz_nxt[i][k];
Ex_nxt[i][k] = Ex_pre[i][k]*exp(-1*deltaz[k]*dt/EPSILON0)-(1-exp(-1*deltaz[k]*dt/EPSILON0))/dz/deltaz[k]*(Hy_nxt[i][k]-Hy_nxt[i][k-1]);
Ez_nxt[i][k] = Ez_pre[i][k]+0.5*Z0*(Hy_nxt[i][k]-Hy_nxt[i-1][k]);
}
}
//上边
for (i=8;i<1498;i++)
{
for (k=399;k<406;k++)
{
Hyx_nxt[i][k] = Hyx_pre[i][k]+0.5/Z0*(Ez_pre[i+1][k]-Ez_pre[i][k]);
Hyz_nxt[i][k] = Hyz_pre[i][k]*exp(-1*deltamz[k-398]*dt/MU0)-(1-exp(-1*deltamz[k-398]*dt/MU0))/dz/deltamz[k-398]*(Ex_pre[i][k+1]-Ex_pre[i][k]);
Hy_nxt[i][k] = Hyz_nxt[i][k]+Hyz_nxt[i][k];
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -