source.h

基于电磁波方程

void AddSourceToFDTD(int flag,int n,int is=10,int js=10)
{
	float temp;
	switch(flag)
	{
	case 0:
		temp=10*sin(OMIGA*n);
		if(n<=WL) temp*=0.5*(1.0-cos(OMIGA*n));
		break;
	case 1:
		temp=1000*(6-0.083*n)*expf(-(6-0.083*n)*(6-0.083*n)/4.0);
		break;
	case 2:
		{
			float pi=3.1415926;//设敏感介电常数为9，v*dt=0.5dx，v=1e8(m/s),设置dx=1/20波长
			float freq=900e6;//而最大频率为2700MHz，而波长为1/27（=0.037m),
			float dt=2.25e-12;//则dx=0.00135(m),则dt=2.25e-12(s)
			float rickth=2.0*pi*pi*freq*freq;
			float rickper=1.0/freq;
			float time=n*dt;
			float ricksc=sqrt(exp(1.0)/(2.0*rickth));
			float delay=(time-rickper);
			float tempRicker=exp(-rickth*delay*delay);
			temp=1000*ricksc*tempRicker*(-2.0)*rickth*delay;
		}
		break;
	default:
		break;
	}
	Ez[is][js]-=temp;
}

float InE(float i,float j)
{
	int II,III;
	float t1;
	t1=i*cp+j*sp-offset;
	II=int(t1);
	if(t1>=0.0)
	{
		III=II+1;
		t1=float(III)-t1;
	}
	else
	{
		III=II-1;
		t1=t1-float(III);
	}
	return t1*(Ein[II]-Ein[III])+Ein[III];
}
//----------------------------
float InH(float i,float j)
{
	int II,III;
	float t1;
	t1=i*cp+j*sp-offset;
	if(t1>0.5)
	{
		II=int(t1-0.5);
		III=II+1;
		t1=float(III)+0.5-t1;
	}
	else
	{
		II=int(t1+0.5);
		III=II-1;
		t1=t1-0.5-float(III);
	}
	return t1*(Hin[II]-Hin[III])+Hin[III];
}
//-------------------------------------------------------
void AddSqareWave(int flag,int N)
{
	//设置入射波
	
	float temp;
	int InStart,InEnd,ISource;
	InStart=0,InEnd=900,ISource=20;
	for(int i=InStart;i<=InEnd-1;i++)
	{
		Hin[i]=Hin[i]-0.5*(Ein[i+1]-Ein[i]);
	}
	for(i=InStart+1;i<=InEnd-1;i++)
	{
		Ein[i]=Ein[i]+0.5*(Hin[i-1]-Hin[i]);
	}
	Ein[InStart]=EBin[3];
	EBin[3]=EBin[2];
	EBin[2]=Ein[InStart+1];
	Ein[InEnd]=EBin[1];
	EBin[1]=EBin[0];
	EBin[0]=Ein[InEnd-1];
	
	switch(flag)
	{
	case 0:
		temp=10*sin(OMIGA*N);
		if(N<=WL) temp*=0.5*(1.0-cos(OMIGA*N));
		break;
	case 1:
		temp=1000*(6-0.083*N)*expf(-(6-0.083*N)*(6-0.083*N)/4.0);
		break;
	case 2:
		{
			float pi=3.1415926;//设敏感介电常数为9，v*dt=0.5dx，v=1e8(m/s),设置dx=1/20波长
			float freq=900e6;//而最大频率为2700MHz，而波长为1/27（=0.037m),
			float dt=2.25e-12;//则dx=0.00135(m),则dt=2.25e-12(s)
			float rickth=2.0*pi*pi*freq*freq;
			float rickper=1.0/freq;
			float time=N*dt;
			float ricksc=sqrt(exp(1.0)/(2.0*rickth));
			float delay=(time-rickper);
			float tempRicker=exp(-rickth*delay*delay);
			temp=1000*ricksc*tempRicker*(-2.0)*rickth*delay;
		}
		break;
	default:
		temp=1000*(6-0.083*N)*expf(-(6-0.083*N)*(6-0.083*N)/4.0);
		break;
	}
	
	Ein[ISource]=temp;

	
}
void IntroduceInE()
{
	//引入入射波E
	for(int i=ItMin+1;i<=ItMax-1;i++)
	{
		Ez[i][JtMax]-=0.5*sp*InH(i,JtMax+0.5);
		Ez[i][JtMin]+=0.5*sp*InH(i,JtMin-0.5);
	}
	for(int j=JtMin+1;j<=JtMax-1;j++)
	{
		Ez[ItMax][j]-=0.5*cp*InH(ItMax+0.5,j);
		Ez[ItMin][j]+=0.5*cp*InH(ItMin-0.5,j);
	}//角点的处理
	Ez[ItMax][JtMax]+=0.5*(-sp*InH(ItMax,JtMax+0.5)-cp*InH(ItMax+0.5,JtMax));
	Ez[ItMin][JtMin]+=0.5*(sp*InH(ItMin,JtMin-0.5)+cp*InH(ItMin-0.5,JtMin));
	Ez[ItMin][JtMax]+=0.5*(-sp*InH(ItMin,JtMax+0.5)+cp*InH(ItMin-0.5,JtMax));
	Ez[ItMax][JtMin]+=0.5*(sp*InH(ItMax,JtMin-0.5)-cp*InH(ItMax+0.5,JtMin));
}
//---------------------------------------------
void IntroduceInH()
{
	for(int i=ItMin;i<=ItMax;i++)
	{
		Hx[i][JtMax]-=0.5*InE(i,JtMax);
		Hx[i][JtMin-1]+=0.5*InE(i,JtMin);
	}
	for(int j=JtMin;j<=JtMax;j++)
	{
		Hy[ItMax][j]+=0.5*InE(ItMax,j);
		Hy[ItMin-1][j]-=0.5*InE(ItMin,j);
	}
}

/*
float gprmaxso(char type[],float amp,float freq,float dt,float total_time)
{
	if(isstr(type)!=1)
	printf('First argument should be a source type');
}

RAMPD=0.25;
if(nargin < 5)
error('GPRMAXSO requires all five arguments ');
end;

if(isstr(type)~=1)
error('First argument should be a source type');
end;
if(freq==0)
error(['Frequency is zero']);
end;
iter=total_time/dt;
time=0;


if(strcmp(type,'ricker')==1) 
rickth=2.0*pi*pi*freq*freq;
rickper=1.0/freq;
ricksc=sqrt(exp(1.0)/(2.0*rickth));
i=1;
while(time<=total_time)
 delay=(time-rickper);
 temp=exp(-rickth*delay*delay);
 excitation(i)=ricksc*temp*(-2.0)*rickth*delay;
 time=time+dt;
 i=i+1;
 end;
end;

if(strcmp(type,'gaussian')==1)  
rickper=1.0/freq;
rickth=2.0*pi*pi*freq*freq;
i=1;
while(time<=total_time)
 delay=(time-rickper);
 excitation(i)=exp((-rickth)*delay*delay);
 time=time+dt;
 i=i+1;
 end;
end;


if(strcmp(type,'cont_sine')==1)  
i=1;
while(time<=total_time)
 ramp=time*RAMPD*freq;
  if(ramp>1.0) 
   ramp=1.0;
   end;
 excitation(i)=ramp*sin(2.0*pi*freq*time);
 time=time+dt;
 i=i+1;
 end;
end;

if(strcmp(type,'sine')==1)  
i=1;
while(time<=total_time)
 excitation(i)=sin(2.0*pi*freq*time);
 if(time*freq>1.0) 
   excitation(i)=0.0;
   end;
 time=time+dt;
 i=i+1;
 end;
end;

excitation=excitation.*amp;
*/