c_geometry.cpp

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

	
	ax.resize(maxcon);
	bx.resize(maxcon);
	cx.resize(maxcon);
}

/* arc generates segment geometry data for an arc of segment_count segments */
void c_geometry::arc( int tag_id, int segment_count, nec_float rada,
    nec_float ang1, nec_float ang2, nec_float rad )
{
	int istart = n;
	n += segment_count;
	np= n;
	mp= m;
	m_ipsym=0;
	
	if ( segment_count < 1)
		return;
	
	if ( fabs(ang2 - ang1) > 360.0)
	{
		throw new nec_exception("ERROR -- ARC ANGLE EXCEEDS 360 DEGREES");
	}
	
	/* Reallocate tags buffer */
	segment_tags.resize(n+m);

	/* Reallocate wire buffers */
	x.resize(n);
	y.resize(n);
	z.resize(n);
	x2.resize(n);
	y2.resize(n);
	z2.resize(n);
	segment_radius.resize(n);

	nec_float ang = degrees_to_rad(ang1);
	nec_float dang = degrees_to_rad(ang2- ang1) / segment_count;
	nec_float xs1= rada * cos(ang);
	nec_float zs1= rada * sin(ang);

	for(int i = istart; i < n; i++ )
	{
		ang += dang;
		nec_float xs2 = rada * cos(ang);
		nec_float zs2 = rada * sin(ang);
		x[i]= xs1;
		y[i]=0.;
		z[i]= zs1;
		x2[i]= xs2;
		y2[i]=0.;
		z2[i]= zs2;
		xs1= xs2;
		zs1= zs2;
		segment_radius[i]= rad;
		segment_tags[i]= tag_id;
	} /* for( i = ist; i < n; i++ ) */
}


/*-----------------------------------------------------------------------*/

/*! \brief patch generates and modifies patch geometry data.
*/
void c_geometry::patch( int nx, int ny,
    nec_float ax1, nec_float ay1, nec_float az1,
    nec_float ax2, nec_float ay2, nec_float az2,
    nec_float ax3, nec_float ay3, nec_float az3,
    nec_float ax4, nec_float ay4, nec_float az4 )
{
	int mi, ntp, iy, ix;
	nec_float s1x=0., s1y=0., s1z=0., s2x=0., s2y=0., s2z=0., xst=0.;
	nec_float znv, xnv, ynv, xa, xn2, yn2, zn2, salpn, xs, ys, zs, xt, yt, zt;
	
	/* new patches.  for nx=0, ny=1,2,3,4 patch is (respectively) */;
	/* arbitrary, rectagular, triangular, or quadrilateral. */
	/* for nx and ny  > 0 a rectangular surface is produced with */
	/* nx by ny rectangular patches. */
	
	m++;
	mi= m-1;
	
	/* Reallocate patch buffers */
	px.resize(m);
	py.resize(m);
	pz.resize(m);
	t1x.resize(m);
	t1y.resize(m);
	t1z.resize(m);
	t2x.resize(m);
	t2y.resize(m);
	t2z.resize(m);
	pbi.resize(m);
	psalp.resize(m);
	
	if ( nx > 0)
		ntp=2;
	else
		ntp= ny;
	
	if ( ntp <= 1)
	{
		px[mi]= ax1;
		py[mi]= ay1;
		pz[mi]= az1;
		pbi[mi]= az2;
		znv= cos( ax2);
		xnv= znv* cos( ay2);
		ynv= znv* sin( ay2);
		znv= sin( ax2);
		xa= sqrt( xnv* xnv+ ynv* ynv);
	
		if ( xa >= 1.0e-6)
		{
			t1x[mi]=- ynv/ xa;
			t1y[mi]= xnv/ xa;
			t1z[mi]=0.;
		}
		else
		{
			t1x[mi]=1.;
			t1y[mi]=0.;
			t1z[mi]=0.;
		}
	
	} /* if ( ntp <= 1) */
	else
	{
		s1x= ax2- ax1;
		s1y= ay2- ay1;
		s1z= az2- az1;
		s2x= ax3- ax2;
		s2y= ay3- ay2;
		s2z= az3- az2;
	
		if ( nx != 0)
		{
			s1x= s1x/ nx;
			s1y= s1y/ nx;
			s1z= s1z/ nx;
			s2x= s2x/ ny;
			s2y= s2y/ ny;
			s2z= s2z/ ny;
		}
	
		xnv= s1y* s2z- s1z* s2y;
		ynv= s1z* s2x- s1x* s2z;
		znv= s1x* s2y- s1y* s2x;
		xa= sqrt( xnv* xnv+ ynv* ynv+ znv* znv);
		xnv= xnv/ xa;
		ynv= ynv/ xa;
		znv= znv/ xa;
		xst= sqrt( s1x* s1x+ s1y* s1y+ s1z* s1z);
		t1x[mi]= s1x/ xst;
		t1y[mi]= s1y/ xst;
		t1z[mi]= s1z/ xst;
	
		if ( ntp <= 2)
		{
			px[mi]= ax1+.5*( s1x+ s2x);
			py[mi]= ay1+.5*( s1y+ s2y);
			pz[mi]= az1+.5*( s1z+ s2z);
			pbi[mi]= xa;
		}
		else
		{
			if ( ntp != 4)
			{
				px[mi]=( ax1+ ax2+ ax3)/3.;
				py[mi]=( ay1+ ay2+ ay3)/3.;
				pz[mi]=( az1+ az2+ az3)/3.;
				pbi[mi]=.5* xa;
			}
			else
			{
				s1x= ax3- ax1;
				s1y= ay3- ay1;
				s1z= az3- az1;
				s2x= ax4- ax1;
				s2y= ay4- ay1;
				s2z= az4- az1;
				xn2= s1y* s2z- s1z* s2y;
				yn2= s1z* s2x- s1x* s2z;
				zn2= s1x* s2y- s1y* s2x;
				xst= sqrt( xn2* xn2+ yn2* yn2+ zn2* zn2);
				salpn=1./(3.*( xa+ xst));
				px[mi]=( xa*( ax1+ ax2+ ax3)+ xst*( ax1+ ax3+ ax4))* salpn;
				py[mi]=( xa*( ay1+ ay2+ ay3)+ xst*( ay1+ ay3+ ay4))* salpn;
				pz[mi]=( xa*( az1+ az2+ az3)+ xst*( az1+ az3+ az4))* salpn;
				pbi[mi]=.5*( xa+ xst);
				s1x=( xnv* xn2+ ynv* yn2+ znv* zn2)/ xst;
			
				if ( s1x <= 0.9998)
				{
					throw new nec_exception("ERROR -- CORNERS OF QUADRILATERAL PATCH DO NOT LIE IN A PLANE");
				}
			} /* if ( ntp != 4) */
		} /* if ( ntp <= 2) */
	
	} /* if ( ntp <= 1) */
	
	t2x[mi]= ynv* t1z[mi]- znv* t1y[mi];
	t2y[mi]= znv* t1x[mi]- xnv* t1z[mi];
	t2z[mi]= xnv* t1y[mi]- ynv* t1x[mi];
	psalp[mi]=1.;
	
	if ( nx != 0)
	{
		m += nx*ny-1;
	
		/* Reallocate patch buffers */
		px.resize(m);
		py.resize(m);
		pz.resize(m);
		t1x.resize(m);
		t1y.resize(m);
		t1z.resize(m);
		t2x.resize(m);
		t2y.resize(m);
		t2z.resize(m);
		pbi.resize(m);
		psalp.resize(m);
	
		xn2= px[mi]- s1x- s2x;
		yn2= py[mi]- s1y- s2y;
		zn2= pz[mi]- s1z- s2z;
		xs= t1x[mi];
		ys= t1y[mi];
		zs= t1z[mi];
		xt= t2x[mi];
		yt= t2y[mi];
		zt= t2z[mi];
	
		for( iy = 0; iy < ny; iy++ )
		{
			xn2 += s2x;
			yn2 += s2y;
			zn2 += s2z;
		
			for( ix = 1; ix <= nx; ix++ )
			{
				xst= (nec_float)ix;
				px[mi]= xn2+ xst* s1x;
				py[mi]= yn2+ xst* s1y;
				pz[mi]= zn2+ xst* s1z;
				pbi[mi]= xa;
				psalp[mi]=1.;
				t1x[mi]= xs;
				t1y[mi]= ys;
				t1z[mi]= zs;
				t2x[mi]= xt;
				t2y[mi]= yt;
				t2z[mi]= zt;
				mi++;
			} /* for( ix = 0; ix < nx; ix++ ) */
		} /* for( iy = 0; iy < ny; iy++ ) */
	
	} /* if ( nx != 0) */
	
	m_ipsym=0;
	np= n;
	mp= m;
}


/*!\brief Divide a patch into four (was subph).
Used when a patch is connected to a wire. The patch
nx is  divided into 4 patches that become nx, nx+1, nx+2, nx+3. The other
patches are shifted to make room for the three new patches.

\param nx The index of the patch to divide (starting at 1)
*/
void c_geometry::divide_patch(int nx)

{	
	m += 3;

	px.resize(m);
	py.resize(m);
	pz.resize(m);
	t1x.resize(m);
	t1y.resize(m);
	t1z.resize(m);
	t2x.resize(m);
	t2y.resize(m);
	t2z.resize(m);
	pbi.resize(m);
	psalp.resize(m);
	
	/* Shift patches to make room for new ones */
	for (int iy = m-1; iy > nx; iy--)
	{
		int old_index = iy - 3;
		px[iy] = px[old_index];
		py[iy] = py[old_index];
		pz[iy] = pz[old_index];
		pbi[iy] = pbi[old_index];
		psalp[iy] = psalp[old_index];
		t1x[iy] = t1x[old_index];
		t1y[iy] = t1y[old_index];
		t1z[iy] = t1z[old_index];
		t2x[iy] = t2x[old_index];
		t2y[iy] = t2y[old_index];
		t2z[iy] = t2z[old_index];
	}
	
	/* divide patch for connection */
	int patch_index = nx-1;
	nec_float xs = px[patch_index];
	nec_float ys = py[patch_index];
	nec_float zs = pz[patch_index];
	nec_float xa = pbi[patch_index]/4.;
	nec_float xst = sqrt(xa)/2.;
	nec_float s1x = t1x[patch_index];
	nec_float s1y = t1y[patch_index];
	nec_float s1z = t1z[patch_index];
	nec_float s2x = t2x[patch_index];
	nec_float s2y = t2y[patch_index];
	nec_float s2z = t2z[patch_index];
	nec_float saln = psalp[patch_index];
	nec_float xt = xst;
	nec_float yt = xst;
	
	int new_index = patch_index;
	
	/* Generate the four new patches */
	for (int ix = 1; ix <= 4; ix++ )
	{
		px[new_index]= xs+ xt* s1x+ yt* s2x;
		py[new_index]= ys+ xt* s1y+ yt* s2y;
		pz[new_index]= zs+ xt* s1z+ yt* s2z;
		pbi[new_index]= xa;
		t1x[new_index]= s1x;
		t1y[new_index]= s1y;
		t1z[new_index]= s1z;
		t2x[new_index]= s2x;
		t2y[new_index]= s2y;
		t2z[new_index]= s2z;
		psalp[new_index]= saln;
	
		if (2 == ix)
			yt = -yt;
	
		if ( (ix == 1) || (ix == 3) )
			xt = -xt;
	
		new_index++;
	}
	
	/* Readjust the mp patch index to account for the added patches */
	if (nx <= mp)
		mp += 3;
}

/*-----------------------------------------------------------------------

	Read Geometry Data from a Card

-------------------------------------------------------------------------*/

void c_geometry::read_geometry_card(FILE* input_fp,  char *gm,
	int *in_i1, int *in_i2,
	nec_float *in_x1, nec_float *in_y1, nec_float *in_z1,
	nec_float *in_x2, nec_float *in_y2, nec_float *in_z2,
	nec_float *in_rad )
{
	char line_buf[134];
	int i, line_idx;
	int n_integer_params = 2, n_float_params = 7;
	int integer_params[2] = { 0, 0 };
	nec_float real_params[7] = { 0., 0., 0., 0., 0., 0., 0. };
	
	/* read a line from input file */
	load_line( line_buf, input_fp );
	
	/* get line length */
	int line_length = (int)strlen( line_buf );
	
	/* abort if card's mnemonic too short or missing */
	if ( line_length < 2 )
	{
		nec_exception* nex = new nec_exception("GEOMETRY DATA CARD ERROR:");
		nex->append(" CARD'S MNEMONIC CODE TOO SHORT OR MISSING.");
		throw nex;
	}
	
	/* extract card's mnemonic code */
	strncpy( gm, line_buf, 2 );
	gm[2] = '\0';
	
	/* Exit if "XT" command read (for testing) */
	if ( strcmp( gm, "XT" ) == 0 )
	{
			nec_exception* nex = new nec_exception("Exiting after an \"XT\" command in read_geometry_card()");
			throw nex;
	}
	
	/* Return if only mnemonic on card */
	if ( line_length == 2 )
	{
		*in_i1 = *in_i2 = 0;
		*in_x1 = *in_y1 = *in_z1 = *in_x2 = *in_y2 = *in_z2 = *in_rad = 0.0;
		return;
	}
	
	/* read integers from line */
	line_idx = 1;
	for( i = 0; i < n_integer_params; i++ )
	{
		/* Find first numerical character */
		while( ((line_buf[++line_idx] <  '0')  ||
			(line_buf[  line_idx] >  '9')) &&
			(line_buf[  line_idx] != '+')  &&
			(line_buf[  line_idx] != '-') )
		if ( (line_buf[line_idx] == '\0') )
		{
			*in_i1= integer_params[0];
			*in_i2= integer_params[1];
			*in_x1= real_params[0];
			*in_y1= real_params[1];
			*in_z1= real_params[2];
			*in_x2= real_params[3];
			*in_y2= real_params[4];
			*in_z2= real_params[5];
			*in_rad= real_params[6];
			return;
		}
	
		/* read an integer from line */
		integer_params[i] = atoi( &line_buf[line_idx] );
	
		/* traverse numerical field to next ' ' or ',' or '\0' */
		line_idx--;
		while( 	(line_buf[++line_idx] != ' ') &&
				(line_buf[  line_idx] != ',') &&
				(line_buf[  line_idx] != '\0') )
		{
			/* test for non-numerical characters */
			if ( 	((line_buf[line_idx] <  '0')  ||
					(line_buf[line_idx] >  '9')) &&
					(line_buf[line_idx] != '+')  &&
					(line_buf[line_idx] != '-') )
			{
				nec_stop(
					"GEOMETRY DATA CARD \"%s\" ERROR:"
					"\n  NON-NUMERICAL CHARACTER '%c' IN INTEGER FIELD AT CHAR. %d\n",
					gm, line_buf[line_idx], (line_idx+1)  );
			}
		} /* while( (line_buff[++line_idx] ... */
	
		/* Return on end of line */
		if ( line_buf[line_idx] == '\0' )
		{
			*in_i1= integer_params[0];
			*in_i2= integer_params[1];
			*in_x1= real_params[0];
			*in_y1= real_params[1];
			*in_z1= real_params[2];
			*in_x2= real_params[3];
			*in_y2= real_params[4];
			*in_z2= real_params[5];
			*in_rad= real_params[6];
			return;
		}
	
	} /* for( i = 0; i < n_integer_params; i++ ) */
	
	/* read nec_floats from line */
	for( i = 0; i < n_float_params; i++ )
	{
		/* Find first numerical character */
		while( ((line_buf[++line_idx] <  '0')  ||
			(line_buf[  line_idx] >  '9')) &&
			(line_buf[  line_idx] != '+')  &&
			(line_buf[  line_idx] != '-')  &&
			(line_buf[  line_idx] != '.') )
		if ( (line_buf[line_idx] == '\0') )
		{
			*in_i1= integer_params[0];
			*in_i2= integer_params[1];
			*in_x1= real_params[0];
			*in_y1= real_params[1];
			*in_z1= real_params[2];
			*in_x2= real_params[3];
			*in_y2= real_params[4];
			*in_z2= real_params[5];
			*in_rad= real_params[6];
			return;
		}
	
		/* read a nec_float from line */
		real_params[i] = atof( &line_buf[line_idx] );
	
		/* traverse numerical field to next ' ' or ',' or '\0' */
		line_idx--;
		while(	(line_buf[++line_idx] != ' ') &&
				(line_buf[  line_idx] != ',') &&
				(line_buf[  line_idx] != '\0') )
		{
			/* test for non-numerical characters */
			if ( ((line_buf[line_idx] <  '0')  ||
				(line_buf[line_idx] >  '9')) &&
				(line_buf[line_idx] != '.')  &&
				(line_buf[line_idx] != '+')  &&
				(line_buf[line_idx] != '-')  &&
				(line_buf[line_idx] != 'E')  &&
				(line_buf[line_idx] != 'e') )
			{
				nec_stop(
					"\n  GEOMETRY DATA CARD \"%s\" ERROR:"
					"\n  NON-NUMERICAL CHARACTER '%c' IN FLOAT FIELD AT CHAR. %d.\n",
					gm, line_buf[line_idx], (line_idx+1) );
			}
		} /* while( (line_buff[++line_idx] ... */
	
		/* Return on end of line */
		if ( line_buf[line_idx] == '\0' )
		{
			*in_i1= integer_params[0];
			*in_i2= integer_params[1];
			*in_x1= real_params[0];
			*in_y1= real_params[1];
			*in_z1= real_params[2];
			*in_x2= real_params[3];
			*in_y2= real_params[4];
			*in_z2= real_params[5];
			*in_rad= real_params[6];
			return;
		}
	
	} /* for( i = 0; i < n_float_params; i++ ) */
	
	*in_i1  = integer_params[0];
	*in_i2  = integer_params[1];
	*in_x1  = real_params[0];
	*in_y1  = real_params[1];
	*in_z1  = real_params[2];
	*in_x2  = real_params[3];
	*in_y2  = real_params[4];
	*in_z2  = real_params[5];
	*in_rad = real_params[6];
}



/* compute basis function i */
void c_geometry::tbf( int i, int icap )
{
	int jcoxx, njun1=0, njun2, jsnop, jsnox;
	nec_float sdh, cdh, sd, omc, aj, pm=0, cd, ap, qp, qm, xxi;
	
	jsno=0;
	
	nec_float pp = 0.0;
	int ix = i-1;
	int jcox = icon1[ix];
	
	if ( jcox > PCHCON)
		jcox= i;
	
	int jend = -1;
	int iend = -1;
	
	nec_float _sig = -1.0;
	do
	{
		if ( jcox != 0 )
		{
			if ( jcox < 0 )
				jcox=- jcox;
			else
			{
				_sig = -_sig;
				jend=- jend;
			}
		
			jcoxx = jcox-1;
			jsno++;
			jsnox = jsno-1;
			jco[jsnox]= jcox;
			nec_float d = pi() * segment_length[jcoxx];
			sdh = sin(d);
			cdh = cos(d);
			sd = 2.0 * 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= pp- omc/ sd* aj;
			ax[jsnox]= aj/ sd* _sig;
			bx[jsnox]= aj/(2.* cdh);
			cx[jsnox]=- aj/(2.* sdh)* _sig;
		
			if ( jcox != i)
			{
				if ( jend == 1)
					jcox= icon2[jcoxx];
				else
					jcox= icon1[jcoxx];
			
				if ( abs(jcox) != i )
				{
					if ( jcox != 0 )
						continue;
					else
					{
						nec_exception* nex = new nec_exception("TBF - SEGMENT CONNECTION ERROR FOR SEGMENT ");
						nex->append(i);
						throw nex;
					}
				}
			} /* if ( jcox != i) */
			else
			bx[jsnox] =- bx[jsnox];
		
			if ( iend == 1)
				break;
		} /* if ( jcox != 0 ) */
	
		pm=- pp;
		pp=0.;