parameters.h

2维的FDTD C程序

#include<stdio.h>
#include<math.h>
#include<float.h>
#include<stdlib.h>
#include<conio.h>
#include<iomanip.h>
#include<afx.h>
#include<iostream.h>
#include<fstream.h>
#include<memory.h>

const double	pii=3.1415926535898;
const double	Epsilon0=1e-9/(36.*pii);	//Permit in free space
const double	Mu0=1e-7*4.*pii;			//Permeab in free space
const double	LightSpeed=3.0e+8;			//Light Speed in free space(m/s)
const double	Eta0=120*pii;				//wave impendance in free space(Omn)
const double	Epsilon_r=1.0;				//Epsilon/Epsilon0
const double	Mu_r=1.0;					//Mu/Mu0
const double	Epsilon=Epsilon_r*Epsilon0;		//air space
const double	Mu=Mu_r*Mu0;

const double	Frequency=3.0e9;							//frequency(Hz)
const double	Wavelength=LightSpeed/Frequency;			//wavelength(0.1m)

//for fdtd: 2 dimension cartesian coordinates
const double	dx=Wavelength/140.0;//2.5e-4;
const double	dy=Wavelength/140.0;//2.5e-4;//
const double	dt=0.5*__min(dx,dy)/LightSpeed;
const double	tao=6e-11/dt;		

const int		Target_X=18;							//the X/2 thickness of scattering target
const int		Target_Y=18;								//the Y/2 thickness of scattering target
const int		rPlasma=0;//70;								    //the plasma thickness 
const int		PML_Thickness=5;							//the thickness of PML
const int		Interval_PML_scatter=5;					//the interval between PML and total-scattered field boundary 
const int		Interval_scatter_target=5;					//the interval between total-scattered field boundary and target

const int		MX=2*(Target_X+PML_Thickness+Interval_scatter_target+Interval_PML_scatter+rPlasma); 
const int		MY=2*(Target_Y+PML_Thickness+Interval_scatter_target+Interval_PML_scatter+rPlasma);		//y axis maxium length(in deltas)		
const int		dia=int(sqrt(MX*MX+MY*MY))+2*PML_Thickness+20;
const int		ii0=int(MX/2)-(Target_X+Interval_scatter_target);  //total-scattered field boundary
const int		ii1=int(MX/2)+(Target_X+Interval_scatter_target);
const int		jj0=int(MY/2)-(Target_Y+Interval_scatter_target);
const int		jj1=int(MY/2)+(Target_Y+Interval_scatter_target);
	
//for PML in 2 dimension cartesian coordinates
const int		I1=PML_Thickness;
const int		I2=MX-PML_Thickness;
const int		J1=PML_Thickness;
const int		J2=MY-PML_Thickness;
const double	dte=dt/Epsilon;			//	air 
const double	dtm=dt/Mu;
const double	me=Mu/Epsilon;

const double	R0=0.001;	//the reflect coefficent on the ObjectSpace-PML interface
const int		index=3;	//variant index for Berenger's PML
const double	sigma_max=0.0-log(R0)*(index+1)/2.0*Epsilon*LightSpeed/PML_Thickness/__min(dx,dy);

double inline conduct(double L){
	if (L<=0) return 0.0;
	else return sigma_max*pow(L/PML_Thickness,index);
}
double inline mconduct(double  L){
	return me*conduct(L+0.5);
}


//	the part for exciting source

double scp(int t){
	if(t<=0)  return 0.0;
	else
		return exp(0.0-(t-5.0*tao)*(t-5.0*tao)/(tao*tao/2.0));//tao=1e-10/dt
}

double scp0(int t){
	if(0<=t&&t<=10*tao) 
		return -1000*(t-5.0*tao)*exp(0.0-(t-5.0*tao)*(t-5.0*tao)/(tao*tao/2.0));
	else return(0.0);
}

double sch(int t, const double aa=2.0){
	double enevop=1.0,w=2.0*pii*Frequency;
	int    T_max=(int)(0.5+aa/Frequency/dt);
	if (t<0) enevop=0.0;
	else if ((0<=t)&&(t<=T_max)) enevop=(0.5*(1.0-cos(w*t*dt*0.5/aa)));
		else enevop=1.0;		
	return(enevop*sin(w*t*dt));
}

double sins(int t){
	if(t<=0)  return 0.0;
	else      return(sin(2.0*pii*Frequency*t*dt));
}

//	the part for exciting source end


int inline _cdecl round(double x){

	int xx=(int)x;

	if(x>0.0){
		if(x-xx>=0.5)
			return (xx+1);
		else
			return xx;
	}
	else{
		if(fabs(x-xx)>=0.5)
			return(xx-1);
		else
			return xx;
	}
}