fdtd.h

等离子体覆盖目标的电磁特性的FDTD分析

#include<stdio.h>
#include<math.h>
#include<fstream.h>
#include<iomanip.h>


const double	pii=3.1415926535898;
const double	epsn0=1e-9/(36.*pii);	
const double	miu0=1e-7*4.*pii;		
const double	Ltspd=3.0e+8;			
const double	yta0=120*pii;			
const double	relpt=1.0;				
const double	relpb=1.0;				
const double	epsn=relpt*epsn0;
const double	miu=relpb*miu0;
const double	fr=8.0e11;	
						


const double	dz=1.5e-5;//0.125;//2.0e-4;
const double	dt=2.5e-14;//0.208e-9;//0.5e-12;

const int   MZ=4500;			
const int   Nth=5;			
const int	K1=Nth;
const int	K2=MZ-Nth;

double etx1[Nth+1],ety1[Nth+1],etz1[Nth+1],etx2[Nth+1],ety2[Nth+1],etz2[Nth+1];
double mtx1[Nth+1],mty1[Nth+1],mtz1[Nth+1],mtx2[Nth+1],mty2[Nth+1],mtz2[Nth+1];

const double		tem1=dt/epsn;
const double		tem2=dt/miu;
const double		tem3=miu/epsn;
const double		g=1.44;	//the ratio of the geometric progression
const double		R0=0.001;	
const double		ct0=epsn*Ltspd*log(R0)/(2.0*(pow(g,(double)Nth)-1.0));
const double		ce0=ct0*(1.0-sqrt(g));
const double		cm0=tem3*ce0;

const int npml=5;
const double sigmmax=0.0-log(R0)*(npml+1)/2.0*epsn0*Ltspd/Nth;
const double sigma0=sigmmax/((npml+1)*pow(2,npml+1)*pow(Nth,npml));



const int thick=2000;			
const int SN=3;
const int low=MZ/2-thick;
const int high=MZ/2+thick;
const int kk0=low-SN;		
const int kk1=high+SN;		
int index;
double EXi[MZ+1],HYi[MZ];	


double inline conduct(int L){
	if (L<=0) return sigma0;
	else return sigma0*(pow(2.0*L+1.0,npml+1)-pow(2.0*L-1.0,npml+1));
}

double inline mconduct(int  L){
	if (L<=0) return tem3*sigma0;
	else return tem3*sigma0*(pow(2.0*(L+0.5)+1.0,npml+1)-pow(2.0*(L+0.5)-1.0,npml+1));
}



double sch(int t, const double aa=2.0){
	double fr=1e10;
	double enevop=1.0,w=2.0*pii*fr;
	int    T_max=(int)(0.5+aa/fr/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 scp(int t){
	double tao=50*dt;
    double to=2.25/fr; 
	return  -1.0*cos(2.0*pii*fr*t*dt)*exp(0.0-4.0*pii*(t*dt-to)*(t*dt-to)/(25.0*tao*tao));

}



void Zeroinit(void){
	int i,index;

	for(i=0;i<Nth+1;i++) etz1[i]=exp(0.0-tem1*conduct(i)/dz);
	for(i=0;i<Nth+1;i++) mtz1[i]=exp(0.0-tem2*mconduct(i)/dz);
	for(i=0;i<Nth+1;i++) etz2[i]=(1.0-exp(0.0-tem1*conduct(i)/dz))/conduct(i)*dz;
	for(i=0;i<Nth+1;i++) mtz2[i]=(1.0-exp(0.0-tem2*mconduct(i)/dz))/mconduct(i)*dz;
	for(index=0;index<=MZ;index++) EXi[index]=0.0;
	for(index=0;index<MZ;index++) HYi[index]=0.0;

 }