c_geometry.cpp

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

					nex->append(i+1);
					nex->append("LIES IN PLANE OF SYMMETRY");
					throw nex;
				}
			
				px[nx]= px[i];
				py[nx]=- py[i];
				pz[nx]= pz[i];
				t1x[nx]= t1x[i];
				t1y[nx]=- t1y[i];
				t1z[nx]= t1z[i];
				t2x[nx]= t2x[i];
				t2y[nx]=- t2y[i];
				t2z[nx]= t2z[i];
				psalp[nx]=- psalp[i];
				pbi[nx]= pbi[i];
			
			} /* for( i = m2; i <= m; i++ ) */
		
			m= m*2;
		} /* if ( m >= m2) */
	
		} /* if ( iy != 0) */
	
		/* reflect along x axis */
		if ( ix == 0 )
			return;
	
		if ( n > 0 )
		{
		/* Reallocate tags buffer */
		segment_tags.resize(2*n + m);
		// segment_tags.resize((2*n+m));/*????*/
	
		/* Reallocate wire buffers */
		int new_size = 2*n;
		x.resize(new_size);
		y.resize(new_size);
		z.resize(new_size);
		x2.resize(new_size);
		y2.resize(new_size);
		z2.resize(new_size);
		segment_radius.resize(new_size);
	
		for( i = 0; i < n; i++ )
		{
			nx= i+ n;
			e1= x[i];
			e2= x2[i];
		
			if ( (fabs(e1)+fabs(e2) <= 1.0e-5) || (e1*e2 < -1.0e-6) )
			{
				nec_exception* nex = new nec_exception("GEOMETRY DATA ERROR--SEGMENT ");
				nex->append(i+1);
				nex->append("LIES IN PLANE OF SYMMETRY");
				throw nex;
			}
		
			x[nx]=- e1;
			y[nx]= y[i];
			z[nx]= z[i];
			x2[nx]=- e2;
			y2[nx]= y2[i];
			z2[nx]= z2[i];
			itagi= segment_tags[i];
		
			if ( itagi == 0)
			segment_tags[nx]=0;
			if ( itagi != 0)
			segment_tags[nx]= itagi+ iti;
		
			segment_radius[nx]= segment_radius[i];
		}
	
		n= n*2;
	
		} /* if ( n > 0) */
	
		if ( m == 0 )
			return;
	
		/* Reallocate patch buffers */
		int new_size = 2*m;
		px.resize(new_size);
		py.resize(new_size);
		pz.resize(new_size);
		t1x.resize(new_size);
		t1y.resize(new_size);
		t1z.resize(new_size);
		t2x.resize(new_size);
		t2y.resize(new_size);
		t2z.resize(new_size);
		pbi.resize(new_size);
		psalp.resize(new_size);
	
		for( i = 0; i < m; i++ )
		{
			nx= i+m;
			if ( fabs( px[i]) <= 1.0e-10)
			{
				nec_exception* nex = new nec_exception("GEOMETRY DATA ERROR--PATCH ");
				nex->append(i+1);
				nex->append("LIES IN PLANE OF SYMMETRY");
				throw nex;
			}
		
			px[nx]=- px[i];
			py[nx]= py[i];
			pz[nx]= pz[i];
			t1x[nx]=- t1x[i];
			t1y[nx]= t1y[i];
			t1z[nx]= t1z[i];
			t2x[nx]=- t2x[i];
			t2y[nx]= t2y[i];
			t2z[nx]= t2z[i];
			psalp[nx]=- psalp[i];
			pbi[nx]= pbi[i];
		}
	
		m= m*2;
		return;
	
	} /* if ( ix >= 0) */
	
	/* reproduce structure with rotation to form cylindrical structure */
	fnop= (nec_float)nop;
	m_ipsym=-1;
	sam=two_pi() / fnop;
	cs= cos( sam);
	ss= sin( sam);
	
	if ( n > 0)
	{
		n *= nop;
		nx= np;
	
		/* Reallocate tags buffer */
		segment_tags.resize(n + m);
		//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);
	
		for( i = nx; i < n; i++ )
		{
			k= i- np;
			xk= x[k];
			yk= y[k];
			x[i]= xk* cs- yk* ss;
			y[i]= xk* ss+ yk* cs;
			z[i]= z[k];
			xk= x2[k];
			yk= y2[k];
			x2[i]= xk* cs- yk* ss;
			y2[i]= xk* ss+ yk* cs;
			z2[i]= z2[k];
			segment_radius[i]= segment_radius[k];
			itagi= segment_tags[k];
		
			if ( itagi == 0)
				segment_tags[i]=0;
			if ( itagi != 0)
				segment_tags[i]= itagi+ iti;
		}
	
	} /* if ( n >= n2) */
	
	if ( m == 0 )
		return;
	
	m *= nop;
	nx= mp;
	
	/* 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);
	
	for( i = nx; i < m; i++ )
	{
		k = i-mp;
		xk= px[k];
		yk= py[k];
		px[i]= xk* cs- yk* ss;
		py[i]= xk* ss+ yk* cs;
		pz[i]= pz[k];
		xk= t1x[k];
		yk= t1y[k];
		t1x[i]= xk* cs- yk* ss;
		t1y[i]= xk* ss+ yk* cs;
		t1z[i]= t1z[k];
		xk= t2x[k];
		yk= t2y[k];
		t2x[i]= xk* cs- yk* ss;
		t2y[i]= xk* ss+ yk* cs;
		t2z[i]= t2z[k];
		psalp[i]= psalp[k];
		pbi[i]= pbi[k];
	
	} /* for( i = nx; i < m; i++ ) */
}
	
/*-----------------------------------------------------------------------*/

/* connect sets up segment connection data in arrays icon1 and */
/* icon2 by searching for segment ends that are in contact. */
void c_geometry::connect_segments( int ignd )
{
	nscon= -1;
	maxcon = 1;
	
	if (n <= 1)
	{
		throw new nec_exception("GEOMETRY HAS ONE OR FEWER SEGMENTS. Please send bug report. This causes an error that we're trying to fix.");
	}
	
	if ( ignd != 0)
	{
		m_output->nec_printf( "\n\n     GROUND PLANE SPECIFIED." );
	
		if ( ignd > 0)
			m_output->nec_printf(
				"\n     WHERE WIRE ENDS TOUCH GROUND, CURRENT WILL"
				" BE INTERPOLATED TO IMAGE IN GROUND PLANE.\n" );
	
		if ( m_ipsym == 2)
		{
			np=2* np;
			mp=2* mp;
		}
	
		if ( abs( m_ipsym) > 2 )
		{
			np= n;
			mp= m;
		}
	
		if ( np > n)
		{
			throw new nec_exception("ERROR: NP > N IN c_geometry::connect_segments()" );
		}
	
		if ( (np == n) && (mp == m) )
			m_ipsym=0;
	} /* if ( ignd != 0) */
	
	if ( n != 0)
	{
		/* Allocate memory to connections */
		icon1.resize((n+m));
		icon2.resize((n+m));
	
		for (int i = 0; i < n; i++ )
		{
			int iz = i+1;
		
			nec_float zi1 = z[i];
			nec_float zi2 = z2[i];
			
			nec_3vector v1(x[i], y[i], z[i]);
			nec_3vector v2(x2[i], y2[i], z2[i]);
			nec_float slen = norm(v2 - v1) * SMIN;		
			
			/* determine connection data for end 1 of segment. */
			bool segment_on_ground = false;
			if ( ignd > 0)
			{
				if ( zi1 <= -slen)
				{
					nec_exception* nex = new nec_exception("GEOMETRY DATA ERROR--SEGMENT ");
					nex->append(iz);
					nex->append("EXTENDS BELOW GROUND");
					throw nex;
				}
			
				if ( zi1 <= slen)
				{
					icon1[i]= iz;
					z[i]=0.;
					segment_on_ground = true;	
				} /* if ( zi1 <= slen) */
			} /* if ( ignd > 0) */
		
			if ( false == segment_on_ground )
			{
				int ic= i;
				nec_float sep=0.0;
				for (int j = 1; j < n; j++)
				{
					ic++;
					if ( ic >= n)
						ic=0;
				
					nec_3vector vic(x[ic], y[ic], z[ic]);
					sep = normL1(v1 - vic);
					if ( sep <= slen)
					{
						icon1[i]= -(ic+1);
						break;
					}
				
					nec_3vector v2ic(x2[ic], y2[ic], z2[ic]);
					sep = normL1(v1 - v2ic);
					if ( sep <= slen)
					{
						icon1[i]= (ic+1);
						break;
					}
				
				} /* for( j = 1; j < n; j++) */
			
				if ( ((iz > 0) || (icon1[i] <= PCHCON)) && (sep > slen) )
					icon1[i]=0;
			
			} /* if ( ! jump ) */
		
			/* determine connection data for end 2 of segment. */
			if ( (ignd > 0) || segment_on_ground )
			{
				if ( zi2 <= -slen)
				{
					nec_exception* nex = new nec_exception("GEOMETRY DATA ERROR--SEGMENT ");
					nex->append(iz);
					nex->append("EXTENDS BELOW GROUND");
					throw nex;
				}
			
				if ( zi2 <= slen)
				{
					if ( icon1[i] == iz )
					{
						nec_exception* nex = new nec_exception("GEOMETRY DATA ERROR--SEGMENT ");
						nex->append(iz);
						nex->append("LIES IN GROUND PLANE");
						throw nex;
					}
				
					icon2[i]= iz;
					z2[i]=0.;
					continue;
				} /* if ( zi2 <= slen) */	
			} /* if ( ignd > 0) */
		
			int ic= i;
			nec_float sep=0.0;
			for (int j = 1; j < n; j++ )
			{
				ic++;
				if ( ic >= n)
					ic=0;
			
				nec_3vector vic(x[ic], y[ic], z[ic]);
				sep = normL1(v2 - vic);
				if (sep <= slen)
				{
					icon2[i]= (ic+1);
					break;
				}
			
				nec_3vector v2ic(x2[ic], y2[ic], z2[ic]);
				sep = normL1(v2 - v2ic);
				if (sep <= slen)
				{
					icon2[i]= -(ic+1);
					break;
				}
			} /* for( j = 1; j < n; j++ ) */
		
			if ( ((iz > 0) || (icon2[i] <= PCHCON)) && (sep > slen) )
				icon2[i]=0;
		
		} /* for( i = 0; i < n; i++ ) */
		

		/* find wire-surface connections for new patches */
		for (int ix=0; ix <m; ix++)
		{
//			DEBUG_TRACE("i: " << ix+1 << " ix: " << ix << " m: " << m);		
			nec_3vector vs(px[ix], py[ix], pz[ix]);
		
			for (int iseg = 0; iseg < n; iseg++ )
			{
				nec_3vector v1(x[iseg], y[iseg], z[iseg]);
				nec_3vector v2(x2[iseg], y2[iseg], z2[iseg]);
			
				/* for first end of segment */
				nec_float slen = normL1(v2 - v1) * SMIN;
				
				nec_float sep = normL1(v1 - vs);
				/* connection - divide patch into 4 patches at present array loc. */
				if ( sep <= slen)
				{
					icon1[iseg]=PCHCON + ix + 1;
					divide_patch(ix + 1);
					break;
				}
			
				sep = normL1(v2 - vs);
				
				if ( sep <= slen)
				{
					icon2[iseg]=PCHCON+ ix + 1;
					divide_patch(ix + 1);
					break;
				}
			}
		}	
	} /* if ( n != 0) */
	
	m_output->nec_printf( "\n\n"
		"     TOTAL SEGMENTS USED: %d   SEGMENTS IN A"
		" SYMMETRIC CELL: %d   SYMMETRY FLAG: %d",
		n, np, m_ipsym );
	
	if ( m > 0)
		m_output->nec_printf(	"\n"
		"       TOTAL PATCHES USED: %d   PATCHES"
		" IN A SYMMETRIC CELL: %d",  m, mp );
	
	
	if (0 == np + mp)
		throw new nec_exception("connect_segments Geometry has zero wires and zero patches.");
	
	int symmetry = (n+m)/(np+mp); /* was iseg */
	if ( symmetry != 1)
	{
		/*** may be error condition?? ***/
		if ( m_ipsym == 0 )
		{
			nec_error_mode nem(*m_output);
			m_output->endl();
			m_output->line("ERROR: IPSYM=0 IN connect_segments()" );
			throw new nec_exception("ERROR: IPSYM=0 IN connect_segments()");
		}
	
		if ( m_ipsym < 0 )
			m_output->nec_printf(
				"\n  STRUCTURE HAS %d FOLD ROTATIONAL SYMMETRY\n", symmetry );
		else
		{
			int sym_planes = symmetry/2;
			if ( symmetry == 8)
				sym_planes=3;
			m_output->nec_printf(
				"\n  STRUCTURE HAS %d PLANES OF SYMMETRY\n", sym_planes );
		} /* if ( m_ipsym < 0 ) */
	
	} /* if ( symmetry == 1) */
	
	if ( n == 0)
		return;
	
	/* Allocate to connection buffers */
	jco.resize(maxcon);
	
	/*
		Adjust connected seg. ends to exactly coincide.  print junctions
		of 3 or more seg.  also find old seg. connecting to new seg.
	*/
	int junction_counter = 0; // used just to print the junction number out if there are 3 or more segments
	bool header_printed = false; // Have we printed the header
	for (int j = 0; j < n; j++ )
	{
		int jx = j+1;
		int iend = -1;
		int jend = -1;
		int ix= icon1[j];
		int ic=1;
		jco[0]= -jx;
		nec_float xa = x[j];
		nec_float ya = y[j];
		nec_float za = z[j];
	
		while( true )
		{
		if ( (ix != 0) && (ix != (j+1)) && (ix <= PCHCON) )
		{
		int nsflg=0;
	
		bool jump = false;
		
		do
		{
			if ( ix == 0 )
			{
				nec_exception* nex = new nec_exception("CONNECT - SEGMENT CONNECTION ERROR FOR SEGMENT: ");
				nex->append(ix);
				throw nex;
			}
		
			if ( ix < 0 )
				ix= -ix;
			else
				jend= -jend;
		
			jump = false;
		
			if ( ix == jx )
				break;
		
			if ( ix < jx )
			{
				jump = true;
				break;
			}
		
			/* Record max. no. of connections */
			ic++;
			if ( ic >= maxcon )
			{
				maxcon = ic+1;
				jco.resize(maxcon);
			}
			jco[ic-1]= ix* jend;
		
			if ( ix > 0)
				nsflg=1;
		
			int ixx = ix-1;
			if ( jend != 1)
			{
				xa= xa+ x[ixx]; // dies here if n == 1. ix is totally fried.
				ya= ya+ y[ixx];
				za= za+ z[ixx];
				ix= icon1[ixx];
				continue;
			}
		
			xa= xa+ x2[ixx];
			ya= ya+ y2[ixx];
			za= za+ z2[ixx];
			ix= icon2[ixx];
		} /* do */
		while( ix != 0 );
	
		if ( jump && (iend == 1) )
			break;
		else
			if ( jump )
			{
				iend=1;
				jend=1;
				ix= icon2[j];
				ic=1;
				jco[0]= jx;
				xa= x2[j];
				ya= y2[j];
				za= z2[j];
				continue;
			}
	
		nec_float ic_f = (nec_float)ic;
		xa = xa / ic_f;
		ya = ya / ic_f;
		za = za / ic_f;
	
		for (int i = 0; i < ic; i++ )
		{
			ix= jco[i];
			if ( ix <= 0)
			{
				ix=- ix;
				int ixx = ix-1;
				x[ixx]= xa;
				y[ixx]= ya;
				z[ixx]= za;
				continue;
			}
		
			int ixx = ix-1;
			x2[ixx]= xa;
			y2[ixx]= ya;
			z2[ixx]= za;	
		} /* for( i = 0; i < ic; i++ ) */
	
		if ( ic >= 3)
		{
			if ( false == header_printed )
			{
				m_output->nec_printf( "\n\n"
					"    ---------- MULTIPLE WIRE JUNCTIONS ----------\n"
					"    JUNCTION  SEGMENTS (- FOR END 1, + FOR END 2)" );
				header_printed = true;
			}
		
			junction_counter++;
			m_output->nec_printf( "\n   %5d      ", junction_counter );
		
			for (int i = 1; i <= ic; i++ )
			{
				m_output->nec_printf( "%5d", jco[i-1] );
				if ( !(i % 20) )
					m_output->nec_printf( "\n              " );
			}
		} /* if ( ic >= 3) */
	
		} /*if ( (ix != 0) && (ix != j) && (ix <= PCHCON) ) */
	
		if ( iend == 1)
			break;
	
		iend=1;
		jend=1;
		ix= icon2[j];
		ic=1;
		jco[0]= jx;
		xa= x2[j];
		ya= y2[j];
		za= z2[j];
	
		} /* while( true ) */
	
	} /* for( j = 0; j < n; j++ ) */