nec_context.cpp

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

		} /* if ( m_excitation_type == EXCITATION_VOLTAGE_DISC) */
	
		voltage_source_count++;
		source_segment_array.resize(voltage_source_count );
		source_voltage_array.resize(voltage_source_count );
	
		{
			int indx = voltage_source_count-1;
		
			int seg_number = m_geometry->get_segment_number( itmp2, itmp3);
			if (seg_number > m_geometry->segment_length.size())
			{
				nec_exception* nex = new nec_exception("CHECK DATA, PARAMETER SPECIFYING EXCITATION SOURCE SEGMENT [");
				nex->append(seg_number);
				nex->append("] IS TOO LARGE" );
				throw nex;
			}
			source_segment_array[indx] = seg_number;
			
			DEBUG_TRACE("Voltage Source: " << nec_complex( tmp1, tmp2));
			source_voltage_array[indx]= nec_complex( tmp1, tmp2);
			if ( abs( source_voltage_array[indx]) < 1.e-20)
				source_voltage_array[indx] = cplx_10();
		
			iped= itmp4- masym*10;
			impedance_norm_factor= tmp3;
			if ( (iped == 1) && (impedance_norm_factor > 0.0) )
				iped = 2;
			return; /* continue card input loop */
		}	
	} /* if ( (m_excitation_type == 0) || (m_excitation_type == 5) ) */

	nthi= itmp2;
	nphi= itmp3;
	xpr1= tmp1;
	xpr2= tmp2;
	xpr3= tmp3;
	xpr4= tmp4;
	xpr5= tmp5;
	xpr6= tmp6;
	// xpr7= tmp7; Put this in here once we are parsing NEC4 excitation stuff.
	voltage_source_count=0;
	nvqd=0;
	thetis= xpr1;
	phiss= xpr2;
}


/* 5: "tl" cards, network parameters

TL TRANSMISSION LINE 
	I1- PORT 1 TAG #, BLANK/0, USE I2 AS ABSOLUTE
	I2- SEGMENT#, OR ABSOLUTE END 1 SEGMENT, -1=CANCEL NETS/LINES
	I3- AS I1 FOR PORT 2
	I4- AS I2 FOR PORT 2
	F1- LINE Zo, -=CROSSED LINE
	F2- LINE LENGTH METERS, BLANK=STRAIGHT LINE P1 TO P2
	F3- REAL SHUNT ADM., END 1 MHOS
	F4- IMAG SHUNT ADM., END 1
	F5- REAL SHUNT ADM., END 2
	F6- IMAG SHUNT ADM., END 2

*/
void nec_context::tl_card(int itmp1, int itmp2, int itmp3, int itmp4,
	nec_float characteristic_impedance, nec_float tmp2, nec_float tmp3, nec_float tmp4, nec_float tmp5, nec_float tmp6)
{
	if ( iflow != 6)
	{
		network_count=0;
		ntsol=0;
		iflow=6;
	
		if ( processing_state > 3)
			processing_state=3;
	
		if ( itmp2 == -1 )
			return; /* continue card input loop */
	}

	/* Re-allocate network buffers */
	network_count++;
	ntyp.resize(network_count);
	iseg1.resize(network_count);
	iseg2.resize(network_count);
	x11r.resize(network_count);
	x11i.resize(network_count);
	x12r.resize(network_count);
	x12i.resize(network_count);
	x22r.resize(network_count);
	x22i.resize(network_count);

	int idx = network_count-1;
	ntyp[idx] = 2; // TL card

	iseg1[idx]= m_geometry->get_segment_number( itmp1, itmp2);
	iseg2[idx]= m_geometry->get_segment_number( itmp3, itmp4);
	x11r[idx]= characteristic_impedance;
	x11i[idx]= tmp2;
	x12r[idx]= tmp3;
	x12i[idx]= tmp4;
	x22r[idx]= tmp5;
	x22i[idx]= tmp6;

	if (characteristic_impedance <= 0.0)
	{
		// Negative characteristic impedance implies a crossed line
		ntyp[idx] = 3;
		x11r[idx] = -characteristic_impedance;
	}
		
}

/* 4:
	NT	NETWORKS
		I1- PORT 1 TAG #, BLANK/0, USE I2 AS ABSOLUTE
		I2- SEGMENT#, OR ABSOLUTE END 1 SEGMENT, -1=CANCEL NETS/LINES
		I3- AS I1 FOR PORT 2
		I4- AS I2 FOR PORT 2
		F1- REAL OF Y(11), MHOS
		F2- IMAG OF Y(11)
		F3- REAL OF Y(12)
		F4- IMAG OF Y(12)
		F5- REAL OF Y(22)
		F6- IMAG OF Y(22)

*/
void nec_context::nt_card(int itmp1, int itmp2, int itmp3, int itmp4, nec_float tmp1, nec_float tmp2, nec_float tmp3, nec_float tmp4, nec_float tmp5, nec_float tmp6)
{
	if ( iflow != 6)
	{
		network_count=0;
		ntsol=0;
		iflow=6;
	
		if ( processing_state > 3)
			processing_state=3;
	
		if ( itmp2 == -1 )
			return; /* continue card input loop */
	}

	/* Re-allocate network buffers */
	network_count++;
	ntyp.resize(network_count);
	iseg1.resize(network_count);
	iseg2.resize(network_count);
	x11r.resize(network_count);
	x11i.resize(network_count);
	x12r.resize(network_count);
	x12i.resize(network_count);
	x22r.resize(network_count);
	x22i.resize(network_count);

	int idx = network_count-1;
	ntyp[idx]=1; // NT card

	iseg1[idx]= m_geometry->get_segment_number( itmp1, itmp2);
	iseg2[idx]= m_geometry->get_segment_number( itmp3, itmp4);
	x11r[idx]= tmp1;
	x11i[idx]= tmp2;
	x12r[idx]= tmp3;
	x12i[idx]= tmp4;
	x22r[idx]= tmp5;
	x22i[idx]= tmp6;
}

/* "xq" execute card - calc. including radiated fields

XQ	EXECUTE ACCUMULATED CARD DECK
	itmp1-
		0=NO PATTERN,
		1=XY PATTERN,
		2= YZ PATTERN,
		3=BOTH
	(DO NOT USE FOR RADIAL GND SCREEN OR 2ND GND MEDIUM)

	NOTES: FOR A SINGLE FREQUENCY, XQ, NE, NH, RP CAUSE IMMEDIATE EXECUTION
	FOR MULTIPLE FREQS, ONLY XQ, RP CAUSE EXECUTION
*/
void nec_context::xq_card(int itmp1)
{
	DEBUG_TRACE("xq_card(" << itmp1 << ")");
	DEBUG_TRACE("iflow =" << iflow);
	if ( 	((iflow == 10) && (itmp1 == 0)) ||
		((nfrq  ==  1) && (itmp1 == 0) && (iflow > 7)) )
		return; /* continue card input loop */

	if ( itmp1 == 0)
	{
		if ( iflow > 7)
			iflow=11;
		else
			iflow=7;
	}
	else
	{
		ifar=0;
		rfld=0.;
		ipd=0;
		iavp=0;
		m_rp_normalization=0;
		m_rp_output_format=0;
		nth=91;
		nph=1;
		thets=0.0;
		phis=0.0;
		dth=1.0;
		dph=0.0;
	
		if ( itmp1 == 2)
			phis=90.0;
	
		if ( itmp1 == 3)
		{
			nph=2;
			dph=90.0;
		}
	} /* if ( itmp1 == 0) */
	
	simulate(true);
}

/* "gd" card, ground representation */
void nec_context::gd_card(nec_float tmp1, nec_float tmp2, nec_float tmp3, nec_float tmp4)
{
	DEBUG_TRACE("gd_card(" << tmp1 << ")");
	ground.setup_cliff(tmp1, tmp2, tmp3, tmp4);
	iflow=9;
}

/*! \brief Standard radiation pattern parameters 

\param calc_mode
\param n_theta
\param n_phi
\param output_format The output format (0 major axis, minor axis and total gain printed. 
		1 vertical, horizontal ant total gain printed.)
\param normalization Controls the type of normalization of the radiation pattern
N = 0 no normalized gain. 
= 1 major axis gain normalized. 
= 2 minor axis gain normalized. 
= 3 vertical axis gain normalized. 
= 4 horizontal axis gain normalized. 
= 5 total gain normalized.

\param D Selects either power gain or directive gain for both standard printing and normalization. If the structure excitation is an incident plane wave, the quantities printed under the heading "gain" will actually be the scattering cross section (a/lambda 2 ) and will not be affected by the value of d. The column heading for the output will still read "power" or "directive gain," however. 
<ul>
<li> D = 0 power gain. 
<li> D = 1 directive gain.
</ul>

\param A - Requests calculation of average power gain over the region covered by field points. 
<ul>
<li>A = 0 no averaging. 
<li>A = 1 average gain computed. 
<li>A = 2 average gain computed, printing of gain at the field points used for averaging is suppressed. If NTH or NPH is equal to one, average gain will not be computed for any value of A since the area of the region covered by field points vanishes.
</ul>

\param theta0 - Initial theta angle in degrees (initial z coordinate in meters if I1 = 1).

\param phi0 - Initial phi angle in degrees.

\param delta_theta - Increment for theta in degrees (increment for z in meters if I1 = 1).

\param delta_phi - Increment for phi in degrees.

\param radial_distance - Radial distance (R) of field point from the origin in meters. radial_distance is optional. If it is zero, the radiated electric field will have the factor exp(-jkR)/R omitted. If a value of R is specified, it should represent a point in the far-field region since near components of the field cannot be obtained with an RP card. (If I1 = 1, then radial_distance represents the cylindrical coordinate phi in meters and is not optional. It must be greater than about one wavelength.)

\param gain_norm - Determines the gain normalization factor if normalization has been requested in the normalization parameter. If gain_norm is zero, the gain will be normalized to its maximum value. If gain_norm is not zero, the gain w111 be normalized to the value of gain_norm.
*/
void nec_context::rp_card(int calc_mode, int n_theta, int n_phi,
	int output_format, int normalization, int D, int A,
	nec_float theta0, nec_float phi0,
	nec_float delta_theta, nec_float delta_phi,
	nec_float radial_distance, nec_float gain_norm)
{
	DEBUG_TRACE("rp_card(" << calc_mode << ")");
	ifar= calc_mode;
	nth = n_theta;
	nph = n_phi;

	if ( nth == 0)
		nth=1;
	if ( nph == 0)
		nph=1;

	m_rp_output_format = output_format;
	m_rp_normalization = normalization;
	ipd = D;
	iavp = A;
	
	DEBUG_TRACE(" xnda = (" << m_rp_output_format << m_rp_normalization << ipd << iavp << ")");
	if ( m_rp_output_format != 0)
		m_rp_output_format=1;
	if ( ipd != 0)
		ipd=1;
		
	// sanity check. Not point normalizing if there are too few data points to normalize on
	if ( (nth < 2) || (nph < 2) || (ifar == 1) )
		iavp=0;

	thets = theta0;
	phis = phi0;
	dth = delta_theta;
	dph = delta_phi;
	rfld = radial_distance;
	gnor = gain_norm;
	iflow=10;
	
	simulate(true);
}

	/* "pt" card, print control for current */
void nec_context::pt_card(int itmp1, int itmp2, int itmp3, int itmp4)
{
	iptflg= itmp1;
	iptag= itmp2;
	iptagf= itmp3;
	iptagt= itmp4;

	if ( (itmp3 == 0) && (iptflg != -1) )
		iptflg=-2;
	if ( itmp4 == 0)
		iptagt= iptagf;
		
}


	/* "pq" card, print control for charge */
void nec_context::pq_card(int itmp1, int itmp2, int itmp3, int itmp4)
{
	iptflq= itmp1;
	iptaq= itmp2;
	iptaqf= itmp3;
	iptaqt= itmp4;

	if ( (itmp3 == 0) && (iptflq != -1) )
		iptflq=-2;
	if ( itmp4 == 0)
		iptaqt= iptaqf;
}



/* "kh" card, matrix integration limit */
void nec_context::kh_card(nec_float tmp1)
{
	rkh = tmp1;
	if ( processing_state > 2)
		processing_state=2;
	iflow=1;
}

void nec_context::ne_card(int itmp1, int itmp2, int itmp3, int itmp4, nec_float tmp1, nec_float tmp2, nec_float tmp3, nec_float tmp4, nec_float tmp5, nec_float tmp6)
{
	ne_nh_card(0, itmp1, itmp2, itmp3, itmp4, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6);
}
void nec_context::nh_card(int itmp1, int itmp2, int itmp3, int itmp4, nec_float tmp1, nec_float tmp2, nec_float tmp3, nec_float tmp4, nec_float tmp5, nec_float tmp6)
{
	ne_nh_card(1, itmp1, itmp2, itmp3, itmp4, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6);
}

/* \brief Near field calculation parameters 
*/
void nec_context::ne_nh_card(int in_nfeh, int itmp1, int itmp2, int itmp3, int itmp4, nec_float tmp1, nec_float tmp2, nec_float tmp3, nec_float tmp4, nec_float tmp5, nec_float tmp6)
{
	nfeh = in_nfeh;

	if ( (iflow == 8) && (nfrq != 1) )
	{
		m_output.endl(2);
		m_output.line("WHEN MULTIPLE FREQUENCIES ARE REQUESTED, "
			"ONLY ONE NEAR FIELD CARD CAN BE USED -");
		m_output.line("  LAST CARD READ WILL BE USED" );
	}

	m_near= itmp1;
	nrx= itmp2;
	nry= itmp3;
	nrz= itmp4;
	xnr= tmp1;
	ynr= tmp2;
	znr= tmp3;
	dxnr= tmp4;
	dynr= tmp5;
	dznr= tmp6;
	iflow=8;

	if ( nfrq == 1)
		simulate();
}


/* "ek" card,  extended thin wire kernel option */
void nec_context::set_extended_thin_wire_kernel(bool ek_flag)
{
	m_use_exk = ek_flag;
	
	if ( processing_state > 2)
		processing_state=2;
	iflow=1;	
}


/* "cp" card, maximum coupling between antennas */
void nec_context::cp_card(int itmp1, int itmp2, int itmp3, int itmp4)
{
	if ( iflow != 2)
	{
		ncoup=0;
		nctag.resize(0);
		ncseg.resize(0);
		y11a.resize(0);
		y12a.resize(0);
	}

	icoup=0;
	iflow=2;

	if ( itmp2 == 0)
		return; /* continue card input loop */

	ncoup++;
	nctag.resize(ncoup);
	ncseg.resize(ncoup);
	nctag[ncoup-1]= itmp1;
	ncseg[ncoup-1]= itmp2;

	if ( itmp4 == 0)
		return; /* continue card input loop */

	ncoup++;
	nctag.resize(ncoup);
	ncseg.resize(ncoup);
	nctag[ncoup-1]= itmp3;
	ncseg[ncoup-1]= itmp4;
}


/* "pl" card, plot flags 
	throws int on error.
*/
void nec_context::pl_card(const char* ploutput_filename, int itmp1, int itmp2, int itmp3, int itmp4)
{
	std::string fname(ploutput_filename);
	plot_card = c_plot_card(itmp1,itmp2,itmp3,itmp4, fname);
}



/*!	\brief Start a simulation

	This function will trigger a calculation. In the traditional NEC
	world, This signals the end of the main input section and the
	beginning of the frequency do loop.
	
	\param far_field_flag is true if last card was XQ or RP
	\warning far_field_flag is should never be specified as true
	because both the xq_card() and rp_card() functions will call
	this function automatically.
*/
void nec_context::simulate(bool far_field_flag)
{
	DEBUG_TRACE("simulate(" << far_field_flag << ")");
	/* Allocate the normalization buffer */
	{
		int mreq1 = 0;
		int mreq2 = 0;
	
		if ( iped )
			mreq1 = 4*nfrq;
		if ( iptflg >= 2 )
			mreq2 = nthi*nphi;
	
		int newsize = std::max(mreq1,mreq2);
		fnorm.resize(newsize);
	}

	/* igox is a state variable that is used to change from
		one processing mode to another. The processing to be 
		performed are as follows:
			
		1: Memory allocation & Initialization
		2: Structure segment loading
		3: Excitation set up (right hand side, -e inc.)
		4: ?
		5: Near field calculation
		6: standard far field calculation
	*/
	int igox;
	int mhz = 0;
	
	if ( (far_field_flag == true)
		&& (processing_state == 5) )
		igox = 6;
	else
		igox = processing_state;
	
	
	try
	{
	int iresrv = 0;
	bool in_freq_loop = false;
	
	do
	{	
		switch( igox )
		{
		case 1: /* Memory allocation for primary interacton matrix. */

			if (false == in_freq_loop)
			{
				iresrv = m_geometry->n_plus_2m * (m_geometry->np+2*m_geometry->mp);
				cm.resize(iresrv);
			
				/* Memory allocation for symmetry array */
				nop = neq/npeq;