matrixinverse.cpp

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

#include "MatrixInverse.h"

/*double complex vector */

doubleComplexMatrix *doubleComplexMatrixAlloc(int row,int col)
{
	int i;
	doubleComplexMatrix *r;

	r=(doubleComplexMatrix *)malloc(sizeof(doubleComplexMatrix));
	r->a=(doubleComplex **)malloc(row*sizeof(doubleComplex *));
	for(i=0;i<row;i++)
		r->a[i]=(doubleComplex *)malloc(col*sizeof(doubleComplex));
	r->row=row;
	r->col=col;

	return r;
}

void doubleComplexMatrixFree(doubleComplexMatrix *mat)
{
	int i;

	for(i=0;i<mat->row;i++)
		free(mat->a[i]);

	free(mat->a);
	free(mat);
}
void doublePrintComplexMatrix(doubleComplexMatrix *mat)
{
	int i,j;

	for(i=0;i<mat->row;i++)
	{
		for(j=0;j<mat->col;j++)
		{
			printf("%e+i(%e)  ",mat->a[i][j].re,mat->a[i][j].im);
		}
		printf("\n");
	}
}
/*
*Multiplication for complex matrix
*in)two matrixs:m*p,p*n
*out)matrix:m*n
*/
doubleComplexMatrix *doubleComplexMatrixMpy(doubleComplexMatrix *mat1,doubleComplexMatrix *mat2)
{
	int n;
	int i,j,k;
	int row,col;
	doubleComplexMatrix *r;
	doubleComplex temp;
	
	row=mat1->row;
	col=mat2->col;
	n=mat1->col;
	r=doubleComplexMatrixAlloc(row,col);
	
	for(i=0;i<row;i++)
	{
		for(j=0;j<col;j++)
		{
			temp=doubleComplexSet(0,0);
			for(k=0;k<n;k++)
				temp=doubleCAdd(temp,doubleCMpy(mat1->a[i][k],mat2->a[k][j]));
			r->a[i][j]=temp;
		}
	}
	
	return r;
}
/*
* Get the Inverse matrix of a matrix
*/
doubleComplexMatrix *doubleComplexMatrixInverse(doubleComplexMatrix *mat)
{
	int n;
	int i,j,k;
	double Abs,max;
	doubleComplexMatrix *p;
	doubleComplexMatrix *q;
	doubleComplex *rowTemp,*rowTemp1,ptemp;
	
	n=mat->col;
	p=doubleComplexMatrixAlloc(n,n);/* used for calculating the temp matrix after several kinds of preliminary transform */
	q=doubleComplexMatrixAlloc(n,n);/* used for calculating the inverse */
	rowTemp=(doubleComplex *)malloc(n*sizeof(doubleComplex));
	rowTemp1=(doubleComplex *)malloc(n*sizeof(doubleComplex));
	
	for(i=0;i<n;i++)
	{	
		for(j=0;j<n;j++)
		{
			p->a[i][j]=mat->a[i][j];
			if(i==j)
				q->a[i][j]=doubleComplexSet(1,0);
			else
				q->a[i][j]=doubleComplexSet(0,0);
		}
	}	
	for(k=0;k<n;k++)
	{//printf("k  %d\n",k);
		//行主元数为0，则与同列中不为0的行，进行行互换
		Abs=doubleComplexAbs(p->a[k][k]);
		if((Abs==0)||(Abs<1E-3))
		{
			for(i=k,max=0,j=k+1;j<n;j++)
			{
				Abs=doubleComplexAbs(p->a[j][k]);
				if((Abs>max) && (Abs>1E-3))
				{
					max=Abs;
					i=j;
				}
				if(max==0)
					return 0;
			}
			for(j=k;j<n;j++)
			{
				rowTemp[j]=p->a[k][j];
				p->a[k][j]=p->a[i][j];
				p->a[i][j]=rowTemp[j];
				
				rowTemp[j]=q->a[k][j];
				q->a[k][j]=q->a[i][j];
				q->a[i][j]=rowTemp[j];
			}
		}
		
		//将k行各元素除以该行主元数，即规一化
		ptemp=doubleCDiv(p->a[k][k]);
		
		for(j=0;j<k;j++)
		{
			q->a[k][j]=doubleCMpy(q->a[k][j],ptemp);
		}
		for(j=k;j<n;j++)
		{
			p->a[k][j]=doubleCMpy(p->a[k][j],ptemp);
			q->a[k][j]=doubleCMpy(q->a[k][j],ptemp);
		}
		
		for(i=0;i<k;i++)
		{
			rowTemp1[i]=q->a[k][i];
		}
		for(i=k;i<n;i++)
		{
			rowTemp[i]=p->a[k][i];
			rowTemp1[i]=q->a[k][i];
		}
		for(i=0;i<n;i++)
		{
			if(i!=k)
			{
				ptemp=p->a[i][k];
				for(j=0;j<k;j++)
				{
					q->a[i][j]=doubleCSub(q->a[i][j],doubleCMpy(ptemp,rowTemp1[j]));
				//	q->a[i][j]=doubleCMpy(rowTemp1[j],doubleCSub(q->a[i][j],ptemp));
				}
				for(j=k;j<n;j++)
				{
					p->a[i][j]=doubleCSub(p->a[i][j],doubleCMpy(ptemp,rowTemp[j]));
					q->a[i][j]=doubleCSub(q->a[i][j],doubleCMpy(ptemp,rowTemp1[j]));
				//	p->a[i][j]=doubleCMpy(rowTemp[j],doubleCSub(p->a[i][j],ptemp));
				//	q->a[i][j]=doubleCMpy(rowTemp1[j],doubleCSub(q->a[i][j],ptemp));
				}
			}
		}	
	}
	doubleComplexMatrixFree(p);
	free(rowTemp);
	free(rowTemp1);
	return q;
}

doubleComplex doubleComplexSet(double real,double imag)
{
	doubleComplex a;
	a.re=real;
	a.im=imag;
	return a;
}
/*
* Get the modulus of a complex
*/
double doubleComplexAbs(doubleComplex x)
{
	double r;
	r=sqrt(x.re*x.re+x.im*x.im);
	return r;
}
/* ADDITION FOR COMPLEX
* parameters:
* in)x--input complex x
* y--intput complex y
* out)z--output complex z
*/

doubleComplex doubleCAdd(doubleComplex x,doubleComplex y)
{
	doubleComplex z;
	z.re=x.re+y.re;
	z.im=x.im+y.im;
	return z;
}

doubleComplex doubleCSub(doubleComplex x,doubleComplex y)
{
	doubleComplex z;
	z.re=x.re-y.re;
	z.im=x.im-y.im;
	return z;
}
/*
* MULTIPLICATION FOR COMPLEX
* parameters:
* in)x--input complex x
* y--intput complex y
* out)z--output complex z
*/

doubleComplex doubleCMpy(doubleComplex x,doubleComplex y)
{
	doubleComplex z;
	z.re=x.re*y.re-x.im*y.im;
	z.im=x.re*y.im+x.im*y.re;
	return z;
}

doubleComplex doubleCDiv(doubleComplex x)
{
	doubleComplex z;
	double r;
	
	r=x.re*x.re+x.im*x.im;
	z.re=x.re/r;
	z.im=-x.im/r;
	
	return z;
}

doubleComplex doubleDMpy(double c,doubleComplex y)
{
	doubleComplex z;
	
	z.re=c*y.re;
	z.im=c*y.im;
	
	return z;
}

//对复数开根号
doubleComplex sqrtC(doubleComplex z)
{
	doubleComplex sqrtC;
	double x,y,r;
	double re,im;
	
	x=z.re;
	y=z.im;
	r=sqrt(x*x + y*y);
	
	re=sqrt((r + x)/2);
	im=sqrt((r - x)/2);
	
	sqrtC.re = re;
	
	if(y > 0)
		sqrtC.im = im;
	else 
		if(y < 0)
			sqrtC.im = - im;
		else
		{
			if(x > 0)
				sqrtC.im = 0;
			else
			{
				sqrtC.re = 0;
				sqrtC.im = im;
			}
		}
		
		return sqrtC;
}

doubleComplex expC(doubleComplex z)
{
	doubleComplex expC;
	double x,y;
	
	x=z.re;
	y=z.im;
	expC.re=exp(x)*cos(y);
	expC.im=exp(x)*sin(y);
	
	return expC;
}

doubleComplex SinhC(doubleComplex z)
{
	doubleComplex a,b,SinhC;
	
	a=expC(z);
	b=doubleDMpy(-1,expC(doubleDMpy(-1,z)));
	z=doubleCAdd(a,b);
	SinhC=doubleDMpy(0.5,z);
	
	return SinhC;
}
doubleComplex CoshC(doubleComplex z)
{
	doubleComplex a,b,CoshC;
	
	a=expC(z);
	b=expC(doubleDMpy(-1,z));
	z=doubleCAdd(a,b);
	CoshC=doubleDMpy(0.5,z);
	
	return CoshC;
}

doubleComplex TanhC(doubleComplex z)
{
	doubleComplex TanhC;
	
	TanhC=doubleCMpy(SinhC(z),doubleCDiv(CoshC(z)));
	
	return TanhC;
}

doubleComplex SinC(doubleComplex z)
{
	doubleComplex x,a,b,SinC;
	
	x.re = -z.im;
	x.im = z.re;
	a=expC(x);
	b=doubleDMpy(-1,expC(doubleDMpy(-1,x)));
	x=doubleCAdd(a,b);
	x=doubleDMpy(0.5,x);
	SinC.re = x.im;
	SinC.im = -x.re;
	
	return SinC;
}
doubleComplex CosC(doubleComplex z)
{
	doubleComplex x,a,b,CosC;
	
	x.re = -z.im;
	x.im = z.re;
	a=expC(x);
	b=expC(doubleDMpy(-1,x));
	x=doubleCAdd(a,b);
	CosC=doubleDMpy(0.5,x);
	
	return CosC;
}

void doubleCprint(doubleComplex z)
{
	printf("%e+i(%e)  \n",z.re,z.im);
}

doubleComplex doubleCintegrat(double a,double b,doubleComplex (*f)(double zi))
{
	int n = 1;          //赋初值
	double D;
	doubleComplex Tn,T2n,In,I2n,integrat;
	
	D = b - a;
	
	T2n.re = I2n.re = D * (f(a).re + f(b).re)/2;
	Tn.re = 0;
	T2n.im = I2n.im = D * (f(a).im + f(b).im)/2;
	Tn.im = 0;
	
	double sigma;
	double siga;
	//积分运算函数
	while (( fabs(I2n.re - In.re) >= eps)&&(fabs(I2n.im - In.im) >= eps))
	{
		sigma = 0.0;
		siga = 0.0;
		
		//实部
		double x;
		Tn.re = T2n.re;
		In.re = I2n.re;
		//虚部
		Tn.im = T2n.im;
		In.im = I2n.im;
		
		for( int i = 0; i < n; i++)
		{
			x = a + ( i + 0.5 ) * D;
			sigma+= f(x).re;
		}
		for( int j = 0; j < n; j++)
		{
			x = a + ( j + 0.5 ) * D;
			siga+= f(x).im;
		}
		
		T2n.re = ( Tn.re + D * sigma )/2.0;
		I2n.re = ( 4 * T2n.re - Tn.re) / 3.0;
		
		T2n.im = ( Tn.im + D * siga )/2.0;
		I2n.im = ( 4 * T2n.im - Tn.im) / 3.0;
		
		n*=2;
		D/=2;
	}
	while(fabs(I2n.re - In.re) >= eps)
	{
		//实部
		Tn.re = T2n.re;
		In.re = I2n.re;
		
		sigma = 0.0;
		
		for( int i = 0; i < n; i++)
		{
			double x = a + ( i + 0.5 ) * D;
			sigma+= f(x).re;
		}
		
		T2n.re = ( Tn.re + D * sigma )/2.0;
		I2n.re = ( 4 * T2n.re - Tn.re) / 3.0;
		
		n*=2;
		D/=2;
	}
	
	while( fabs(I2n.im - In.im) >= eps)
	{	//虚部
		Tn.im = T2n.im;
		In.im = I2n.im;
		
		siga = 0.0;
		
		for( int j = 0; j < n; j++)
		{
			double x = a + ( j + 0.5 ) * D;
			siga+= f(x).im;
		}
		
		T2n.im = ( Tn.im + D * siga )/2.0;
		I2n.im = ( 4 * T2n.im - Tn.im) / 3.0;
		
		n*=2;
		D/=2;
	}			
	integrat=I2n;
	
	return integrat;
}