integral2.cpp

运用物理光学法和矩量法计算带罩天线的远场辐射。

#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include "complex.h"
#include <stdlib.h>
#include "assert.h"
#include "complex.cpp"
#include "parameter.h"

const long double E0=(1./(36.*PI))*1.e-9;
const long double U0=4.*PI*1.e-7;
int m;

complex integrationgn(double Zp,double Zq,double Pp,double Pq,double h,double sita,double Er,int n)//积分方程表达式,Vp是cosVp
{
	double R;
    double K;
    complex F;
    complex gn;
 
    K=k0*sqrt(Er);
	if(Zp==Zq&&Pp==Pq)
		R=sqrt((h/8.)*(h/8.)+2.*Pp*Pq*(1.-cos(sita)));          //h为延轴线的分段（未定义）Vp为母线的切线方向与z轴的夹角
	else
        R=sqrt((Zp-Zq)*(Zp-Zq)+(Pp-Pq)*(Pp-Pq)+2.*Pp*Pq*(1.-cos(sita)));
    F=exp((-1.)*_JC*K*R)/R;
    gn=cos(n*sita)*F;
  return(gn);
}

complex integrationgcn(double Zp,double Zq,double Pp,double Pq,double h,double sita,double Er,int n)//积分方程表达式
{
	double R;
    double K;
    complex F;
    complex gcn;

    K=k0*sqrt(Er);
	if(Zp==Zq&&Pp==Pq)
		R=sqrt((h/8.)*(h/8.)+2.*Pp*Pp*(1.-cos(sita)));//h为延轴线的分段（未定义）Vp为母线的切线方向与z轴的夹角
	else
        R=sqrt((Zp-Zq)*(Zp-Zq)+(Pp-Pq)*(Pp-Pq)+2.*Pp*Pq*(1.-cos(sita)));
    F=exp(_JC*(-1.)*K*R)/R;
    gcn=cos(sita)*cos(n*sita)*F;
    return (gcn);
} 

complex integrationgsn(double Zp,double Zq,double Pp,double Pq,double h,double sita,double Er,int n)//积分方程表达式
{
	double R;
    double K;
    complex F;
    complex gsn;
    K=k0*sqrt(Er);
	if(Zp==Zq&&Pp==Pq)
		R=sqrt((h/8.)*(h/8.)+2.*Pp*Pp*(1.-cos(sita)));//h为延轴线的分段（未定义）Vp为母线的切线方向与z轴的夹角
	else
        R=sqrt((Zp-Zq)*(Zp-Zq)+(Pp-Pq)*(Pp-Pq)+2.*Pp*Pq*(1.-cos(sita)));
    F=exp(_JC*(-1.)*K*R)/R;
    gsn=cos(sita)*cos(n*sita)*F;
    return(gsn);
}
 
complex integrationhn(double Zp,double Zq,double Pp,double Pq,double h,double sita,double Er,int n)//积分方程表达式
{
	double R;
    double K;
    complex F;
    complex hn;
    K=k0*sqrt(Er);
	if(Zp==Zq&&Pp==Pq)
		R=sqrt((h/8.)*(h/8.)+2.*Pp*Pp*(1.-cos(sita)));//h为延轴线的分段（未定义）Vp为母线的切线方向与z轴的夹角
	else
        R=sqrt((Zp-Zq)*(Zp-Zq)+(Pp-Pq)*(Pp-Pq)+2.*Pp*Pq*(1.-cos(sita)));
    F=(1.+_JC*K*R)*exp(_JC*(-1.)*K*R)/pow(R,3);
    hn=cos(n*sita)*F;
    return(hn);
} 

complex integrationhcn(double Zp,double Zq,double Pp,double Pq,double h,double sita,double Er,int n)//积分方程表达式
{
	double R;
    double K;
    complex F;
    complex hcn;
    K=k0*sqrt(Er);
	if(Zp==Zq&&Pp==Pq)
		R=sqrt((h/8.)*(h/8.)+2.*Pp*Pp*(1.-cos(sita)));//h为延轴线的分段（未定义）Vp为母线的切线方向与z轴的夹角
	else
        R=sqrt((Zp-Zq)*(Zp-Zq)+(Pp-Pq)*(Pp-Pq)+2.*Pp*Pq*(1.-cos(sita)));
    F=(1.+_JC*K*R)*exp(_JC*(-1.)*K*R)/pow(R,3);
    hcn=F*cos(sita)*cos(n*sita);
    return(hcn);
} 

complex integrationhsn(double Zp,double Zq,double Pp,double Pq,double h,double sita,double Er,int n)//积分方程表达式
{
	double R;
    double K;
    complex F;
    complex hsn;
    K=2.*PI*sqrt(Er);
	if(Zp==Zq&&Pp==Pq)
		R=sqrt((h/8.)*(h/8.)+2.*Pp*Pp*(1.-cos(sita)));//h为延轴线的分段（未定义）Vp为母线的切线方向与z轴的夹角
	else
        R=sqrt((Zp-Zq)*(Zp-Zq)+(Pp-Pq)*(Pp-Pq)+2.*Pp*Pq*(1.-cos(sita)));
    F=(1.+_JC*K*R)*exp(_JC*(-1.)*K*R)/pow(R,3);
    hsn=F*sin(sita)*sin(n*sita);
    return(hsn);
}
 
complex gn(double Zp,double Zq,double Pp,double Pq,double h,double Er,int n)
{
	complex x=complex(0.,0.);
	if(n<0)
		m=-n;
	else
	{
		if(n==0)
			m=1;
		else
			m=n;
	}
	for(int k=0;k<m;k++)
	{
   	    for(int i=0;i<4;i++)
		{
		    double t;
	        t=PI/(2.*m)*A[i]+(2*k+1)*PI/(2.*m);
		    x+=integrationgn(Zp,Zq,Pp,Pq,h,t,Er,n)*B[i];
		}
	}
	x=x*PI/(2.*m);
	return(x);
}

complex gcn(double Zp,double Zq,double Pp,double Pq,double h,double Er,int n)
{
	
    complex x=complex(0.,0.);
	if(n<0)
		m=-n;
	else
	{
		if(n==0)
			m=1;
		else
			m=n;
	}
	for(int k=0;k<m;k++)
	{
	    for(int i=0;i<4;i++)
		{
		    double t;
		    t=PI/(2.*m)*A[i]+(2*k+1)*PI/(2.*m);
		    x+=integrationgcn(Zp,Zq,Pp,Pq,h,t,Er,n)*B[i];
		}
	}
	x=x*PI/(2.*m);
	return(x);
}
complex gsn(double Zp,double Zq,double Pp,double Pq,double h,double Er,int n)
{
    complex x=complex(0.,0.);
	if(n<0)
		m=-n;
	else
	{
		if(n==0)
			m=1;
		else
			m=n;
	}
	for(int k=0;k<m;k++)
	{
	    for(int i=0;i<4;i++)
		{
		    double t;
		    t=PI/(2.*m)*A[i]+(2*k+1)*PI/(2.*m);
		    x+=integrationgsn(Zp,Zq,Pp,Pq,h,t,Er,n)*B[i];
		}
	}
	x=x*PI/(2.*m);
	return(x);	
}
complex hn(double Zp,double Zq,double Pp,double Pq,double h,double Er,int n)
{
    complex x=complex(0.,0.);
	if(n<0)
		m=-n;
	else
	{
		if(n==0)
			m=1;
		else
			m=n;
	}
	for(int k=0;k<m;k++)
	{
	    for(int i=0;i<4;i++)
		{
		    double t;
		    t=PI/(2.*m)*A[i]+(2*k+1)*PI/(2.*m);
		    x+=integrationhn(Zp,Zq,Pp,Pq,h,t,Er,n)*B[i];
		}
	}
	x=x*PI/(2.*m);
	return(x);
}
complex hcn(double Zp,double Zq,double Pp,double Pq,double h,double Er,int n)
{
    complex x=complex(0.,0.);
	if(n<0)
		m=-n;
	else
	{
		if(n==0)
			m=1;
		else
			m=n;
	}
	for(int k=0;k<m;k++)
	{
	    for(int i=0;i<4;i++)
		{
		    double t;
		    t=PI/(2.*m)*A[i]+(2*k+1)*PI/(2.*m);
		    x+=integrationhcn(Zp,Zq,Pp,Pq,h,t,Er,n)*B[i];
		}
	}
	x=x*PI/(2.*m);
	return(x);
}
complex hsn(double Zp,double Zq,double Pp,double Pq,double h,double Er,int n)
{
    complex x=complex(0.,0.);
	if(n<0)
		m=-n;
	else
	{
		if(n==0)
			m=1;
		else
			m=n;
	}
	for(int k=0;k<m;k++)
	{
	    for(int i=0;i<4;i++)
		{
		    double t;
		    t=PI/(2.*m)*A[i]+(2*k+1)*PI/(2.*m);
		    x+=integrationhsn(Zp,Zq,Pp,Pq,h,t,Er,n)*B[i];
		}
	}
	x=x*PI/(2.*m);
	return(x);
}
double Z1p[4];
double X1p[4];
double zh1p[4];
double yu1p[4];
double h1p[4];

double Z1q[4];
double X1q[4];
double zh1q[4];
double yu1q[4];
double h1q[4];

complex Ltt(double Z1,double Z2,double l1,double l2,double R1,double R2,double Er,int n)
{
	complex L=complex(0.,0.);
	for(int p=0;p<4;p++)
	{
		 Z1p[p]=Z1-(2*p+1)*l1/4.;                       //内母线
		 X1p[p]=sqrt(Rn*Rn-Z1p[p]*Z1p[p])-O;
		 zh1p[p]=(-1.)*Z1p[p]/R1;                       //正弦值
		 yu1p[p]=sqrt(R1*R1-Z1p[p]*Z1p[p])/R1;          //余弦值
		 h1p[p]=l1/yu1p[p];
		 Z1q[p]=Z2-(2*p+1)*l2/4.;
		 X1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])-O;
		 zh1q[p]=(-1.)*Z1q[p]/R2;
		 yu1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])/R2;
		 h1q[p]=l2/yu1q[p];
	}
		for(p=0;p<4;p++)
		{
			for(int q=0;q<4;q++)
			{
				L+=_JC*omiga*U0*T1[p]*T1[q]*(zh1p[p]*zh1q[q]*(h1p[p]/2.)*(h1q[q]/2.)*gcn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n)
						+yu1p[p]*yu1q[q]*(h1p[p]/2.)*(h1q[q]/2.)*gn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n))
						-(_JC/omiga/E0/Er)*T2[p]*T2[q]*(h1p[p]/2.)*(h1q[q]/2.)*gn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n);
			}

		}
		return(L);
}

complex Ltf(double Z1,double Z2,double l1,double l2,double R1,double R2,double Er,int n)
{
	complex L=complex(0.,0.);
	for(int p=0;p<4;p++)
	{
		 Z1p[p]=Z1-(2*p+1)*l1/4.;                       //内母线
		 X1p[p]=sqrt(Rn*Rn-Z1p[p]*Z1p[p])-O;
		 zh1p[p]=(-1.)*Z1p[p]/R1;                       //正弦值
		 yu1p[p]=sqrt(R1*R1-Z1p[p]*Z1p[p])/R1;          //余弦值
		 h1p[p]=l1/yu1p[p];
		 Z1q[p]=Z2-(2*p+1)*l2/4.;
		 X1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])-O;
		 zh1q[p]=(-1.)*Z1q[p]/R2;
		 yu1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])/R2;
		 h1q[p]=l2/yu1q[p];
	}
	for(p=0;p<4;p++)
	{
		for(int q=0;q<4;q++)
		{	
			L+=omiga*U0*zh1p[p]*T1[p]*T1[q]*(h1p[p]/2.)*(h1q[q]/2.)*gsn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n)
				+n/omiga/Er/E0/X1q[q]*T2[p]*T1[q]*(h1p[p]/2.)*(h1q[q]/2.)*gn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n);
		}
	}
	return(L);
}

complex Lft(double Z1,double Z2,double l1,double l2,double R1,double R2,double Er,int n)
{
	complex L=complex(0.,0.);
	for(int p=0;p<4;p++)
	{
	     Z1p[p]=Z1-(2*p+1)*l1/4.;                       //内母线
		 X1p[p]=sqrt(Rn*Rn-Z1p[p]*Z1p[p])-O;
		 zh1p[p]=(-1.)*Z1p[p]/R1;                       //正弦值
		 yu1p[p]=sqrt(R1*R1-Z1p[p]*Z1p[p])/R1;          //余弦值
		 h1p[p]=l1/yu1p[p];
		 Z1q[p]=Z2-(2*p+1)*l2/4.;
		 X1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])-O;
		 zh1q[p]=(-1.)*Z1q[p]/R2;
		 yu1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])/R2;
		 h1q[p]=l2/yu1q[p];
	}
	for(p=0;p<4;p++)
	{
		for(int q=0;q<4;q++)
		{	
			L+=(-1.)*(omiga*U0*zh1q[q]*T1[p]*T1[q]*(h1p[p]/2.)*(h1q[q]/2.)*gsn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n)
				+n/omiga/Er/E0/X1p[p]*T1[p]*T2[q]*(h1p[p]/2.)*(h1q[q]/2.)*gn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n));
		}
	}
	return(L);
}
complex Lff(double Z1,double Z2,double l1,double l2,double R1,double R2,double Er,int n)
{
	complex L=complex(0.,0.);
	for(int p=0;p<4;p++)
	{
		 Z1p[p]=Z1-(2*p+1)*l1/4.;                       //内母线
		 X1p[p]=sqrt(Rn*Rn-Z1p[p]*Z1p[p])-O;
		 zh1p[p]=(-1.)*Z1p[p]/R1;                       //正弦值
		 yu1p[p]=sqrt(R1*R1-Z1p[p]*Z1p[p])/R1;          //余弦值
		 h1p[p]=l1/yu1p[p];
		 Z1q[p]=Z2-(2*p+1)*l2/4.;
		 X1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])-O;
		 zh1q[p]=(-1.)*Z1q[p]/R2;
		 yu1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])/R2;
		 h1q[p]=l2/yu1q[p];
	}
	for(p=0;p<4;p++)
	{
		for(int q=0;q<4;q++)
		{	
			L+=_JC*omiga*U0*T1[p]*T1[q]*(h1p[p]/2.)*(h1q[q]/2.)*gcn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n)
				+n*n*T1[p]*T1[q]/_JC/omiga/Er/E0/X1p[p]/X1q[q]*(h1p[p]/2.)*(h1q[q]/2.)*gn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n);
		}
	}
	return(L);
}

complex Ktt(double Z1,double Z2,double l1,double l2,double R1,double R2,double Er,int n)
{
	complex L=complex(0.,0.);
	for(int p=0;p<4;p++)
	{
		 Z1p[p]=Z1-(2*p+1)*l1/4.;                       //内母线
		 X1p[p]=sqrt(Rn*Rn-Z1p[p]*Z1p[p])-O;
		 zh1p[p]=(-1.)*Z1p[p]/R1;                       //正弦值
		 yu1p[p]=sqrt(R1*R1-Z1p[p]*Z1p[p])/R1;          //余弦值
		 h1p[p]=l1/yu1p[p];
		 Z1q[p]=Z2-(2*p+1)*l2/4.;
		 X1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])-O;
		 zh1q[p]=(-1.)*Z1q[p]/R2;
		 yu1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])/R2;
		 h1q[p]=l2/yu1q[p];
	}
	for(p=0;p<4;p++)
	{
		for(int q=0;q<4;q++)
		{	
			L-=_JC*T1[p]*T1[q]*((Z1p[p]-Z1q[q])*zh1p[p]*zh1q[q]-X1p[p]*zh1q[q]*yu1p[p]+X1q[q]*zh1p[p]*yu1q[q])
				*(h1p[p]/2.)*(h1q[q]/2.)*hsn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n);
		}
	}
	return(L);
}

complex Ktf(double Z1,double Z2,double l1,double l2,double R1,double R2,double Er,int n)
{
	complex L=complex(0.,0.);
	for(int p=0;p<4;p++)
	{
		 Z1p[p]=Z1-(2*p+1)*l1/4.;                       //内母线
		 X1p[p]=sqrt(Rn*Rn-Z1p[p]*Z1p[p])-O;
		 zh1p[p]=(-1.)*Z1p[p]/R1;                       //正弦值
		 yu1p[p]=sqrt(R1*R1-Z1p[p]*Z1p[p])/R1;          //余弦值
		 h1p[p]=l1/yu1p[p];
		 Z1q[p]=Z2-(2*p+1)*l2/4.;
		 X1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])-O;
		 zh1q[p]=(-1.)*Z1q[p]/R2;
		 yu1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])/R2;
		 h1q[p]=l2/yu1q[p];
	}
	for(p=0;p<4;p++)
	{
		for(int q=0;q<4;q++)
		{	
			L-=T1[p]*T1[q]*(zh1p[p]*(Z1p[p]-Z1q[q])*(h1p[p]/2.)*(h1q[q]/2.)*hcn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n)
				+yu1p[p]*(X1q[q]*(h1p[p]/2.)*(h1q[q]/2.)*hn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n)
				-X1p[p]*(h1p[p]/2.)*(h1q[q]/2.)*hcn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n)));
		}
	}
	return(L);
}

complex Kft(double Z1,double Z2,double l1,double l2,double R1,double R2,double Er,int n)
{
	complex L=complex(0.,0.);
	for(int p=0;p<4;p++)
	{
		 Z1p[p]=Z1-(2*p+1)*l1/4.;                       //内母线
		 X1p[p]=sqrt(Rn*Rn-Z1p[p]*Z1p[p])-O;
		 zh1p[p]=(-1.)*Z1p[p]/R1;                       //正弦值
		 yu1p[p]=sqrt(R1*R1-Z1p[p]*Z1p[p])/R1;          //余弦值
		 h1p[p]=l1/yu1p[p];
		 Z1q[p]=Z2-(2*p+1)*l2/4.;
		 X1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])-O;
		 zh1q[p]=(-1.)*Z1q[p]/R2;
		 yu1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])/R2;
		 h1q[p]=l2/yu1q[p];
	}
	for(p=0;p<4;p++)
	{
		for(int q=0;q<4;q++)
		{	
			L-=T1[p]*T1[q]*(zh1q[q]*(Z1q[q]-Z1p[p])*(h1p[p]/2.)*(h1q[q]/2.)*hcn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n)
				+yu1q[q]*(X1p[p]*(h1p[p]/2.)*(h1q[q]/2.)*hn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n)
				-X1q[q]*(h1p[p]/2.)*(h1q[q]/2.)*hcn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n)));
		}
	}
	return(L);
}

complex Kff(double Z1,double Z2,double l1,double l2,double R1,double R2,double Er,int n)
{
	complex L=complex(0.,0.);
	for(int p=0;p<4;p++)
	{
		 Z1p[p]=Z1-(2*p+1)*l1/4.;                       //内母线
		 X1p[p]=sqrt(Rn*Rn-Z1p[p]*Z1p[p])-O;
		 zh1p[p]=(-1.)*Z1p[p]/R1;                       //正弦值
		 yu1p[p]=sqrt(R1*R1-Z1p[p]*Z1p[p])/R1;          //余弦值
		 h1p[p]=l1/yu1p[p];
		 Z1q[p]=Z2-(2*p+1)*l2/4.;
		 X1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])-O;
		 zh1q[p]=(-1.)*Z1q[p]/R2;
		 yu1q[p]=sqrt(R2*R2-Z1q[p]*Z1q[p])/R2;
		 h1q[p]=l2/yu1q[p];
	}
	for(p=0;p<4;p++)
	{
		for(int q=0;q<4;q++)
		{	
			L-=_JC*T1[p]*T1[q]*(Z1p[p]-Z1q[q])*(h1p[p]/2.)*(h1q[q]/2.)*hsn(Z1p[p],Z1q[q],X1p[p],X1q[q],h1p[p],Er,n);
		}
	}
	return(L);
}