fieldu.cpp

二维FDTD算法

/////////////////////////////////////////////
           //高阶FDTD二维程序
/////////////////////////////////////////////
#include "stdafx.h"
#include "stdlib.h"
#include "math.h"
#include <memory.h> 
#include "FieldU1.h"
#include "Def.h"
#include "PMLU.h"
#include "engine.h"


//迭代电磁场
void FieldIterate()
{
	//定义MATLAB变量
	Engine  *ep;			//引擎指针
	mxArray *show = NULL;   //显示矩阵
	double exchange[DOMAINX][DOMAINY];  //交换变量

	//打开MATLAB引擎
	if (!(ep = engOpen("\0"))) 
	{
		fprintf(stderr, "\nCan't start MATLAB engine\n");
		return;
	}
	else
	{
		printf("Already connect ot Matlab Engine, Please Wait!\n");
	}
	
	engEvalString(ep, "clear;clc");  //倾空内存
	show=mxCreateDoubleMatrix(DOMAINX, DOMAINY, mxREAL);  //分配空间
	/////////////////////////////////////////////////////////////////
	double Source=0.;        //源
	int    ipA=IC-18;    //测试点A点               
    int    jpA=JC;
    int    ipB=IC-18;    //测试点B点              
    int    jpB=JC-18;
	
    /////////////////////////////////////////////////////////////////
	//迭代主程序
	for (int time=0; time<T; time++)  //迭代次数(整数时间)
	{
	
			/**********************E平面波源的加入********************/
			Source=sin(2*PI*FREQUENCY*time*DELTAT);				//正弦
			if (T<=TIME_period/2)	
            {				   
			   Source=Source*0.5*(1-cos(PI*2*time/TIME_period)); //升余弦平滑
			}
			Source=exp(-0.5*pow((T0-time)/SPREAD,2.0)); //脉冲源

            /*************************E系数**************************/
			//H->E
			/**************************E入射源**************************/
			//IterateSourceEz(Coefficient1,Coefficient2);  //迭代一维源Ez_inc
			//GetEzValue(Source);							 //Ez赋值
			//PerfectAbsorbEz();                           //电场完美吸收

			/*********************E场************************/

			EzField(Coefficient1);  //Ez场
		    PMLEz(Coefficient1);    //Ez的PML
			//SourceEz(Coefficient1,Coefficient2); //Ez源		

			/*******************E软源条件**********************/
		    Ez[IC][JC]=Source;  
			
			/**********************H系数**********************/
			//E->H
			/*********************H源***********************/
			//IterateSourceHx(Coefficient1,Coefficient2); //迭代一维源Hx_inc
		   //PerfectAbsorbHx();                          //磁场完美吸收
			
			/**********************H场************************/
			HxField(Coefficient1);				//Hx场
			PMLHx(Coefficient1);               //Hx的PML
			//SourceHx(Coefficient1,Coefficient2);			//Hx源
		

			HyField(Coefficient1);				//Hy场
			PMLHy(Coefficient1);               //Hy的PML
			//SourceHy(Coefficient1,Coefficient2);			//Hy源
        /************************PML测试***************************/
            EzA[time]=Ez[ipA][jpA];
			EzB[time]=Ez[ipB][jpB];			
		
			if (maxA<=EzA[time])
			{
			   maxA=EzA[time];
			   maxB=EzB[time];
			}

		/**************************显示****************************/
		//赋值
		for (int i=0; i<DOMAINX; i++)
		{
			for (int j=0; j<DOMAINY; j++)
			{				
				exchange[i][j]=Ez[i][j];
			}
		}
		//拷贝变量
		memcpy((void *)mxGetPr(show), (void*)(exchange), DOMAINX*DOMAINY*sizeof(double));
		//放入matlab,命名
		engPutVariable(ep, "show", show);
		//3维场绘图
		//engEvalString(ep, "mesh(show),axis([0,120,0,120,-1,1]);");
		engEvalString(ep, "mesh(show);");
		//engEvalString(ep, "contour(show);");
		printf("第%d个时间步\n",time);
    }
}


//主程序
int main(int argc, char* argv[])
{
	InitialField();       //场初始化
	InitialPML();         //初始化PML
    FieldIterate();       //场迭代
    
	FILE *fp1,*fp2;
	fp1=fopen("EzA","w");
    fp2=fopen("EzB","w");
	for(int i=0;i<T;i++)
	{
	  fprintf(fp1,"%f ",EzA[i]);
      fprintf(fp2,"%f ",EzB[i]);	  
	}
	fclose(fp1);
	fclose(fp2);
	printf("迭代结束!\n");
    printf("max1= %f  max2= %f \n",Ez[23][45],Ez[45][23]);
	return 0;
}