nec_context.cpp

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

					int idx4, idx5;

					itmp4= iseg1[j];
					itmp5= iseg2[j];
					idx4 = itmp4-1;
					idx5 = itmp5-1;

#if 0
					This code has moved to calculate_network_data()
					
					if ( (itmp2 >= 2) && (x11i[j] <= 0.0) )
					{
						nec_float xx = m_geometry->x[idx5]- m_geometry->x[idx4];
						nec_float yy = m_geometry->y[idx5]- m_geometry->y[idx4];
						nec_float zz = m_geometry->z[idx5]- m_geometry->z[idx4];
						
						// set the length of the transmission line to be the 
						// straight line distance.
						x11i[j]= wavelength*sqrt(xx*xx+yy*yy+zz*zz);
					}
#endif
					m_output.endl();
					m_output.nec_printf(
						" %4d %5d %4d %5d  "
						"%11.4E %11.4E  %11.4E %11.4E  "
						"%11.4E %11.4E %s",
						m_geometry->segment_tags[idx4], itmp4, m_geometry->segment_tags[idx5], itmp5,
						x11r[j], x11i[j], x12r[j], x12i[j],
						x22r[j], x22i[j], pnet[itmp2-1] );

				} /* if (( itmp2/ itmp1) == 1) */

			} /* for( j = 0; j < network_count; j++) */
			if ( itmp3 == 0)
				break;
			itmp1= itmp3;
		} /* for( j = 0; j < network_count; j++) */
	} /* if ( (network_count != 0) && (inc <= 1) ) */
} /* print_network_data */

void nec_context::print_norm_rx_pattern(int iptflg, int nthi, int nphi,
	nec_float thetis, nec_float phiss)
{
	if ( iptflg < 2)
		return; // do not print
		
	{
		// Call the new nec_results class that handles
		// a normalized receiving pattern.
		
		nec_float theta_step = xpr4;
		nec_float phi_step = xpr5;
		nec_float eta = xpr3;
		nec_float axial_ratio = xpr6;
		string pol_type( hpol[m_excitation_type-1]);
		int segment_number = isave;
		
		nec_norm_rx_pattern* rx_pattern = new nec_norm_rx_pattern(nthi, nphi,
			fnorm,
			thetis, theta_step,
			phiss, phi_step,
			eta, 
			axial_ratio, 
			segment_number, 
			pol_type);
			
		m_results.add(rx_pattern);
		
		// write the restuls to our output.
		std::stringstream ss;
		rx_pattern->write_to_file(ss);
		m_output.line(ss.str().c_str());
	}
	
#if 0
	nec_float phi = phiss;
	int itmp1 = nthi * nphi;

	nec_float norm_factor = fnorm[0];
	
	for(int j = 1; j < itmp1; j++ )
		if ( fnorm[j] > norm_factor)
			norm_factor = fnorm[j];

	m_output.endl(3);
	m_output.nec_printf(
		"                     "
		"---- NORMALIZED RECEIVING PATTERN ----\n"
		"                      "
		"NORMALIZATION FACTOR: %11.4E\n"
		"                      "
		"ETA: %7.2f DEGREES\n"
		"                      "
		"TYPE: %s\n"
		"                      AXIAL RATIO: %6.3f\n"
		"                      SEGMENT No: %d\n\n"
		"                      "
		"THETA     PHI       ---- PATTERN ----\n"
		"                      "
		"(DEG)    (DEG)       DB     MAGNITUDE",
		norm_factor, xpr3, hpol[m_excitation_type-1], xpr6, isave );
	
	for(int j = 0; j < nphi; j++ )
	{
		int itmp2 = nthi*j;

		for(int i = 0; i < nthi; i++ )
		{
			int itmp3 = i + itmp2;

			if ( itmp3 < itmp1)
			{
				nec_float _tmp2 = fnorm[itmp3] / norm_factor;
				nec_float _tmp3 = db20( _tmp2);

				m_output.endl();
				m_output.nec_printf("                    %7.2f  %7.2f   %7.2f  %11.4E",
					xpr1, phi, _tmp3, _tmp2 );

				xpr1 += xpr4;
			}

		} /* for( i = 0; i < nthi; i++ ) */

		xpr1= thetis;
		phi += xpr5;
	} /* for( j = 0; j < nphi; j++ ) */
#endif
}

void nec_context::print_power_budget(void)
{
	if ( (m_excitation_type == EXCITATION_VOLTAGE) || (m_excitation_type == EXCITATION_VOLTAGE_DISC) )
	{
		nec_float  radiated_power = input_power- network_power_loss - structure_power_loss;
		nec_float  efficiency = 100.0 * radiated_power/input_power;

		m_output.endl(3);
		m_output.nec_printf(
		    "                               "
		    "---------- POWER BUDGET ---------\n"
		    "                               "
		    "INPUT POWER   = %11.4E Watts\n"
		    "                               "
		    "RADIATED POWER= %11.4E Watts\n"
		    "                               "
		    "STRUCTURE LOSS= %11.4E Watts\n"
		    "                               "
		    "NETWORK LOSS  = %11.4E Watts\n"
		    "                               "
		    "EFFICIENCY    = %7.2f Percent",
		    input_power, radiated_power, structure_power_loss, network_power_loss, efficiency );
	} /* if ( (m_excitation_type == EXCITATION_VOLTAGE) || (m_excitation_type == EXCITATION_VOLTAGE_DISC) ) */
}

void nec_context::print_input_impedance(int iped, int ifrq, int nfrq, nec_float delfrq)
{
	nec_float tmp1, tmp2, tmp3, tmp4, tmp5;
	int i, itmp1, itmp2;

	if ( iped != 0)
	{
		int iss;

		if ( nvqd > 0)
			iss = ivqd[nvqd-1];
		else
			iss = source_segment_array[voltage_source_count-1];

		m_output.endl(3);
		m_output.nec_printf(
		    "                            "
		    " -------- INPUT IMPEDANCE DATA --------\n"
		    "                                     "
		    " SOURCE SEGMENT No: %d\n"
		    "                                  "
		    " NORMALIZATION FACTOR:%12.5E\n\n"
		    "              ----------- UNNORMALIZED IMPEDANCE ----------  "
		    "  ------------ NORMALIZED IMPEDANCE -----------\n"
		    "      FREQ    RESISTANCE    REACTANCE    MAGNITUDE    PHASE  "
		    "  RESISTANCE    REACTANCE    MAGNITUDE    PHASE\n"
		    "       MHz       OHMS         OHMS         OHMS     DEGREES  "
		    "     OHMS         OHMS         OHMS     DEGREES",
		    iss, impedance_norm_factor );

		itmp1= nfrq;
		if ( 0 == ifrq )
			tmp1= freq_mhz-( nfrq-1)* delfrq;
		else
			tmp1= freq_mhz/( pow(delfrq, (nfrq-1)) );

		for( i = 0; i < itmp1; i++ )
		{
			itmp2= 4*i;
			tmp2= fnorm[itmp2  ]/ impedance_norm_factor;
			tmp3= fnorm[itmp2+1]/ impedance_norm_factor;
			tmp4= fnorm[itmp2+2]/ impedance_norm_factor;
			tmp5= fnorm[itmp2+3];

			m_output.endl();
			m_output.nec_printf(
			    " %9.3f   %11.4E  %11.4E  %11.4E  %7.2f  "
			    " %11.4E  %11.4E  %11.4E  %7.2f",
			    tmp1, fnorm[itmp2], fnorm[itmp2+1], fnorm[itmp2+2],
			    fnorm[itmp2+3], tmp2, tmp3, tmp4, tmp5 );

			if ( ifrq == 0)
				tmp1 += delfrq;
			else
				tmp1 *= delfrq;
		} /* for( i = 0; i < itmp1; i++ ) */
		m_output.end_section();
	} /* if ( iped != 0) */
}


void nec_context::structure_segment_loading()
{
	nec_float tim1, tim, tim2;

	m_output.end_section();
	m_output.line("                          ------ STRUCTURE IMPEDANCE LOADING ------" );

	if ( nload != 0)
		load();

	if ( nload == 0 )
	{	
		m_output.line("                                 THIS STRUCTURE IS NOT LOADED" );
	}
	antenna_env();

	/* fill and factor primary interaction matrix */
	secnds( &tim1 );
	cmset( neq, cm, rkh );
	secnds( &tim2 );
	tim= tim2- tim1;
	factrs(m_output, npeq, neq, cm, ip );
	secnds( &tim1 );
	tim2= tim1- tim2;
	
	m_output.end_section();
	m_output.line("                             ---------- MATRIX TIMING ----------");
	m_output.string("                               FILL= ");
	m_output.integer(int(tim)); 
	m_output.string(" msec  FACTOR: "); 
	m_output.integer(int(tim2)); m_output.string(" msec");
}


void nec_context::setup_excitation(int iptflg)
{
	nec_float tmp1, tmp2, tmp3, tmp4, tmp5, tmp6;

	tmp1=tmp2=tmp3=tmp4=tmp5=tmp6=0.0;

	if ( (m_excitation_type != 0) && (m_excitation_type != 5) )
	{
		if ( (iptflg <= 0) || (m_excitation_type == 4) )
		{
			m_output.endl(3);
			m_output.line("                             ---------- EXCITATION ----------" );
		}
		
		tmp5= degrees_to_rad(xpr5);
		tmp4= degrees_to_rad(xpr4);

		if ( m_excitation_type == 4)
		{
			tmp1= xpr1/ wavelength;
			tmp2= xpr2/ wavelength;
			tmp3= xpr3/ wavelength;
			tmp6= xpr6/( wavelength* wavelength);

			m_output.endl();
			m_output.line(		"                                      CURRENT SOURCE");
			m_output.line(		"                     -- POSITION (METERS) --       ORIENTATION (DEG)");
			m_output.line(		"                     X          Y          Z       ALPHA        BETA   DIPOLE MOMENT");
			m_output.nec_printf("               %10.5f %10.5f %10.5f  %7.2f     %7.2f    %8.3f",
				xpr1, xpr2, xpr3, xpr4, xpr5, xpr6 );
		}
		else
		{
			tmp1= degrees_to_rad(xpr1);
			tmp2= degrees_to_rad(xpr2);
			tmp3= degrees_to_rad(xpr3);
			tmp6= xpr6;

			if ( iptflg <= 0)
			{
				m_output.endl();
				m_output.nec_printf(
					"PLANE WAVE - THETA: %7.2f deg, PHI: %7.2f deg,"
					" ETA=%7.2f DEG, TYPE - %s  AXIAL RATIO: %6.3f",
					xpr1, xpr2, xpr3, hpol[m_excitation_type-1], xpr6 );
			}
		} /* if ( m_excitation_type == 4) */
	}

	nec_float incident_amplitude = 1.0; // the incident electric field amplitude...
	if (xpr7 != 0.0)
		incident_amplitude = xpr7;
		

	/* Fill E-field right-hand matrix */
	etmns( tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, incident_amplitude, m_excitation_type, current_vector);
}

enum excitation_return nec_context::excitation_loop(int in_freq_loop_state, int mhz, 
	int iptflg, int iptflq, int iptag, int iptagf, int iptagt,
	int iptaq, int iptaqf, int iptaqt, 
	nec_float thetis, 
	int nfrq, int iflow, int nthi, int nphi, int iped)
{
	int itmp1;

	do
	{
		if (in_freq_loop_state < 4)
		{
			setup_excitation(iptflg);

			/* matrix solving  (netwk calls solves) */
			calculate_network_data();
			print_network_data();

			if ( (inc > 1) && (iptflg > 0) )
				nprint=1;

			// c_network::net_solve( cm, &cm[ib11], &cm[ic11], &cm[id11], ip, current_vector );
			netwk( cm, ip, current_vector );
			ntsol=1;

			if ( iped != 0)
			{
				itmp1= 4*( mhz-1);

				fnorm[itmp1  ] = real( zped);
				fnorm[itmp1+1] = imag( zped);
				fnorm[itmp1+2] = abs( zped);
				fnorm[itmp1+3] = arg_degrees( zped);

				if ( iped != 2 )
				{
			 		 if ( fnorm[itmp1+2] > impedance_norm_factor)
			 			impedance_norm_factor = fnorm[itmp1+2];
				}
			} /* if ( iped != 0) */

			/* printing structure currents */
			print_structure_currents(hpol[m_excitation_type-1], iptflg, iptflq,
				iptag, iptagf, iptagt, iptaq, iptaqf, iptaqt);

			print_power_budget();
	     
			processing_state = 4;

			if ( ncoup > 0)
				couple( current_vector, wavelength );

			if ( iflow == 7)
			{
				if ( (m_excitation_type > 0) && (m_excitation_type < 4) )
				{
					nthic++;
					inc++;
					xpr1 += xpr4;

					if ( nthic <= nthi )
						continue; /* continue excitation loop */

					nthic=1;
					xpr1= thetis;
					xpr2= xpr2+ xpr5;
					nphic++;

					if ( nphic <= nphi )
						continue; /* continue excitation loop */

					return FREQ_PRINT_NORMALIZATION;
				} /* if ( (m_excitation_type >= 1) && (m_excitation_type <= 3) ) */

				if ( nfrq != 1) 
				{
					return FREQ_LOOP_CONTINUE; /* continue the freq loop */
				}

				m_output.end_section();
				return FREQ_LOOP_CARD_CONTINUE; /* continue card input loop */

			} /*if ( iflow == 7) */

		}
		
		if (in_freq_loop_state < 5)
			processing_state = 5;

		
		/* near field calculation */
		if (in_freq_loop_state < 6)
		{
			if ( m_near != -1)
			{
				nfpat();

				if ( mhz == nfrq)
					m_near=-1;

				if ( nfrq == 1)
				{
					m_output.end_section();
				  	return FREQ_LOOP_CARD_CONTINUE; /* continue card input loop */
				}
			} /* if ( near != -1) */
		}
		
		/* Standard far-field calculation */
		if ( ifar != -1)
		{
			nec_radiation_pattern* rad_pat =
				new nec_radiation_pattern(nth, nph,
						thets, phis, dth, dph,
						rfld, ground,
						ifar, wavelength,
						input_power, network_power_loss,
						m_rp_output_format, m_rp_normalization, ipd, iavp,
						gnor, plot_card);

			rad_pat->analyze(this);
						
			/* Here we write the largest gain to the standard output
				if the user has specified this on the command line.
			*/
			if (m_output_flags.get_gain_flag())
			{
				rad_pat->write_gain_normalization();
				delete rad_pat;
			}
			else
			{
				m_results.add(rad_pat);
				
				// write the results to the ordinary NEC output file.
				
				std::stringstream ss;
				rad_pat->write_to_file(ss);
				m_output.line(ss.str().c_str());
				// print_radiation_pattern(input_power, network_power_loss);
			}
		}

		if ( (m_excitation_type == 0) || (m_excitation_type >= 4) )
		{
			if ( mhz == nfrq )
				ifar=-1;

			if ( nfrq != 1)
			{
				return FREQ_LOOP_CONTINUE;  /* continue the freq loop */
			}

			m_output.end_section();
			return FREQ_LOOP_CARD_CONTINUE;  /* continue card input loop */
		}

		nthic++;
		inc++;
		xpr1 += xpr4;

		if ( nthic <= nthi )
			continue; /* continue excitation loop */

		nthic = 1;
		xpr1  = thetis;
		xpr2 += xpr5;
		nphic++;

		if ( nphic > nphi )
			return FREQ_PRINT_NORMALIZATION;

	}
	while( true );
} /* excitation_loop */



/* load calculates the impedance of specified */
/* segments for various types of loading */
void nec_context::load()
{
	int istep, istepx, l1, l2, ldtags, ichk;
	bool iwarn = false;
	nec_complex zt, tpcj;
	
	int n = m_geometry->n;
	int np = m_geometry->np;
	
	tpcj = nec_complex(0.0,1.883698955e+9);
	m_output.endl();
	m_output.line("  LOCATION        RESISTANCE  INDUCTANCE  CAPACITANCE     IMPEDANCE (OHMS)   CONDUCTIVITY  CIRCUIT");
	m_output.line("  ITAG FROM THRU     OHMS       HENRYS      FARADS       REAL     IMAGINARY   MHOS/METER      TYPE");
	
	/* Initialize impedance array, used for temporary */
	/* storage of loading information. */
	zarray.resize(n);
	zarray.fill(cplx_00());
	
	istep = 0;
	
	/* Surely This should be rewritten as...
		for (int istepx=0; istepx<nload; istepx++)
		{
			if ( ldtyp[istepx] > 5 )
			{
				m_output.nec_printf(
					"\n  IMPROPER LOAD TYPE CHOSEN,"
					" REQUESTED TYPE IS %d", ldtyp[istepx] );