c_geometry.cpp

矩量法仿真电磁辐射和散射的源代码（c++）

		njun1= jsno;
	
		jcox= icon2[ix];
		if ( jcox > PCHCON)
			jcox= i;
	
		jend=1;
		iend=1;
		_sig=-1.;
	
	} /* do */
	while( jcox != 0 );
	
	njun2= jsno- njun1;
	jsnop= jsno;
	jco[jsnop]= i;
	
	nec_float d = pi()* segment_length[ix];
	sdh = sin(d);
	cdh = cos(d);
	sd = 2.0 * sdh * cdh;
	cd= cdh*cdh - sdh*sdh;
	
	if ( d <= 0.015)
	{
		omc = 4.0* d*d;
		omc = ((1.3888889e-3* omc-4.1666666667e-2)* omc+.5)* omc;
	}
	else
		omc = 1.0 - cd;
	
	ap=1./( log(1./( pi()* segment_radius[ix]))-.577215664);
	aj= ap;
	
	if ( njun1 == 0)
	{
		if ( njun2 == 0)
		{
			bx[jsnop]=0.;
		
			if ( icap == 0)
				xxi=0.;
			else
			{
				qp= pi()* segment_radius[ix];
				xxi= qp* qp;
				xxi= qp*(1.-.5* xxi)/(1.- xxi);
			}
		
			cx[jsnop]=1./( cdh- xxi* sdh);
			jsno= jsnop+1;
			ax[jsnop]=-1.;
			return;
		} /* if ( njun2 == 0) */
	
		if ( icap == 0)
			xxi=0.;
		else
		{
			qp= pi()* segment_radius[ix];
			xxi= qp* qp;
			xxi= qp*(1.-.5* xxi)/(1.- xxi);
		}
	
		qp=-( omc+ xxi* sd)/( sd*( ap+ xxi* pp)+ cd*( xxi* ap- pp));
		d= cd- xxi* sd;
		bx[jsnop]=( sdh+ ap* qp*( cdh- xxi* sdh))/ d;
		cx[jsnop]=( cdh+ ap* qp*( sdh+ xxi* cdh))/ d;
	
		for( iend = 0; iend < njun2; iend++ )
		{
			ax[iend]=-ax[iend]* qp;
			bx[iend]= bx[iend]* qp;
			cx[iend]=- cx[iend]* qp;
		}
	
		jsno= jsnop+1;
		ax[jsnop]=-1.;
		return;
	} /* if ( njun1 == 0) */
	
	if ( njun2 == 0)
	{
		if ( icap == 0)
			xxi=0.;
		else
		{
			qm= pi()* segment_radius[ix];
			xxi= qm* qm;
			xxi= qm*(1.-.5* xxi)/(1.- xxi);
		}
	
		qm=( omc+ xxi* sd)/( sd*( aj- xxi* pm)+ cd*( pm+ xxi* aj));
		d= cd- xxi* sd;
		bx[jsnop]=( aj* qm*( cdh- xxi* sdh)- sdh)/ d;
		cx[jsnop]=( cdh- aj* qm*( sdh+ xxi* cdh))/ d;
	
		for( iend = 0; iend < njun1; iend++ )
		{
			ax[iend]= ax[iend]* qm;
			bx[iend]= bx[iend]* qm;
			cx[iend]= cx[iend]* qm;
		}
	
		jsno= jsnop+1;
		ax[jsnop]=-1.;
		return;
	
	} /* if ( njun2 == 0) */
	
	qp= sd*( pm* pp+ aj* ap)+ cd*( pm* ap- pp* aj);
	qm=( ap* omc- pp* sd)/ qp;
	qp=-( aj* omc+ pm* sd)/ qp;
	bx[jsnop]=( aj* qm+ ap* qp)* sdh/ sd;
	cx[jsnop]=( aj* qm- ap* qp)* cdh/ sd;
	
	for( iend = 0; iend < njun1; iend++ )
	{
		ax[iend]= ax[iend]* qm;
		bx[iend]= bx[iend]* qm;
		cx[iend]= cx[iend]* qm;
	}
	
	jend= njun1;
	for( iend = jend; iend < jsno; iend++ )
	{
		ax[iend]=- ax[iend]* qp;
		bx[iend]= bx[iend]* qp;
		cx[iend]=- cx[iend]* qp;
	}
	
	jsno= jsnop+1;
	ax[jsnop]=-1.;
}



/* compute the components of all basis functions on segment j */
void c_geometry::trio( int j )
{
	int jcox, jcoxx, jsnox, jx, jend=0, iend=0;
	
	jsno=0;
	jx = j-1;
	jcox= icon1[jx];
	jcoxx = jcox-1;
	
	if ( jcox <= PCHCON)
	{
		jend=-1;
		iend=-1;
	}
	
	if ( (jcox == 0) || (jcox > PCHCON) )
	{
		jcox= icon2[jx];
		jcoxx = jcox-1;
	
		if ( jcox <= PCHCON)
		{
			jend=1;
			iend=1;
		}
	
		if ( jcox == 0 || (jcox > PCHCON) )
		{
			jsnox = jsno;
			jsno++;
		
			/* Allocate to connections buffers */
			if ( jsno >= maxcon )
			{
				maxcon = jsno +1;
				jco.resize(maxcon);
				ax.resize(maxcon);
				bx.resize(maxcon);
				cx.resize(maxcon);
			}
		
			sbf( j, j, &ax[jsnox], &bx[jsnox], &cx[jsnox]);
			jco[jsnox]= j;
			return;
		}
	} /* if ( (jcox == 0) || (jcox > PCHCON) ) */
	
	do
	{
		if ( jcox < 0 )
			jcox=- jcox;
		else
			jend=- jend;
		jcoxx = jcox-1;
	
		if ( jcox != j)
		{
			jsnox = jsno;
			jsno++;
		
			/* Allocate to connections buffers */
			if ( jsno >= maxcon )
			{
				maxcon = jsno +1;
				jco.resize(maxcon );
				ax.resize(maxcon);
				bx.resize(maxcon);
				cx.resize(maxcon);
			}
		
			sbf( jcox, j, &ax[jsnox], &bx[jsnox], &cx[jsnox]);
			jco[jsnox]= jcox;
		
			if ( jend != 1)
				jcox= icon1[jcoxx];
			else
				jcox= icon2[jcoxx];
		
			if ( jcox == 0 )
			{
				nec_exception* nex = new nec_exception("TRIO - SEGMENT CONNENTION ERROR FOR SEGMENT ");
				nex->append(j);
				throw nex;
			}
			else
				continue;
		} /* if ( jcox != j) */
	
		if ( iend == 1)
			break;
	
		jcox= icon2[jx];
	
		if ( jcox > PCHCON)
			break;
	
		jend=1;
		iend=1;
	} /* do */
	while( jcox != 0 );
	
	jsnox = jsno;
	jsno++;
	
	/* Allocate to connections buffers */
	if ( jsno >= maxcon )
	{
		maxcon = jsno +1;
		jco.resize(maxcon );
		ax.resize(maxcon);
		bx.resize(maxcon);
		cx.resize(maxcon);
	}
	
	sbf( j, j, &ax[jsnox], &bx[jsnox], &cx[jsnox]);
	jco[jsnox]= j;
}



/*! \brief compute component of basis function i on segment is.
This is a version of the tbf() method that does not store the basis functions
*/
void c_geometry::sbf( int i, int is, nec_float *aa, nec_float *bb, nec_float *cc )
{
	int local_jsno; // this parameter is renamed because it shadows the member variable of the same name
	int ix, june, jcox, jcoxx, jend, iend, njun1=0, njun2;
	nec_float d, sig, pp, sdh, cdh, sd, omc, aj, pm=0, cd, ap, qp, qm, xxi;
	
	*aa=0.;
	*bb=0.;
	*cc=0.;
	june=0;
	local_jsno=0;
	pp=0.;
	ix=i-1;
	
	jcox= icon1[ix];
	if ( jcox > PCHCON)
		jcox= i;
	jcoxx = jcox-1;
	
	jend=-1;
	iend=-1;
	sig=-1.;
	
	do
	{
		if ( jcox != 0 )
		{
			if ( jcox < 0 )
				jcox=- jcox;
			else
			{
				sig=- sig;
				jend=- jend;
			}
		
			jcoxx = jcox-1;
			local_jsno++;
			d= pi()* segment_length[jcoxx];
			sdh= sin( d);
			cdh= cos( d);
			sd=2.* sdh* cdh;
		
			if ( d <= 0.015)
			{
				omc=4.* d* d;
				omc=((1.3888889e-3* omc -4.1666666667e-2)* omc +.5)* omc;
			}
			else
				omc=1.- cdh* cdh+ sdh* sdh;
		
			aj=1./( log(1./( pi()* segment_radius[jcoxx]))-.577215664);
			pp -= omc/ sd* aj;
		
			if ( jcox == is)
			{
				*aa= aj/ sd* sig;
				*bb= aj/(2.* cdh);
				*cc=- aj/(2.* sdh)* sig;
				june= iend;
			}
		
			if ( jcox != i )
			{
				if ( jend != 1)
					jcox= icon1[jcoxx];
				else
					jcox= icon2[jcoxx];
			
				if ( abs(jcox) != i )
				{
					if ( jcox == 0 )
					{
						nec_exception* nex = new nec_exception("SBF - SEGMENT CONNECTION ERROR FOR SEGMENT ");
						nex->append(i);
						throw nex;
					}
					else
						continue;
				}	
			} /* if ( jcox != i ) */
			else
			if ( jcox == is)
				*bb=- *bb;
		
			if ( iend == 1)
				break;
		} /* if ( jcox != 0 ) */
	
		pm=- pp;
		pp=0.;
		njun1= local_jsno;
	
		jcox= icon2[ix];
		if ( jcox > PCHCON)
			jcox= i;
	
		jend=1;
		iend=1;
		sig=-1.;
	
	} /* do */
	while( jcox != 0 );
	
	njun2= local_jsno- njun1;
	d= pi()* segment_length[ix];
	sdh= sin( d);
	cdh= cos( d);
	sd=2.* sdh* cdh;
	cd= cdh* cdh- sdh* sdh;
	
	if ( d <= 0.015)
	{
		omc=4.* d* d;
		omc=((1.3888889e-3* omc -4.1666666667e-2)* omc +.5)* omc;
	}
	else
		omc=1.- cd;
	
	ap=1./( log(1./( pi()* segment_radius[ix])) -.577215664);
	aj= ap;
	
	if ( njun1 == 0)
	{
		if ( njun2 == 0)
		{
			*aa =-1.;
			qp= pi()* segment_radius[ix];
			xxi= qp* qp;
			xxi= qp*(1.-.5* xxi)/(1.- xxi);
			*cc=1./( cdh- xxi* sdh);
				return;
		}
	
		qp= pi()* segment_radius[ix];
		xxi= qp* qp;
		xxi= qp*(1.-.5* xxi)/(1.- xxi);
		qp=-( omc+ xxi* sd)/( sd*( ap+ xxi* pp)+ cd*( xxi* ap- pp));
	
		if ( june == 1)
		{
			*aa=- *aa* qp;
			*bb=  *bb* qp;
			*cc=- *cc* qp;
			if ( i != is)
				return;
		}
	
		*aa -= 1.;
		d = cd - xxi * sd;
		*bb += (sdh + ap * qp * (cdh - xxi * sdh)) / d;
		*cc += (cdh + ap * qp * (sdh + xxi * cdh)) / d;
		return;
	
	} /* if ( njun1 == 0) */
	
	if ( njun2 == 0)
	{
		qm= pi()* segment_radius[ix];
		xxi= qm* qm;
		xxi= qm*(1.-.5* xxi)/(1.- xxi);
		qm=( omc+ xxi* sd)/( sd*( aj- xxi* pm)+ cd*( pm+ xxi* aj));
	
		if ( june == -1)
		{
			*aa= *aa* qm;
			*bb= *bb* qm;
			*cc= *cc* qm;
			if ( i != is)
				return;
		}
	
		*aa -= 1.;
		d= cd- xxi* sd;
		*bb += ( aj* qm*( cdh- xxi* sdh)- sdh)/ d;
		*cc += ( cdh- aj* qm*( sdh+ xxi* cdh))/ d;
		return;	
	} /* if ( njun2 == 0) */
	
	qp= sd*( pm* pp+ aj* ap)+ cd*( pm* ap- pp* aj);
	qm=( ap* omc- pp* sd)/ qp;
	qp=-( aj* omc+ pm* sd)/ qp;
	
	if ( june != 0 )
	{
		if ( june < 0 )
		{
			*aa= *aa* qm;
			*bb= *bb* qm;
			*cc= *cc* qm;
		}
		else
		{
			*aa=- *aa* qp;
			*bb= *bb* qp;
			*cc=- *cc* qp;
		}
	
		if ( i != is)
			return;
	} /* if ( june != 0 ) */
	
	*aa -= 1.;
	*bb += ( aj* qm+ ap* qp)* sdh/ sd;
	*cc += ( aj* qm- ap* qp)* cdh/ sd;
}


/*
	get_current_coefficients computes coefficients of the:
		constant 	[air, aii]
		sine		[bir, bii]
		cosine		[cir, cii]
	terms in the current interpolation functions for the current vector curx.
 */
void c_geometry::get_current_coefficients(nec_float wavelength, complex_array& curx,
			real_array& air, real_array& aii,
			real_array& bir, real_array& bii,
			real_array& cir, real_array& cii,
			complex_array& vqds, int nqds,
			int_array& iqds)
{
	static nec_complex s_CCJ(0.0,-0.01666666667);
	
	nec_float ar, ai, sh;
	nec_complex cs1, cs2;
	
	air.fill(0,n,0.0);
	aii.fill(0,n,0.0);
	bir.fill(0,n,0.0);
	bii.fill(0,n,0.0);
	cir.fill(0,n,0.0);
	cii.fill(0,n,0.0);
	
	if ( n != 0)
	{	
		for (int i = 0; i < n; i++ )
		{
			ar= real( curx[i]);
			ai= imag( curx[i]);
			tbf( i+1, 1 );
		
			for (int jx = 0; jx < jsno; jx++ )
			{
				int j = jco[jx]-1;
				air[j] += ax[jx]* ar;
				aii[j] += ax[jx]* ai;
				bir[j] += bx[jx]* ar;
				bii[j] += bx[jx]* ai;
				cir[j] += cx[jx]* ar;
				cii[j] += cx[jx]* ai;
			}
		} /* for( i = 0; i < n; i++ ) */
	
		for (int is = 0; is < nqds; is++ )
		{
			int i= iqds[is]-1;
			int jx= icon1[i];
			icon1[i]=0;
			tbf(i+1,0);
			icon1[i]= jx;
			sh = segment_length[i]*.5;
			
			
			nec_complex curd = s_CCJ * 
				vqds[is]/( (log(2.* sh/ segment_radius[i])-1.)*
				(bx[jsno-1]* cos(two_pi() * sh)+ cx[jsno-1]* sin(two_pi() * sh))* wavelength );
			ar = real( curd);
			ai = imag( curd);
		
			for ( jx = 0; jx < jsno; jx++ )
			{
				int j = jco[jx]-1;
				air[j] += ax[jx]* ar;
				aii[j] += ax[jx]* ai;
				bir[j] += bx[jx]* ar;
				bii[j] += bx[jx]* ai;
				cir[j] += cx[jx]* ar;
				cii[j] += cx[jx]* ai;
			}
		
		} /* for( is = 0; is < nqds; is++ ) */
	
		for (int i = 0; i < n; i++ )
			curx[i]= nec_complex( air[i]+cir[i], aii[i]+cii[i] );
	
	} /* if ( n != 0) */
	
	if ( m == 0)
		return;
	
	/* convert surface currents from */
	/* t1,t2 components to x,y,z components */
	int jco1 = n_plus_2m;
	int jco2 = jco1 + m;
	
	for (int i = 1; i <= m; i++ )
	{
		jco1 -= 2;
		jco2 -= 3;
		cs1= curx[jco1];
		cs2= curx[jco1+1];
		curx[jco2]  = cs1* t1x[m-i]+ cs2* t2x[m-i];
		curx[jco2+1]= cs1* t1y[m-i]+ cs2* t2y[m-i];
		curx[jco2+2]= cs1* t1z[m-i]+ cs2* t2z[m-i];
	}
}


void c_geometry::frequency_scale(nec_float freq_mhz)
{
	DEBUG_TRACE("frequency_scale(" << freq_mhz << ")");
	nec_float fr = freq_mhz/ CVEL;

	for (int i = 0; i < n; i++ )
	{
		x[i]= x_unscaled[i]* fr;
		y[i]= y_unscaled[i]* fr;
		z[i]= z_unscaled[i]* fr;
		segment_length[i]= si_unscaled[i]* fr;
		segment_radius[i]= bi_unscaled[i]* fr;
	}

	nec_float fr2 = fr*fr;
	for (int i = 0; i < m; i++ )
	{
		px[i]= px_unscaled[i]* fr;
		py[i]= py_unscaled[i]* fr;
		pz[i]= pz_unscaled[i]* fr;
		pbi[i]= pbi_unscaled[i]* fr2;
	}
}


void c_geometry::fflds(nec_float rox, nec_float roy, nec_float roz,
	complex_array& scur, 
	nec_complex *in_ex, nec_complex *in_ey, nec_complex *in_ez )
{
	static nec_complex _const4(0.0,+188.365);
	
	// From FORTRAN common block 
	// EQUIVALENCE (XS,X), (YS,Y), (ZS,Z), (S,BI), (CONS,CONSX)
	nec_complex ex(cplx_00());
	nec_complex ey(cplx_00());
	nec_complex ez(cplx_00());
	
	for (int i = 0; i < m; i++ )
	{
		nec_float arg = patch_angle(i,rox,roy,roz);
		nec_complex ct = cplx_exp(arg) * pbi[i];
		int k = 3*i;
		ex += scur[k]* ct;
		ey += scur[k+1]* ct;
		ez += scur[k+2]* ct;
	}
	
	nec_complex ct = rox*ex+ roy*ey+ roz*ez;
	
	*in_ex = _const4*(ct*rox - ex);
	*in_ey = _const4*(ct*roy - ey);
	*in_ez = _const4*(ct*roz - ez);
}


int c_geometry::test_ek_approximation(int seg1, int seg2)
{
	nec_float segment_ratio = segment_radius[seg2] / segment_radius[seg1