VC++：无限细半波振子的辐射阻抗（感应电动势法）
//无限细半波振子间的互阻抗（经典结果）#include <iostream.h>
#include <iomanip.h>
#include <math.h>
#define EPS 1E-7 //正(余)弦积分的误差界
const double L = 0.5; //半波振子全长
const double M_PI =3.14159265358979324;
const double TPI =2 * M_PI;
const double Euler=0.57721566490153286; //欧拉常数
double Si(double x) //这是所谓正弦积分
{
double Si;
if(fabs(x)<=15)
{ //以下为泰勒展开
double ai=3,x2=x*x,sgn,t;
sgn=(x<0?1:-1);
t=x2*x/18;
Si=x-t;
while(fabs(t)>EPS)
{
ai+=2;
t*=x2/((ai-1)*ai*ai)*(ai-2);
sgn=-sgn;
Si+=sgn*t;
}
}
else
{ //以下为渐近展开
double x1=1/x,x2=x*x,c,s,ai=0,sgn=1,sii,sij,t;
c=cos(x)*x1;
s=sin(x)*x1;
sii=t=c;
while(fabs(t)>EPS)
{
ai+=2;
t*=ai*(ai-1)/x2;
sgn=-sgn;
sii+=sgn*t;
}
ai=sgn=1;
if(x<0)sgn=-1;
sij=t=s*x1;
while(fabs(t)>EPS)
{
ai+=2;
t*=ai*(ai-1)/x2;
sgn=-sgn;
sij+=sgn*t;
}
Si=M_PI/2-sii-sij;
}
return Si;
}
double Ci(double x) //这是所谓余弦积分
{
double Ci;
if(fabs(x)<=15)
{
double ai=2,x2,sgn=-1,t;
x2=x*x;
t=x2/4;
Ci=Euler+log(fabs(x))-t;
while(fabs(t)>EPS)
{
ai+=2;
t*=x2/((ai-1)*ai*ai)*(ai-2);
sgn=-sgn;
Ci+=sgn*t;
}
}
else
{
double x1=1/x,x2=x*x,c,s,ai=0,sgn=1,cii,cij,t;
c=cos(x)*x1;
s=sin(x)*x1;
cii=t=s;
while(fabs(t)>EPS)
{
ai+=2;
t*=ai*(ai-1)/x2;
sgn=-sgn;
cii+=sgn*t;
}
ai=sgn=1;
if(x<0)sgn=-1;
cij=t=c*x1;
while(fabs(t)>EPS)
{
ai+=2;
t*=ai*(ai-1)/x2;
sgn=-sgn;
cij+=sgn*t;
}
Ci=cii-cij;
}
return Ci;
}
double f(double x,double y)
{
return TPI*(y+sqrt(x*x+y*y));
}
int main(void)
{
double R12,X12,w0,m0,w1,m1,w2,m2,w_1,m_1,w_2,m_2;
double d,H,Bd,BH,BL=TPI*L;
begin:
cout<< "d=? H=?" <<endl;
cin >> d >> H ;
if(d<0)d=-d;
if(H<0)H=-H;
//以上，d和H分别是两振子间的平行距离和纵向距离(纵向错开)，以波长为单位。
//L为振子全长，对于半波振子L=0.5（波长）。无限细半波振子的自辐射阻抗为
//73.1+j42.5Ω，就是在参数 d 和 H 均趋于零情况下的极限。
if(!d && H && H<L)
{
cout << "积分发散" << endl;
return -1;
}
Bd=TPI*d;
BH=TPI*H;
w0=f(d, (H ));
w1=f(d, (H+L/2));
w2=f(d, (H+L ));
w_1=f(d,(H-L/2));
w_2=f(d,(H-L ));
m0=f(d,-(H ));
m1=f(d,-(H+L/2));
m2=f(d,-(H+L ));
m_1=f(d,-(H-L/2));
m_2=f(d,-(H-L ));
if(d==0 && H==0 && L==0.5)//求半波振子的自辐射阻抗(欧)
{
R12=30*(-Ci(TPI)+Euler+log(TPI));
X12=30*(Si(TPI));
}
else if(d==0 && H==L && L==0.5)//准全波振子
{
R12= 15*cos( BH )*(3*Ci(w0)-2*Ci(w1)+2*log(fabs(1-L*L/(4*H*H)))-(log(BL/2)+Euler))
+15*cos(BH+BL)*(Ci(w2)-2*Ci(w1)+Ci(w0)+log(fabs(pow(1+L/H/2,2)/(1+L/H))));
X12=-30*cos( BH )*(-Si(w_1)+2*Si(w0)-Si(w1))-15*cos(BH-BL)*(Si(w_2)-2*Si(w_1)+Si(w0))
-15*cos(BH+BL)*(Si(w2)-2*Si(w1)+Si(w0));
}
else if(d==0 && H>L && L==0.5)//共线半波振子
{
R12= 15*cos( BH )*(3*Ci(w0)-2*Ci(w1)+2*log(fabs(1-pow(L/H/2,2)))
-Ci(w_2)-log(fabs(pow(1-L/H/2,2)/(1-L/H))))
+30*sin( BH )*(-Si(w_1)+2*Si(w0)-Si(w1))
+15*sin(BH-BL)*(Si(w_2)-2*Si(w_1)+Si(w0))
+15*cos(BH+BL)*(Ci(w2 )-2*Ci(w1 )+Ci(w0) +log(fabs(pow(1+L/H/2,2)/(1+L/H))))
+15*sin(BH+BL)*(Si(w2 )-2*Si(w1 )+Si(w0));
X12=-30*cos( BH )*(-Si(w_1)+2*Si(w0)-Si(w1))
+30*sin( BH )*(-Ci(w_1)+2*Ci(w0)-Ci(w1) -log(fabs(1-pow(L/H/2,2))))
-15*cos(BH-BL)*(Si(w_2)-2*Si(w_1)+Si(w0))
+15*sin(BH-BL)*(Ci(w_2)-2*Ci(w_1)+Ci(w0) -log(fabs(pow(1-L/H/2,2)/(1-L/H))))
-15*cos(BH+BL)*(Si(w2 )-2*Si(w1 )+Si(w0))
+15*sin(BH+BL)*(Ci(w2 )-2*Ci(w1 )+Ci(w0) -log(fabs(pow(1+L/H/2,2)/(1+L/H))));
}
else if(d==0 && H>L)//共线对称振子
{
R12= 30*cos( BH )*(-Ci(w_1)+2*Ci(w0)-Ci(w1) +log(fabs(1-pow(L/H/2,2))))
+30*sin( BH )*(-Si(w_1)+2*Si(w0)-Si(w1))
+15*cos(BH-BL)*(Ci(w_2)-2*Ci(w_1)+Ci(w0) +log(fabs(pow(1-L/H/2,2)/(1-L/H))))
+15*sin(BH-BL)*(Si(w_2)-2*Si(w_1)+Si(w0))
+15*cos(BH+BL)*(Ci(w2 )-2*Ci(w1 )+Ci(w0) +log(fabs(pow(1+L/H/2,2)/(1+L/H))))
+15*sin(BH+BL)*(Si(w2 )-2*Si(w1 )+Si(w0));
X12=-30*cos( BH )*(-Si(w_1)+2*Si(w0)-Si(w1))
+30*sin( BH )*(-Ci(w_1)+2*Ci(w0)-Ci(w1) -log(fabs(1-pow(L/H/2,2))))
-15*cos(BH-BL)*(Si(w_2)-2*Si(w_1)+Si(w0))
+15*sin(BH-BL)*(Ci(w_2)-2*Ci(w_1)+Ci(w0) -log(fabs(pow(1-L/H/2,2)/(1-L/H))))
-15*cos(BH+BL)*(Si(w2 )-2*Si(w1 )+Si(w0))
+15*sin(BH+BL)*(Ci(w2 )-2*Ci(w1 )+Ci(w0) -log(fabs(pow(1+L/H/2,2)/(1+L/H))));
}
else
{
R12= 30*cos( BH )*(-Ci(w_1)-Ci(m_1)+2*Ci(w0)+2*Ci(m0)-Ci(w1)-Ci(m1))
+30*sin( BH )*(-Si(w_1)+Si(m_1)+2*Si(w0)-2*Si(m0)-Si(w1)+Si(m1))
+15*cos(BH-BL)*(Ci(w_2)+Ci(m_2)-2*Ci(w_1)-2*Ci(m_1)+Ci(w0)+Ci(m0))
+15*sin(BH-BL)*(Si(w_2)-Si(m_2)-2*Si(w_1)+2*Si(m_1)+Si(w0)-Si(m0))
+15*cos(BH+BL)*(Ci(w2)+Ci(m2)-2*Ci(w1)-2*Ci(m1)+Ci(w0)+Ci(m0))
+15*sin(BH+BL)*(Si(w2)-Si(m2)-2*Si(w1)+2*Si(m1)+Si(w0)-Si(m0));
X12=-30*cos( BH )*(-Si(w_1)-Si(m_1)+2*Si(w0)+2*Si(m0)-Si(w1)-Si(m1))
+30*sin( BH )*(-Ci(w_1)+Ci(m_1)+2*Ci(w0)-2*Ci(m0)-Ci(w1)+Ci(m1))
-15*cos(BH-BL)*(Si(w_2)+Si(m_2)-2*Si(w_1)-2*Si(m_1)+Si(w0)+Si(m0))
+15*sin(BH-BL)*(Ci(w_2)-Ci(m_2)-2*Ci(w_1)+2*Ci(m_1)+Ci(w0)-Ci(m0))
-15*cos(BH+BL)*(Si(w2)+Si(m2)-2*Si(w1)-2*Si(m1)+Si(w0)+Si(m0))
+15*sin(BH+BL)*(Ci(w2)-Ci(m2)-2*Ci(w1)+2*Ci(m1)+Ci(w0)-Ci(m0));
}
if(X12<0)
{
cout << "Z12 = ("<< R12 << "-j" << -X12 << ")欧" << endl;
}
else
{
cout << "Z12 = ("<< R12 << "+j" << +X12 << ")欧" << endl;
}
goto begin;
return 0;
}
[ 本帖最后由 yu_hua 于 2010-10-24 15:03 编辑 ]
