antenna.cpp

计算了在有耗介质中半波天线的近场和原场的电场和磁场分量

#include "antenna.h"
//天线
int N;									  //天线每一臂分为N段
double h;                                 //每一小段的长度为h
double L;                                 //半波天线的臂长
double a;                                 //半波电线的半径
double lcoffie;							  //天线长度与波长的系数比

int m=0,n=0;

//电流元
doubleComplexMatrix *I;                   //电流元的全域系数按级数展开的系数
Mparameter Z;

//将所求的媒质参数赋给天线运算
void setting(Mparameter P)
{
	Z=P;
}	

//////////////////////////////////////////////////////////////////////
//                                  
//                           天       线
//
//////////////////////////////////////////////////////////////////////

//设置天线的电尺寸
void setAntenna(Mparameter P,double l,double radiu)   
{
	setting(P);

	lcoffie=l;
	L=lcoffie*Z.wavelength;
	a=radiu;
}	

double Antennaparameter()
{
	return L;
}

//天线单元剖分
void Antennaincisecell()                        
{
	int x;

	x=(int)10*lcoffie+1;

	if(x>5)  N=x;
	else     N=5;

	h=L/N;
	cout<<"Antenna"<<endl;
	cout<<"length  "<<L<<"   N  "<<N<<"   h "<<h/Z.wavelength<<endl<<endl;
}

////////////////////////////////////////////////////////////////////////
//
//                            海 伦 公 式
//
////////////////////////////////////////////////////////////////////////

//格林函数
doubleComplex greenfun(double z,double zi)
{
	double R,Ri;
	doubleComplex  x,y,greenfun,K;

	K=Z.K;
	
	R=sqrt(a*a+pow((z-zi),2));
	Ri=sqrt(a*a+pow((z+zi),2));
	
//	x = doubleComplexSet(K.im*R,-K.re*R);
	x.re = K.im*R;
	x.im = -K.re*R;
	x=doubleDMpy(1/(4*PI*R),expC(x));

//	y = doubleComplexSet(K.im*Ri,-K.re*Ri);
	y.re = K.im*Ri;
	y.im = -K.re*Ri;
	y=doubleDMpy(1/(4*PI*Ri),expC(y));

	greenfun = doubleCAdd(x,y);

	return greenfun;
}

//被积函数
doubleComplex cellfun(double zi)
{
	doubleComplex cellfun;
	doubleComplex y,zl;
	
	double zm,fabzi;

	if(m == N)
		zm=L;
	else 
		zm=(m+0.5)*h;

//	zl = pow((1-fabs(zi)/L),(n+1));

//	y=doubleDMpy((n+1)*(L-zi),Z.K);
	fabzi=(n+1)*(L-zi);
	y.re=fabzi*Z.K.re;
	y.im=fabzi*Z.K.im;

	zl=SinC(y);

	cellfun=doubleCMpy(zl,greenfun(zm,zi));
	
	return cellfun;
}

//积分总函数
doubleComplex integral()
{
	doubleComplex a,b,c,integral;
	if(m == 0)
	{
		a=integralcell(0,h);
		b=integralcell(h,L);
		integral=doubleCAdd(a,b);
	}
	else
		if(m == N)
		{
			a=integralcell(0,L-h/2);
			b=integralcell(L-h/2,L);
			integral=doubleCAdd(a,b);
		}
		else
		{			
			a=integralcell(0,m*h);
			b=integralcell(m*h,m*h+h);
			c=integralcell(m*h+h,L);
			integral=doubleCAdd(a,b);
			integral=doubleCAdd(integral,c);
		}

	return integral;
}

//辛普生法求积分
doubleComplex integralcell(double a,double b)
{
	return doubleCintegrat(a,b,cellfun);
}

//求解增广矩阵系数，即求电流全域基的系数I[][]
void Antsolvecoffie()
{
	Antennaincisecell();
	double zm;
	doubleComplex x,z,s;

	doubleComplexMatrix *B;
	doubleComplexMatrix *mat;
	doubleComplexMatrix *matInverse;
	
	B=doubleComplexMatrixAlloc(N+1,1);
	mat=doubleComplexMatrixAlloc(N+1,N+1);

	//x = doubleDMpy(Z.V/2,doubleCDiv(Z.ETA));
	z=doubleCDiv(Z.ETA);
	x.re=Z.V/2*z.re;
	x.im=Z.V/2*z.im;

	//增广矩阵系数的确定
	for(m=0;m<=N;m++)
	{
		if(m == N)
			zm=L;
		else 
			zm=(m+0.5)*h;

		for(n=0;n<N;n++)
		{
			mat->a[m][n]=integral();
		}
		mat->a[m][N]=CosC(doubleDMpy(zm,Z.K));        	      //与C相乘的系数
	
	//	z=doubleDMpy(zm,Z.K);
		z.re=zm*Z.K.re;
		z.im=zm*Z.K.im;

		s = doubleCMpy(x,SinC(z));
	//	B->a[m][0]=doubleComplexSet(s.im,-s.re);       //方程组右边B的值
		B->a[m][0].re = s.im;
		B->a[m][0].im = -s.re;
	}
		
	matInverse=doubleComplexMatrixInverse(mat);					//求A[][]的逆，~A[][]
	I=doubleComplexMatrixMpy(matInverse,B);						//电流元的全域基的系数与I=~A*B

/*	cout<<"matrix B"<<endl;
	doublePrintComplexMatrix(B);
	cout<<"matrix mat"<<endl;
	doublePrintComplexMatrix(mat);
	cout<<"matrix matInverse"<<endl;
	doublePrintComplexMatrix(matInverse);*/
	cout<<"matrix I"<<endl;
	doublePrintComplexMatrix(I);
	doubleComplexMatrixFree(mat);
	doubleComplexMatrixFree(matInverse);
}

////////////////////////////////////////////////////////////////////////
//
//                            线  电  流
//
////////////////////////////////////////////////////////////////////////

//画出线电流的分布情况
void PrintIline(FILE* fp)
{
	double IMO,z,x;
	doubleComplex I;
		
	for(int i=0;i<=20;i++)
	{
		z=i*L/20;
		I=Iline(z);
		x=i/20.0;
		z=x*lcoffie;
		IMO=doubleComplexAbs(I);

		if(i == 0|| i == 20)
		{
			cout<<"I("<<right<<setfill(' ')<<setw(10)<<x<<")          "<<I.re<<"+i"<<I.im<<endl;
		//	cout<<"I("<<right<<setfill(' ')<<setw(10)<<x<<")          "<<IMO<<endl;
		
			fprintf(fp,"I = ");
			fprintf(fp,"%e   ",I.re);
			fprintf(fp,"%e   ",I.im);
			fprintf(fp,"IMO = %e\n",IMO);
		}
		else
		{
			cout<<"I("<<left<<setfill('0')<<setw(10)<<x<<")          "<<I.re<<"+i"<<I.im<<endl;
		//	cout<<"I("<<left<<setfill('0')<<setw(10)<<x<<")          "<<IMO<<endl;
			fprintf(fp,"I = ");
			fprintf(fp,"%e   ",I.re);			
			fprintf(fp,"%e   ",I.im);
			fprintf(fp,"IMO = %e\n",IMO);
		}
	}
	cout<<endl;
}     
                                     
//线电流 
doubleComplex Iline(double zi)                              
{
	doubleComplex Iline,x,y,zl;
	double fabzi;

	Iline.re=0;
	Iline.im=0;

//	cout<<L-fabzi<<endl;
	//电流元的全域表达式
	for(int i=0;i<N;i++)
	{		
	//	zl=pow((1-fabzi/L),(i+1));
	//	y=doubleDMpy((i+1)*(L-fabs(zi)),Z.K);
		fabzi=(i+1)*(L-fabs(zi));
		y.re=fabzi*Z.K.re;
		y.im=fabzi*Z.K.im;
		
		zl=SinC(y);
		x=doubleCMpy(zl,I->a[i][0]); 
	
		//	cout<<"x  "<<x.re<<" + i"<<x.im<<endl;
		Iline=doubleCAdd(x,Iline);
	}
   
	return Iline;
}