/////////////////////////////////////////////
           //高階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;
}