banbozhenzi.cpp

微代天线计算中的半波阵子的奇数阶自阻抗计算。

#include <complex>

#include <iostream>
#include <fstream>
#include <math.h>
#include "gauss.h"

#define Z0 50.0
#define Y0 0.02

using namespace std;
#define PI 3.14159265358979323846264
#define GAMMA 0.57721566490153286060651209
#define NN 5  //阵元数为2NN+1
#define MM 21  //无穷大的取值
complex<double>  ZZ[2*MM+1][2*MM+1]={0.0};
#include "inverse.h"
complex<double> i(0.0,1.0);
double lmt=0.238,k=2.0*PI/lmt,l=0.055,d=0.119,a=0.002;//k为波数，l为阵子的半长度，d为阵子间的距离，a为阵子的直径
//double theta0=PI/3, drt=k*d*sin(theta0); //theta为扫角

complex<double> func_ci(int q,int p,double x)
{ return((cos(x)-1)/x); 
}

complex<double> func_si(int q,int p,double x)
{ return(sin(x)/x);
}

complex<double> func_jf1(int q,int p,double z)
{ 
  double r1=sqrt((z-l)*(z-l)+double((q-p)*(q-p))*d*d);
  double r2=sqrt((z+l)*(z+l)+double((q-p)*(q-p))*d*d);
  double r0=sqrt(z*z+double((q-p)*(q-p))*d*d);
  complex<double> y1=sin(k*(l+z))*(exp(-i*k*r1)/r1+exp(-i*k*r2)/r2-2*cos(k*l)*exp(-i*k*r0)/r0); 
  return y1;
}

complex<double> func_jf2(int q,int p,double z)
{ 
  double r1=sqrt((z-l)*(z-l)+double((q-p)*(q-p))*d*d);
  double r2=sqrt((z+l)*(z+l)+double((q-p)*(q-p))*d*d);
  double r0=sqrt(z*z+double((q-p)*(q-p))*d*d);
  complex<double> y2=sin(k*(l-z))*(exp(-i*k*r1)/r1+exp(-i*k*r2)/r2-2*cos(k*l)*exp(-i*k*r0)/r0); 
  return y2;
}

complex<double> func_S_inf(int q,int p,double x)
{
  complex<double> temp_complex=0.0;
	for(int pp=-MM;pp<=MM;pp++)
	 { 
		temp_complex=ZZ[q+MM][pp+MM]*exp(i*(complex<double> ((q-pp)*x)))+temp_complex;
	 }
	   temp_complex=(temp_complex-Z0)/(temp_complex+Z0); 
	   temp_complex=temp_complex*exp(-i*(complex<double> ((q-p)*x)));
	 return temp_complex;
}

complex<double> func_Y_inf(int q,int p,double x)
{
  complex<double> temp_complex=0.0;
	for(int pp=-MM;pp<=MM;pp++)
	 { 
		temp_complex=ZZ[q+MM][pp+MM]*exp(i*(complex<double> ((q-pp)*x)))+temp_complex;
	 }
	    temp_complex=Z0*Y0/temp_complex;
        temp_complex=temp_complex*exp(-i*(complex<double> ((q-p)*x)));
	 return temp_complex;
}


void main()
{  
  /*首先求解有关无限阵的阻抗、导纳、散射参量*/	
	ofstream fout_Z("Z.txt"),fout_ZZ("ZZ.txt"),fout_S_inf("S_inf.txt"),fout_Y_inf("Y_inf.txt");   
    ofstream fout_sigma1("sigma1.txt"),fout_sigma2("sigma2.txt"),fout_sigma3("sigma3.txt"),fout_sigma4("sigma4.txt"),fout_s1("s1.txt"),fout_s2("s2.txt"),fout_s_qp("s_qp.txt");
    ofstream fout_z_qp("z_qp.txt");
	
	/*********************以下求自阻抗矩阵**********************/
	complex<double> ci_2kl,ci_4kl,si_2kl,si_4kl,Rll,ci_kal_1,ci_kal,Xll;

       int N_GAUSS=6;		//待求 gauss 阶数
       double *ZERO_GAUSS = new double [N_GAUSS];				//n阶Legendre 多项式的零点
	   double *XI_GAUSS = new double [N_GAUSS];				//求积的系数 Ai
       IntGaussLegendre(N_GAUSS, ZERO_GAUSS, XI_GAUSS);
  
  double u=0.0,v1=2.0*k*l,v2=4.0*k*l,v3=k*a*a/l;
  ci_2kl=GAMMA+log(v1)+GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_ci, u, v1,2,1);
  ci_4kl=GAMMA+log(v2)+GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_ci, u, v2,2,1);
  si_2kl=GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_si, u, v1,2,1);
  si_4kl=GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_si, u, v2,2,1);
  Rll=60.0*(GAMMA+log(2.0*k*l)-ci_2kl+0.5*sin(2.0*k*l)*(si_4kl-2.0*si_2kl)+0.5*cos(2.0*k*l)*(GAMMA+log(k*l)+ci_4kl-2.0*ci_2kl));
  //ci_kal_1=GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_ci, u, v3,2,1);
  //ci_kal=GAMMA+log(v1)+ci_kal_1;
  Xll=30.0*(2.0*si_2kl+sin(2.0*k*l)*(GAMMA+log(k*l)+ci_4kl-2.0*ci_2kl-2.0*log(l/a))+cos(2.0*k*l)*(2.0*si_2kl-si_4kl));
  //cout<<"Rll="<<Rll<<"\n";
  //cout<<"Xll="<<Xll<<"\n";
  
	   /*******************自阻抗求解结束**********************/

  
  






  /***************************以下为求解互阻抗*************/ 
  complex<double>  R[2*NN+1][2*NN+1]={0.0},RR[2*MM+1][2*MM+1]={0.0},X[2*NN+1][2*NN+1]={0.0},XX[2*MM+1][2*MM+1]={0.0},I[2*NN+1][2*NN+1]={0.0},II[2*MM+1][2*MM+1]={0.0},Z[2*NN+1][2*NN+1]={0.0};
  //complex<double>  S_1[2*NN+1][2*NN+1]={0.0},S_2[2*NN+1][2*NN+1]={0.0};
  int q,p;
  
  for(q=-NN;q<=NN;q++)
  {
   for(p=-NN;p<=NN;p++)
   { 
	   if(q==p) 
	   {
             R[q+NN][p+NN]=Rll;
             X[q+NN][p+NN]=Xll;
			 Z[q+NN][p+NN]=R[q+NN][p+NN]+i*X[q+NN][p+NN];  //自阻抗
			 I[q+NN][p+NN]=1.0;
	   }
       else
	   { 
		     Z[q+NN][p+NN]=i*30.0*(GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_jf1, -l, 0.0, q, p)+GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_jf2, 0.0, l, q, p));
                  //互阻抗
             I[q+NN][p+NN]=0.0;
       }
	 
	 //cout<<"Z="<<Z[q+NN][p+NN]<<"\n";
	 fout_Z<<Z[q+NN][p+NN]<<"\t";
	   
   }
     fout_Z<<"\n";
  }
  /*******************互阻抗求解结束**********************/ 

   for(q=-MM;q<=MM;q++)
   {
     for(p=-MM;p<=MM;p++)
	 { 
		 if(q==p) 
	   {
             RR[q+MM][p+MM]=Rll;
             XX[q+MM][p+MM]=Xll;
			 ZZ[q+MM][p+MM]=RR[q+MM][p+MM]+i*XX[q+MM][p+MM];  //自阻抗
			 II[q+MM][p+MM]=1.0;
	   }
       else
	   { 
		     ZZ[q+MM][p+MM]=i*30.0*(GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_jf1, -l, 0.0, q, p)+GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_jf2, 0.0, l, q, p));
                  //互阻抗
             II[q+MM][p+MM]=0.0;
       }
	   fout_ZZ<<ZZ[q+MM][p+MM]<<"\t";
	 }
	  fout_ZZ<<"\n";   
   }

//以下无限阵的参量计算
  complex<double> S_inf[2*NN+1][2*NN+1]={0.0}, Y_inf[2*NN+1][2*NN+1]={0.0};
  for(q=-NN;q<=NN;q++)
  {
   for(p=-NN;p<=NN;p++)
   { 
     S_inf[q+NN][p+NN]=(1.0/2.0/PI)*GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_S_inf, -PI, -0.5*PI, q, p)+(1.0/2.0/PI)*GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_S_inf, -0.5*PI, 0.0, q, p)+(1.0/2.0/PI)*GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_S_inf, 0.0, 0.5*PI, q, p)+(1.0/2.0/PI)*GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_S_inf, 0.5*PI, PI, q, p);
	 Y_inf[q+NN][p+NN]=(1.0/2.0/PI)*GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_Y_inf, -PI, -0.5*PI, q, p)+(1.0/2.0/PI)*GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_Y_inf, -0.5*PI, 0.0, q, p)+(1.0/2.0/PI)*GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_Y_inf, 0.0, 0.5*PI, q, p)+(1.0/2.0/PI)*GaussInt(N_GAUSS,XI_GAUSS, ZERO_GAUSS, func_Y_inf, 0.5*PI, PI, q, p);
	 
	 fout_S_inf<<S_inf[q+NN][p+NN]<<"\t";  
	 fout_Y_inf<<Y_inf[q+NN][p+NN]<<"\t"; 
   }
     fout_S_inf<<"\n";
	 fout_Y_inf<<"\n";
  }
   
//以下求解有限阵的散射矩阵和阻抗矩阵

complex<double>  sigma1[2*NN+1][2*NN+1]={0.0},sigma2[2*NN+1][2*NN+1]={0.0},sigma3[2*NN+1][2*NN+1]={0.0},sigma4[2*NN+1][2*NN+1]={0.0};
complex<double>  s1[2*NN+1][2*NN+1]={0.0},s2[2*NN+1][2*NN+1]={0.0},s_qp[2*NN+1][2*NN+1]={0.0};

  for(q=-NN;q<=NN;q++)
  {
     for(p=-NN;p<=NN;p++)
	 { 
         for(int n1=-NN;n1<=NN;n1++)
		 {
			 sigma1[q+NN][p+NN]=(I[n1+NN][p+NN]+S_inf[n1+NN][p+NN])*Y_inf[q+NN][n1+NN]/Y0+sigma1[q+NN][p+NN];
		 }
		 
		 fout_sigma1<<sigma1[q+NN][p+NN]<<"\t";
	 }
	     fout_sigma1<<"\n";//输出sigma1到文本文档	 
  }

  
  for(q=-NN;q<=NN;q++)
  {
     for(p=-NN;p<=NN;p++)
	 { 

			 s1[q+NN][p+NN]=I[q+NN][p+NN]/2.0-S_inf[q+NN][p+NN]/2.0-sigma1[q+NN][p+NN]/2.0;
		 
		 fout_s1<<s1[q+NN][p+NN]<<"\t";
	 }
	     fout_s1<<"\n";//输出s1到文本文档	 
  }


  for(q=-NN;q<=NN;q++)
  {
     for(p=-NN;p<=NN;p++)
	 { 
         for(int n=-NN;n<=NN;n++)
		 {
			 sigma2[q+NN][p+NN]=(I[n+NN][p+NN]+S_inf[n+NN][p+NN])*S_inf[q+NN][n+NN]+sigma2[q+NN][p+NN];
		 }
		 
		 fout_sigma2<<sigma2[q+NN][p+NN]<<"\t";
	 }
	     fout_sigma2<<"\n";//输出sigma2到文本文档	 
  }


  for(int n=-NN;n<=NN;n++)
  {
     for(p=-NN;p<=NN;p++)
	 { 
         for(int n1=-NN;n1<=NN;n1++)
		 {
			 sigma3[n+NN][p+NN]=(I[n1+NN][p+NN]+S_inf[n1+NN][p+NN])*Y_inf[n+NN][n1+NN]/Y0+sigma3[n+NN][p+NN];
		 }
		 
		 fout_sigma3<<sigma3[n+NN][p+NN]<<"\t";
	 }
	     fout_sigma3<<"\n";//输出sigma3到文本文档	 
  }


  for(q=-NN;q<=NN;q++)
  {
     for(p=-NN;p<=NN;p++)
	 { 
         for(int n=-NN;n<=NN;n++)
		 {
			 sigma4[q+NN][p+NN]=S_inf[q+NN][n+NN]*sigma3[n+NN][p+NN]+sigma4[q+NN][p+NN];
		 }
		 
		 fout_sigma4<<sigma4[q+NN][p+NN]<<"\t";
	 }
	     fout_sigma4<<"\n";//输出sigma4到文本文档	 
  }


  for(q=-NN;q<=NN;q++)
  {
     for(p=-NN;p<=NN;p++)
	 { 
		s2[q+NN][p+NN]=S_inf[q+NN][p+NN]-sigma2[q+NN][p+NN]/2.0-sigma4[q+NN][p+NN]/2.0;
		s_qp[q+NN][p+NN]=S_inf[q+NN][p+NN]+s1[q+NN][p+NN]+s2[q+NN][p+NN];
		 fout_s2<<s2[q+NN][p+NN]<<"\t";
		 fout_s_qp<<s_qp[q+NN][p+NN]<<"\t";
	 }
	     fout_s2<<"\n";//输出s2到文本文档	 
	     fout_s_qp<<"\n";//输出s_qp到文本文档	
  }
/*求得有限阵的散射与导纳矩阵*/

complex<double>  z1[2*NN+1][2*NN+1]={0.0},z2[2*NN+1][2*NN+1]={0.0},z_1[2*NN+1][2*NN+1]={0.0},z_qp[2*NN+1][2*NN+1]={0.0};

  for(q=-NN;q<=NN;q++)
  {
     for(p=-NN;p<=NN;p++)
	 { 
        z1[q+NN][p+NN]=I[q+NN][p+NN]+s_qp[q+NN][p+NN];
		z2[q+NN][p+NN]=I[q+NN][p+NN]-s_qp[q+NN][p+NN];
		cout<<"z1="<<z1[q+NN][p+NN]<<"\t";
        cout<<"z2="<<z2[q+NN][p+NN]<<"\n";
	 }	 
  }

  Inverse(z2);//求z2的逆，并将其求逆的值赋值于z2


  for(q=-NN;q<=NN;q++)
  {
     for(p=-NN;p<=NN;p++)
	 { 
		for(int m=-NN;m<=NN;m++)
		 {
			 z_1[q+NN][p+NN]=z1[q+NN][m+NN]*z2[m+NN][p+NN]+z_1[q+NN][p+NN];//得到阻抗矩阵
             z_qp[q+NN][p+NN]=z_1[q+NN][p+NN]*Z0;
		 }
	    cout<<"z2="<<z2[q+NN][p+NN]<<"\n";
	    fout_z_qp<<z_qp[q+NN][p+NN]<<"\t";
	 }	 
	    fout_z_qp<<"\n";//输出z_qp到文本文档	
  }

  complex<double> Zin[2*NN+1]={0.0};//输入阻抗求解
  ofstream fout_Zin("Zin.txt");
  
    for(q=-NN;q<=NN;q++)
	{	
	  for(int p=-NN;p<=NN;p++)
	  { 
		Zin[q+NN]=z_qp[q+NN][p+NN]+Zin[q+NN];
	  }
    fout_Zin<<Zin[q+NN]<<"\t";
	cout<<"Zin="<<Zin[q+NN]<<"\n";
	}
  fout_Zin<<"\n";


 complex<double> Zin_inf[2*MM+1]={0.0};//无限阵的输入阻抗
 ofstream fout_Zin_inf("Zin_inf.txt");
 for(q=-MM;q<=MM;q++)
   {
     for(p=-MM;p<=MM;p++)
	 { 
	Zin_inf[q+MM]=ZZ[q+MM][p+MM]+Zin_inf[q+MM];	
       }
	   fout_Zin_inf<<Zin_inf[q+MM]<<"\t";
	 }
	  fout_Zin_inf<<"\n";   

 double Z_qp_abs[2*NN+1][2*NN+1]={0.0};//有限阵阻抗的模值
 ofstream fout_Z_qp_abs("Z_qp_abs.txt");
 for(q=-NN;q<=NN;q++)
   {
     for(p=-NN;p<=NN;p++)
	 { 
	   Z_qp_abs[q+NN][p+NN]=abs(z_qp[q+NN][p+NN]);	
	   fout_Z_qp_abs<<Z_qp_abs[q+NN][p+NN]<<"\t";
	  }
       fout_Z_qp_abs<<"\n";
 }

double Z_qp_abs1[2*NN+1]={0.0};
double Z_qp_abs2[2*NN+1]={0.0}; //取出有限阵的一列
ofstream fout_Z_qp_abs1("Z_qp_abs1.txt");
ofstream fout_Z_qp_abs2("Z_qp_abs2.txt");
 for(q=-NN;q<=NN;q++)
   {
	   Z_qp_abs1[q+NN]=Z_qp_abs[q+NN][0];
	   Z_qp_abs2[q+NN]=Z_qp_abs[q+NN][5];
	   fout_Z_qp_abs1<<Z_qp_abs1[q+NN]<<"\t";
	   fout_Z_qp_abs2<<Z_qp_abs2[q+NN]<<"\t";
	  }
       fout_Z_qp_abs1<<"\n";
       fout_Z_qp_abs2<<"\n";


	     

double Z_abs[2*NN+1][2*NN+1]={0.0};//无限阵阻抗元素的模值
 ofstream fout_Z_abs("Z_abs.txt");
 for(q=-NN;q<=NN;q++)
   {
     for(p=-NN;p<=NN;p++)
	 { 
	   Z_abs[q+NN][p+NN]=abs(Z[q+NN][p+NN]);	
	   fout_Z_abs<<Z_abs[q+NN][p+NN]<<"\t";
	  }
	   fout_Z_abs<<"\n";
	 }
	  
double Z_abs1[2*NN+1]={0.0};//取出其中两列
double Z_abs2[2*NN+1]={0.0};
ofstream fout_Z_abs1("Z_abs1.txt");
ofstream fout_Z_abs2("Z_abs2.txt");
 for(q=-NN;q<=NN;q++)
   {
	   Z_abs1[q+NN]=Z_abs[q+NN][0];
	   Z_abs2[q+NN]=Z_abs[q+NN][5];
	   fout_Z_abs1<<Z_abs1[q+NN]<<"\t";
	   fout_Z_abs2<<Z_abs2[q+NN]<<"\t";
	  }
       fout_Z_abs1<<"\n";
       fout_Z_abs2<<"\n";




}