C语言实现半波全宽算法,VC++：无限细半波振子的辐射阻抗（感应电动势法）

已结贴√

问题点数：20 回复次数：2

VC++：无限细半波振子的辐射阻抗(感应电动势法)

//无限细半波振子间的互阻抗(经典结果)

#include

#include

#include

#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=?" <

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

{

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 编辑]