📄 far.c
字号:
for(j=0;j<MinusTrngleIndex[i];j++)
fscanf(InF1,"%d",&MinusTrngleDetect[i][j]);
}
fscanf(InF1, " %c", &c);
if (c!='#')
{
fprintf(stderr, "Check your inner edge table\n");
exit(1);
}
fgets(buff, sizeof(buff)-1, InF1);
InnerEdgeStat = INT_Vector(TotInnerEdgeNum);
for(i=0;i<TotInnerEdgeNum;i++)
fscanf(InF1,"%d",&InnerEdgeStat[i]);
/* read boundary edge table */
fscanf(InF1, " %c", &c);
if (c!='#')
{
fprintf(stderr, "Check your boundary edge table\n");
exit(1);
}
fgets(buff, sizeof(buff)-1, InF1);
BoundEdgeStat = INT_Vector(TotBoundEdgeNum);
for(i=0;i<TotBoundEdgeNum;i++)
fscanf(InF1,"%d",&BoundEdgeStat[i]);
fscanf(InF1, " %c", &c);
if (c!='#')
{
fprintf(stderr, "Check your source information table\n");
exit(1);
}
fgets(buff, sizeof(buff)-1, InF1);
fscanf(InF1," %c", &c); /* skip the type */
fscanf(InF1,"%lf",&OperateFreq); /* obtain the working frequency */
MVector = CMPLX_Vector(TotBoundEdgeNum);
JVector = CMPLX_Vector(TotBoundEdgeNum);
for(j=0;j<HybrdBoundEdgeNum; j++)
if (fscanf(InF2,"%lf %lf",&MVector[j].x,&MVector[j].y)!=2)
{
fprintf(stderr, "Surface field file damaged\n");
exit(1);
}
for(j=0; j<TotBoundEdgeNum-HybrdBoundEdgeNum; j++)
MVector[j+HybrdBoundEdgeNum]=COMplex_Cmplx(0.0, 0.0);
for(j=0;j<TotBoundEdgeNum;j++)
if (fscanf(InF2,"%lf %lf",&JVector[j].x,&JVector[j].y)!=2)
{
fprintf(stderr, "Surface field file damaged\n");
exit(1);
}
}
/********************** ComputeFarEField() *******************************/
void ComputeFarEField(double *ObservPoint, double cth, double cph,
double sth, double sph)
{
int IEdge, SrcEdgeNum, SrcNodeNum, EdgeEnd1, EdgeEnd2, RowNum,
SourceEdge, SourceTrngle, TrngleNodeArray[3] ;
double tmp, Sgn;
complex ETheta, EPhi, ATheta, APhi, FTheta, FPhi,
AVectorX, AVectorY, AVectorZ, AVectorPQ[3],
FVectorX, FVectorY, FVectorZ,
V_x, V_y, V_z;
ETheta = COMplex_Null(); EPhi = COMplex_Null();
AVectorX = COMplex_Null(); FVectorX = COMplex_Null();
AVectorY = COMplex_Null(); FVectorY = COMplex_Null();
AVectorZ = COMplex_Null(); FVectorZ = COMplex_Null();
for(SourceTrngle=0;SourceTrngle<TotTrngleNum;SourceTrngle++)
{
InoPQ = COMplex_Null() ;
IetaPQ = COMplex_Null() ;
IxiPQ = COMplex_Null() ; IzetaPQ = COMplex_Null() ;
for(IEdge=0;IEdge<=2;IEdge++)
TrngleNodeArray[IEdge] = TrngleNode[SourceTrngle][IEdge]-1;
ComputeNonSingul(TrngleNodeArray,cth,cph,sth,sph,RadDisObserv);
V_x = COMplex_Add2(Real_Mul(NodeCord[TrngleNodeArray[0]][0], IxiPQ),
Real_Mul(NodeCord[TrngleNodeArray[1]][0], IetaPQ),
Real_Mul(NodeCord[TrngleNodeArray[2]][0],IzetaPQ));
V_y = COMplex_Add2(Real_Mul(NodeCord[TrngleNodeArray[0]][1], IxiPQ),
Real_Mul(NodeCord[TrngleNodeArray[1]][1], IetaPQ),
Real_Mul(NodeCord[TrngleNodeArray[2]][1],IzetaPQ));
V_z = COMplex_Add2(Real_Mul(NodeCord[TrngleNodeArray[0]][2], IxiPQ),
Real_Mul(NodeCord[TrngleNodeArray[1]][2], IetaPQ),
Real_Mul(NodeCord[TrngleNodeArray[2]][2],IzetaPQ));
for(SourceEdge=0;SourceEdge<=2;SourceEdge++)
{
double SrcTrngleEdgeLen;
SrcEdgeNum = abs(TrngleEdge[SourceTrngle][SourceEdge]);
SrcNodeNum = TrngleNode[SourceTrngle][SourceEdge] - 1;
Sgn = Sign(TrngleEdge[SourceTrngle][SourceEdge]);
EdgeEnd1 = GlobalEdgeEnds[SrcEdgeNum-1][0]-1;
EdgeEnd2 = GlobalEdgeEnds[SrcEdgeNum-1][1]-1;
SrcTrngleEdgeLen = VTXmag(NodeCord[EdgeEnd1],NodeCord[EdgeEnd2]);
for(IEdge=0;IEdge<TotBoundEdgeNum;IEdge++)
{
if (BoundEdgeStat[IEdge]==SrcEdgeNum)
{
RowNum=IEdge;
break;
}
else RowNum=-1;
}
if (RowNum!=-1)
{
AVectorPQ[0] = COMplex_Sub(V_x,
Real_Mul(NodeCord[SrcNodeNum][0],InoPQ));
AVectorPQ[1] = COMplex_Sub(
V_y,Real_Mul(NodeCord[SrcNodeNum][1],InoPQ));
AVectorPQ[2] = COMplex_Sub(V_z,
Real_Mul(NodeCord[SrcNodeNum][2],InoPQ));
CFactor1= Real_Mul((Sgn * SrcTrngleEdgeLen),JVector[RowNum]);
AVectorX = COMplex_Add(AVectorX,
COMplex_Mul( CFactor1,AVectorPQ[0]));
AVectorY = COMplex_Add(AVectorY,
COMplex_Mul( CFactor1,AVectorPQ[1]));
AVectorZ = COMplex_Add(AVectorZ,
COMplex_Mul( CFactor1,AVectorPQ[2]));
CFactor1 = Real_Mul((Sgn * SrcTrngleEdgeLen),MVector[RowNum]);
FVectorX = COMplex_Add(FVectorX,
COMplex_Mul( CFactor1,AVectorPQ[0]));
FVectorY = COMplex_Add(FVectorY,
COMplex_Mul( CFactor1,AVectorPQ[1]));
FVectorZ = COMplex_Add(FVectorZ,
COMplex_Mul( CFactor1,AVectorPQ[2]));
} /* if(RowNum!=1) */
}
}
ATheta = COMplex_Sub( COMplex_Add(Real_Mul((cth*cph),AVectorX),
Real_Mul((cth*sph),AVectorY)),
Real_Mul(sth,AVectorZ));
APhi = COMplex_Sub( Real_Mul( cph , AVectorY),
Real_Mul( sph , AVectorX) ) ;
FTheta = COMplex_Sub( COMplex_Add(Real_Mul((cth*cph),FVectorX),
Real_Mul((cth*sph),FVectorY)),
Real_Mul( sth,FVectorZ) );
FPhi = COMplex_Sub(Real_Mul(cph,FVectorY), Real_Mul(sph,FVectorX)) ;
FTheta = Real_Mul(1.0/ImpeDance, FTheta);
FPhi = Real_Mul(1.0/ImpeDance, FPhi);
ETheta = COMplex_Mul(CFactor,COMplex_Add(ATheta,FPhi));
EPhi = COMplex_Mul(CFactor,COMplex_Sub(APhi,FTheta));
tmp=1.0/(4.0*M_PI*RadDisObserv);
ETheta = Real_Mul(tmp,ETheta );
EPhi = Real_Mul(tmp, EPhi );
/********E filed computation ***********/
fprintf(OutF, "%e\t%e\n", COMplex_Abs(EPhi),COMplex_Abs(ETheta));
} /* end of ComputeFarEField( ) */
/********************* ComputeNonSingul() **********************/
void ComputeNonSingul(int *TrngleNode, double cth, double cph, double sth,
double sph, double ObsDist)
{
int QuadPoint;
double Rcp,SrcPoint[3];
complex V1, V2, Eprj;
for(QuadPoint=0;QuadPoint<=TotQuadPoint-1;QuadPoint++)
{
ComputeGaussQuadPoint(QuadPoint,TrngleNode,SrcPoint);
Rcp = (SrcPoint[0]*cph + SrcPoint[1]*sph)*sth + SrcPoint[2]*cth;
V1 = COMplex_Cmplx(0.0,(Rcp*WaveNumber));
V1 = COMplex_Expon(1.0,V1);
/* Following codes are quoted by Ali*/
V2 = COMplex_Expon(-1.0,COMplex_Cmplx(0.0,(ObsDist*WaveNumber)));
V1 = COMplex_Mul(V1,V2);
Eprj = Real_Mul(Wght[QuadPoint],V1);
InoPQ = COMplex_Add(InoPQ, Eprj);
IxiPQ = COMplex_Add(IxiPQ, Real_Mul(Qpnt[QuadPoint][0],Eprj) );
IetaPQ = COMplex_Add(IetaPQ, Real_Mul(Qpnt[QuadPoint][1],Eprj) );
}
IzetaPQ = COMplex_Sub(InoPQ,COMplex_Add(IxiPQ ,IetaPQ) );
}
/************** ComputeGaussQuadPoint() ********************************/
void ComputeGaussQuadPoint(int QuadPoint, int *TrngleNode, double *SrcPointCol)
{
SrcPointCol[0] = Qpnt[QuadPoint][0]*NodeCord[TrngleNode[0]][0]
+ Qpnt[QuadPoint][1]*NodeCord[TrngleNode[1]][0]
+ Qpnt[QuadPoint][2]*NodeCord[TrngleNode[2]][0];
SrcPointCol[1] = Qpnt[QuadPoint][0]*NodeCord[TrngleNode[0]][1]
+ Qpnt[QuadPoint][1]*NodeCord[TrngleNode[1]][1]
+ Qpnt[QuadPoint][2]*NodeCord[TrngleNode[2]][1];
SrcPointCol[2] = Qpnt[QuadPoint][0]*NodeCord[TrngleNode[0]][2]
+ Qpnt[QuadPoint][1]*NodeCord[TrngleNode[1]][2]
+ Qpnt[QuadPoint][2]*NodeCord[TrngleNode[2]][2];
} /* ComputeGaussQuadPoint( ) */
/***************** GetGaussQuadData() ***************************************/
void GetGaussQuadData()
{
int QuadPoint;
Wght[0]=0.112500000000000;
Qpnt[0][0]= Qpnt[0][1]=0.333333333333333;
Wght[1]=0.062969590272413;
Qpnt[1][0]=0.797426985353087;
Qpnt[1][1]=0.101286507323456;
Wght[2]=0.062969590272413;
Qpnt[2][0]=0.101286507323456;
Qpnt[2][1]=0.797426985353087;
Wght[3]=0.062969590272413;
Qpnt[3][0]= Qpnt[3][1]=0.101286507323456;
Wght[4]=0.066197076394253;
Qpnt[4][0]= Qpnt[4][1]=0.470142064105115;
Wght[5]=0.066197076394253;
Qpnt[5][0]=0.470142064105115;
Qpnt[5][1]=0.059715871789770;
Wght[6]=0.066197076394253;
Qpnt[6][0]=0.059715871789770;
Qpnt[6][1]=0.470142064105115;
for(QuadPoint=0;QuadPoint<=TotQuadPoint-1;QuadPoint++)
Qpnt[QuadPoint][2] = 1.0 - Qpnt[QuadPoint][0] - Qpnt[QuadPoint][1];
} /* end of GetGaussQuadData() */
double Sign(int Value)
{
return (Value>0 ? 1.0:-1.0);
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -