nec_context.cpp

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

/*
	Copyright (C) 2004-2005  Timothy C.A. Molteno
	[e-mail]
	
	This program is free software; you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation; either version 2 of the License, or
	(at your option) any later version.
	
	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.
	
	You should have received a copy of the GNU General Public License
	along with this program; if not, write to the Free Software
	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#include "nec_context.h"
#include "c_geometry.h"
#include "nec_exception.h"
	
static char* hpol[3] = { "LINEAR", "RIGHT", "LEFT" };

nec_context::nec_context()
{
	m_output_fp=NULL;
	m_results.set_stdout(false);
	m_geometry = new c_geometry();
	inc=0;
	isave=0;
	nthic=0;
	nphic=0;
	
	impedance_norm_factor = 0.0; // was zpnorm
	xpr1=0.0;
	xpr2=0.0;
	xpr3=0.0;
	xpr4=0.0;
	xpr5=0.0;
	xpr7=0.0;
	
	// allocate the ground grid
	ggrid.initialize();
}

nec_context::~nec_context()
{
	delete m_geometry;
}

/*! \brief
Initialize everything, called after construction so that we
can tell the geometry object what nec_context object to
point to.
*/
void nec_context::initialize()
{
	DEBUG_TRACE("initialize()");
	nthi=0;
	nphi=0;
	iflow = 1;
	thetis=0.0;
	phiss=0.0;
	
	iptag=0;
	iptagf=0;
	iptagt=0;
	iptaq=0;
	iptaqf=0;
	iptaqt=0;
	
	init_voltage_sources();
	m_geometry->set_context(this);

	imat=0; // this should go away! It appears to be related to some very old code for handling out-of-core situations */
}

/*! \brief Private method to initialize the voltage source buffers
*/
void nec_context::init_voltage_sources()
{
	/* Free vsource buffers */
	ivqd.resize(0);
	iqds.resize(0);
	vqd.resize(0);
	vqds.resize(0);
	source_segment_array.resize(0);
	source_voltage_array.resize(0);

	voltage_source_count=0;
	nvqd=0;
	iped=0;
}

/*! \brief After the geometry has been specified, this function prepares for calculations
*/
void nec_context::calc_prepare()
{
	DEBUG_TRACE("calc_prepare()");
	iflow=1;		
	
	int n_plus_m = m_geometry->n_plus_m;
	/* Allocate some buffers */
	air.resize(n_plus_m);
	aii.resize(n_plus_m);
	bir.resize(n_plus_m);
	bii.resize(n_plus_m);
	cir.resize(n_plus_m);
	cii.resize(n_plus_m);
	
	ip.resize(m_geometry->n_plus_2m);		
	current_vector.resize(m_geometry->n_plus_3m);
	
	/* Matrix parameters */
	if ( imat == 0)
	{
		neq= m_geometry->n_plus_2m;
		neq2=0;
	}

	/* default values for input parameters and flags */
	npeq = m_geometry->np + 2*m_geometry->mp;
	processing_state=1;
	rkh=1.;
	m_use_exk=false;
	m_excitation_type = EXCITATION_VOLTAGE;
	nload=0;
	network_count=0;
	m_near=-1;
	ifar=-1;
	ncoup=0;
	icoup=0;
	freq_mhz= CVEL;
	ground.default_values();
	nfrq=1;
	iptflg=-2;
	iptflq=-1;
	iped=0;
}


/*!\brief Benchmark the libnecpp engine. A score of 100 is roughly an Athlon XP 1800. */
nec_float nec_context::benchmark()
{
	nec_float start_timer, stop_timer;
	
	secnds( &start_timer );
	for (int i=0; i<2; i++)
	{
		{
			/*
			CMEXAMPLE 2. CENTER FED LINEAR ANTENNA. 
			CM           CURRENT SLOPE DISCONTINUITY SOURCE. 
			CM           1. THIN PERFECTLY CONDUCTING WIRE 
			CE           2. THIN ALUMINUM WIRE 
			GW 0 8 0. 0. -.25 0. 0. .25 .00001 
			GE 
			FR 0 3 0 0 200. 50. 
			EX 5 0 5 1 1. 0. 50. 
			XQ 
			LD 5 0 0 0 3.720E+07 
			FR 0 1 0 0 300. 
			EX 5 0 5 0 1. 
			GN 1
			XQ 
			EN
			*/

			nec_context nec;
			nec.set_results_stdout(false);
			nec.set_gain_only(true);
			nec.initialize();
			
			c_geometry* geo = nec.get_geometry();
			geo->wire(0,8, 0.0, 0.0, -0.25,
					0.0, 0.0, 0.25,
					0.00001,
					1.0, 1.0);
					
			nec.geometry_complete(0,0);
			
			nec.fr_card(0, 3, 200.0, 50.0);
			nec.ex_card(EXCITATION_VOLTAGE_DISC, 0, 5, 1, 1.0, 0.0, 50.0, 0.0, 0.0, 0.0);
			nec.xq_card(0);
			nec.ld_card(5, 0, 0,0, 3.72e7, 0.0, 0.0);
			nec.fr_card(0, 1, 300.0, 0.0);
			nec.ex_card(EXCITATION_VOLTAGE_DISC, 0, 5, 0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0);
			nec.gn_card(1, 0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0);
			
			// TODO Keep testing here...
			nec.xq_card(0);
	/*	This line kills the version on Wine! Very strange! I think that the gn card is bad
			fixme:msvcrt:MSVCRT_signal (11 (nil)):stub
			err:seh:EXC_DefaultHandling Unhandled exception code c0000005 flags 0 addr 0x404e2380
			Wine failed with return code 1
	*/
		}
		{
			/*
				An example showing how to evaluate for maximum gain.
				
			CM GA - NEC FILE
			CE go blue ! 
			GW 0 36 0 0 0 -0.042 0.008 0.017 0.001
			GW 0 21 -0.042 0.008 0.017 -0.048 0.021 -0.005 0.001
			GW 0 70 -0.048 0.021 -0.005 0.039 0.032 -0.017 0.001
			GW 0 70 -0.048 0.021 -0.005 0.035 0.043 0.014 0.001
			GW 0 50 -0.042 0.008 0.017 0.017 -0.015 0.014 0.001
			GW 0 66 0.017 -0.015 0.014 -0.027 0.04 -0.031 0.001
			GW 0 85 -0.027 0.04 -0.031 0.046 -0.01 0.028 0.001
			GW 0 47 0.046 -0.01 0.028 -0.013 -0.005 0.031 0.001
			GW 0 70 0.017 -0.015 0.014 -0.048 -0.038 -0.04 0.001
			GW 0 77 -0.048 -0.038 -0.04 0.049 -0.045 -0.04 0.001
			GE 0
			GN -1
			LD 5 0 0 0 3.720E+07 
			FR 0 1 0 0 2400
			PT -1
			EX 1 1 1 0 0 0 0 0 0 0 0
			RP 0 1 1 0500 90 90 0 0
			EN
			
			*/
			nec_context nec;
			nec.set_results_stdout(false);
			nec.set_gain_only(false);
			nec.initialize();
			
			c_geometry* geo = nec.get_geometry();
			geo->wire(0, 36, 0, 0, 0, -0.042, 0.008, 0.017, 0.001, 1.0, 1.0);
			geo->wire(0, 21, -0.042, 0.008, 0.017, -0.048, 0.021, -0.005, 0.001, 1.0, 1.0);
			geo->wire(0, 70, -0.048, 0.021, -0.005, 0.039, 0.032, -0.017, 0.001, 1.0, 1.0);
			geo->wire(0, 70, -0.048, 0.021, -0.005, 0.035, 0.043, 0.014, 0.001, 1.0, 1.0);
			geo->wire(0, 50, -0.042, 0.008, 0.017, 0.017, -0.015, 0.014, 0.001, 1.0, 1.0);
			geo->wire(0, 66, 0.017, -0.015, 0.014, -0.027, 0.04, -0.031, 0.001, 1.0, 1.0);
			geo->wire(0, 85, -0.027, 0.04, -0.031, 0.046, -0.01, 0.028, 0.001, 1.0, 1.0);
			geo->wire(0, 47, 0.046, -0.01, 0.028, -0.013, -0.005, 0.031, 0.001, 1.0, 1.0);
			geo->wire(0, 70, 0.017, -0.015, 0.014, -0.048, -0.038, -0.04, 0.001, 1.0, 1.0);
			geo->wire(0, 77, -0.048, -0.038, -0.04, 0.049, -0.045, -0.04, 0.001, 1.0, 1.0);
			nec.geometry_complete(0,0);
			
			nec.gn_card(-1,0,0.0, 0.0, 0.0,0.0, 0.0, 0.0);
			nec.ld_card(5,0,0,0,3.72e7,0.0,0.0);
			nec.pt_card(-1, 0, 0, 0);
			nec.ex_card(EXCITATION_LINEAR, 1, 1, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
			nec.fr_card(0, 2, 2400.0, 100.0);
			nec.rp_card(0, 1, 1, 0,5,0,0, 90.0, 90.0, 0.0, 0.0, 0.0, 0.0);
			double g = 0;
			g = nec.get_maximum_gain();
		}
	}
	/* time the process */
	secnds( &stop_timer );
	stop_timer -= start_timer;
	cout << endl << endl;
	nec_float sec = stop_timer / 1000.0;
	
	nec_float bench = 7.0 / sec;
	
	return bench;
}


/*! \brief Signal the end of a geometry description.

This function prepares for a calculation by calling calc_prepare().
*/
void nec_context::geometry_complete(int card_int_1, int card_int_2)
{
	DEBUG_TRACE("geometry_complete()");
	m_geometry->geometry_complete(this, card_int_1, card_int_2);
	calc_prepare();
}

/*! Add a wire to the geometry,

All co-ordinates are in meters.

	\param tag_id The tag ID.
	\param segment_count The number of segments.
	\param xw1 The x coordinate of the wire starting point.
	\param yw1 The y coordinate of the wire starting point.
	\param zw1 The z coordinate of the wire starting point.
	\param xw2 The x coordinate of the wire ending point.
	\param yw2 The y coordinate of the wire ending point.
	\param zw2 The z coordinate of the wire ending point.
	\param rad The wire radius (meters)
	\param rdel For tapered wires, the. Otherwise set to 1.0
	\param rrad For tapered wires, the. Otherwise set to 1.0
*/
void nec_context::wire(int tag_id, int segment_count,
	nec_float xw1, nec_float yw1, nec_float zw1,
	nec_float xw2, nec_float yw2, nec_float zw2,
	nec_float rad, nec_float rdel, nec_float rrad)
{
	m_geometry->wire(tag_id, segment_count, xw1, yw1, zw1, xw2, yw2, zw2, rad, rdel, rrad);
}

/*! Add an arc to the geometry,

All co-ordinates are in meters.

	\param tag_id The tag ID.
	\param segment_count The number of segments.
	\param rada The radius.
	\param ang1 The angle of the arc starting point.
	\param ang2 The angle of the arc end point.
	\param rad The wire radius.
*/
void nec_context::arc( int tag_id, int segment_count, nec_float rada,
		nec_float ang1, nec_float ang2, nec_float rad )
{
	m_geometry->arc(tag_id, segment_count, rada, ang1, ang2, rad);
}
		

/*! \brief Add an helix to the geometry,

\remark The helix is a versatile m_geometry->element. For example, to generate a spiral printed circuit antenna, use a helix of zero height.

All co-ordinates are in meters.

	\param tag_id The tag ID.
	\param segment_count The number of segments.
	\param s The turn spacing.
	\param h1 The total length of the helix (negative for a left-handed helix).
	\param a1 x-start radius.
	\param b1 y-start radius.
	\param a2 x-end radius.
	\param b2 y-end radius.
	\param rad The wire radius.
*/
void nec_context::helix(int tag_id, int segment_count, nec_float s, nec_float hl, nec_float a1, nec_float b1,
		nec_float a2, nec_float b2, nec_float rad)
{
	m_geometry->helix(s, hl, a1, b1, a2, b2, rad, segment_count, tag_id);
}

/* "fr" card, frequency parameters

FREQUENCY
I1- O= LINEAR STEP, 1=MULTIPLICATIVE
I2- NO. STEPS, BLANK=1
I3- BLANK -- not used in this function
I4- BLANK -- not used in this function
F1- FREQUENCY OR START FREQUENCY
F2- FREQ INCREMENT, ADD OR MULTIPLY
*/
void nec_context::fr_card(int in_ifrq, int in_nfrq, nec_float in_freq_mhz, nec_float in_del_freq)
{
	DEBUG_TRACE("fr_card()");
	ifrq = in_ifrq;
	nfrq = in_nfrq;
	if ( nfrq == 0)
		nfrq=1;
		
	freq_mhz = in_freq_mhz;
	delfrq = in_del_freq;
	if ( iped == 1)
		impedance_norm_factor = 0.0;
		
	processing_state = 1;
	iflow = 1;
}



void nec_context::ld_card(int itmp1, int itmp2, int itmp3, int itmp4, nec_float tmp1, nec_float tmp2, nec_float tmp3)
{
	DEBUG_TRACE("ld_card()");
	if ( iflow != 3 )
	{
		iflow=3;
		/* Free loading buffers */
		nload=0;
		ldtyp.resize(0);
		ldtag.resize(0);
		ldtagf.resize(0);
		ldtagt.resize(0);
		zlr.resize(0);
		zli.resize(0);
		zlc.resize(0);
	
		if ( processing_state > 2 )
			processing_state=2;
		if ( itmp1 == -1 )
			return; // continue card input loop
	}

	/* Reallocate loading buffers */
	nload++;
	ldtyp.resize(nload);
	ldtag.resize(nload);
	ldtagf.resize(nload);
	ldtagt.resize(nload);
	zlr.resize(nload);
	zli.resize(nload);
	zlc.resize(nload);

	int idx = nload-1;
	ldtyp[idx]= itmp1;
	ldtag[idx]= itmp2;
	if ( itmp4 == 0)
		itmp4= itmp3;
	ldtagf[idx]= itmp3;
	ldtagt[idx]= itmp4; 

	if ( itmp4 < itmp3 )
	{
		nec_stop("DATA FAULT ON LOADING CARD No: %d: ITAG "
			"STEP1: %d IS GREATER THAN ITAG STEP2: %d",
			nload, itmp3, itmp4 );
	}

	zlr[idx]= tmp1;
	zli[idx]= tmp2;
	zlc[idx]= tmp3;
}



/* "gn" card, ground parameters under the antenna
GN		NEAR GROUND, GROUND SCREEN, ADDED GROUND
	I1- -1=SET FREE SPACE (A), 0=REFL COEFF, 1=IDEAL (B), 2-SOMMERFIELD
	I2- (A) BLANK), NO WIRES IN GND SCREEN (C), 0= NO WIRES (D)
	I3- BLANK
	I4- BLANK
	F1- (A,B) BLANK, DIELECTRIC OF NEAR GROUND
	F2- (A,B) BLANK, CONDUCTIVITY OF NEAR GROUND
	F3- (A,B) BLANK, (C) RADIUS OF SCREEN, (D) DIELECTRIC 2ND MEDIUM
	F4- (A,B) BLANK, (C) RADII SCREEN WIRES, (D) CONDUCT. 2ND MEDIUM
	F5- (A,B) BLANK, (C) BLANK, (D) DIST TO 2ND MEDIUM, SEE RP
	F6- (A,B) BLANK, (C) BLANK, (D) HEIGHT 2ND MEDIUM (AS IN GD)
*/
void nec_context::gn_card(int ground_type, int rad_wire_count, nec_float tmp1, nec_float tmp2, nec_float tmp3, nec_float tmp4, nec_float tmp5, nec_float tmp6)
{
	DEBUG_TRACE("gn_card(" << ground_type << ")");
	ground.parse_gn(ground_type, rad_wire_count, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6);
	
	iflow=4;

	if ( processing_state > 2)
		processing_state=2;
}



/*
		EX	EXCITE STRUCTURE, LAST ENCOUNTERED=USED
			I1- 0=E VOLTAGE (A), 1=LINEAR WAVE (B), 2= R CIRC WAVE (B)
			3=L CIRC WAVE (B), 4= CURRENT (C), 5= VOLTAGE DISC. (A)
			I2- (A) SOURCE TAG#, (B) # TH ANGLS, (C) BLANK
			I3- (A) SOURCE SEG#, (B) # PH ANGLS, (C) BLANK
			I4- (A) XX= ADMIT.,IMPED. PRINT, X=0 NO/1 DO, (BC), 1= ADM. PRINT
			F1- (A) EREAL, (B) TH ANGL, (C) X OF SOURCE
			F2- (A) EIMAG, (B) PH ANGL, (C) Y OF SOURCE
			F3- (A) NORM FOR I4, (B) ET ANGL, Z OF SOURCE
			F4- (A) BLANK, (B) TH INC, (C) ALPHA ANGLE FROM XY
			F5- (A) BLANK, (B) PH INC, (C) BETA ANGLE FROM X
			F6- (A) BLANK, (B) MIN/MAJ AXIS, PRODUCT AMPS X LENGTH
			
			// NOT YET DONE... F7- (A) BLANK, (B) INCIDENT AMPLITUDE (Volts/m)
*/
void nec_context::ex_card(enum excitation_type 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)
{
	DEBUG_TRACE("ex_card(" << itmp1
		<< "," << itmp2 
		<< "," << itmp3 
		<< "," << itmp4 
		<< "," << tmp1 
		<< "," << tmp2 
		<< "," << tmp3 
		<< "," << tmp4 
		<< "," << tmp5 
		<< "," << tmp6
		<< ")");

	if ( iflow != 5)
	{
		init_voltage_sources();

		iflow=5;
		
		if ( processing_state > 3)
			processing_state=3;
	}

	masym = itmp4/10;
	ASSERT(itmp1 >= 0);
	
	m_excitation_type = itmp1;
	if ( (m_excitation_type == EXCITATION_VOLTAGE) || (m_excitation_type == EXCITATION_VOLTAGE_DISC) )  // (Voltage excitation)
	{
		ntsol=0;
	
		if ( m_excitation_type == EXCITATION_VOLTAGE_DISC) // Voltage DISC. 
		{
			nvqd++;
			ivqd.resize(nvqd );
			iqds.resize(nvqd );
			vqd.resize( nvqd );
			vqds.resize(nvqd );
		
			int indx = nvqd-1;
	
			ivqd[indx]= m_geometry->get_segment_number( itmp2, itmp3);
			vqd[indx]= nec_complex( tmp1, tmp2);
			if ( abs( vqd[indx]) < 1.e-20)
				vqd[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 */