nec_context.cpp

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

				symmetry_array.resize(nop*nop);	
				mhz = 1;
				
				/* irngf is not used (NGF function not implemented) */
				if ( imat == 0)
					fblock( npeq, neq, iresrv, m_geometry->m_ipsym);
				
				in_freq_loop = true;
			}
			
			if ( mhz != 1)
			{
				if ( ifrq == 1)
					freq_mhz *= delfrq;
				else
					freq_mhz += delfrq;
			}

			wavelength = CVEL / freq_mhz;

			print_freq_int_krnl(freq_mhz, wavelength, rkh, m_use_exk);
				
			m_geometry->frequency_scale(freq_mhz);
			processing_state = 2;

		case 2: /* structure segment loading */
			structure_segment_loading();

			processing_state=3;
			ntsol=0;
		
		case 3: /* excitation set up (right hand side, -e inc.) */
			nthic=1;
			nphic=1;
			inc=1;
			nprint=0;

		default:
			enum excitation_return ret = excitation_loop(igox, mhz, 
				iptflg, iptflq, 
				iptag, iptagf, iptagt, 
				iptaq, iptaqf, iptaqt, 
				thetis, nfrq, iflow, 
				nthi, nphi, iped);
		
			if (FREQ_LOOP_CONTINUE == ret)
			{
				continue; // Continue frequency loop
			}
			if (FREQ_LOOP_CARD_CONTINUE == ret)
			{
				throw 1; // Continue card input
			}
		
			nphic = 1;
	
			/* normalized receiving pattern printed */
			print_norm_rx_pattern(iptflg, nthi, nphi, thetis, phiss);
			
			xpr2  = phiss;
	
			if ( mhz == nfrq)
				ifar=-1;
	
			if ( nfrq == 1)
			{
				m_output.end_section();
				throw 1; // Continue card input
			}
		
			print_input_impedance(iped, ifrq, nfrq, delfrq);
		
			nfrq=1;
			mhz=1;    
		} /* switch( igox ) */
	}
	while( (++mhz <= nfrq) );
	} /* try */
	catch (int excep)
	{
		ASSERT(excep == 1);
		// keep going on the card input. The exception
		// is thrown in order to continue card input.
	}
}

/* ********************************************************************************************************** */




void nec_context::print_freq_int_krnl(
	nec_float f, 
	nec_float lambda, 
	nec_float int_dist, 
	bool using_extended_kernel)
{
	m_output.end_section();
	m_output.set_indent(31);
	m_output.line("--------- FREQUENCY --------");
	m_output.string("FREQUENCY= "); m_output.real_out(11,4,f); m_output.line(" MHZ");
	m_output.string("WAVELENGTH="); m_output.real_out(11,4,lambda); m_output.line(" METERS");

	m_output.endl(2);
	m_output.set_indent(24);
	m_output.line("APPROXIMATE INTEGRATION EMPLOYED FOR SEGMENTS");
	m_output.string("THAT ARE MORE THAN "); m_output.real_out(5,3,int_dist,false); m_output.line(" WAVELENGTHS APART");

	if ( using_extended_kernel )
		m_output.line( "THE EXTENDED THIN WIRE KERNEL WILL BE USED");
		
	m_output.set_indent(0);
}

void nec_context::antenna_env(void)
{
	m_output.end_section();
	m_output.line("                            -------- ANTENNA ENVIRONMENT --------" );
	
	if ( false == ground.present())
	{
	      m_output.line("                            FREE SPACE" );
		  return;
	}
	
	ground.frati=cplx_10();

	if (false == ground.type_perfect()) // if ( ground.iperf != 1)
	{
		if ( ground.sig < 0.)
			ground.sig=- ground.sig/(59.96*wavelength);

		nec_complex epsc = nec_complex( ground.epsr, -ground.sig*wavelength*59.96);
		ground.zrati = 1.0/ sqrt( epsc);
		
		ground_wave.set_u(ground.zrati);

		if (  ground.radial_wire_count != 0)
		{
			ground.scrwl=  ground.radial_wire_length/ wavelength;
			ground.scrwr=  ground.radial_wire_radius/ wavelength;
			ground.m_t1 = cplx_01()*2367.067/ (nec_float) ground.radial_wire_count;
			ground.t2 = ground.scrwr * (nec_float) ground.radial_wire_count;

			m_output.line(
					"                            RADIAL WIRE GROUND SCREEN");
			m_output.nec_printf(
					"                            %d WIRES\n"
					"                            WIRE LENGTH: %8.2f METERS\n"
					"                            WIRE RADIUS: %10.3E METERS",
					ground.radial_wire_count,  ground.radial_wire_length,  ground.radial_wire_radius );

			m_output.endl();
			m_output.line("                            MEDIUM UNDER SCREEN -" );
		}

		if (false == ground.type_sommerfeld_norton())
		{
			m_output.line("                            FINITE GROUND - REFLECTION COEFFICIENT APPROXIMATION" );
		} 
		else
		{
			// calculate the Sommerfeld Norton ground stuff.
			ggrid.sommerfeld( ground.epsr, ground.sig, freq_mhz );
			
			ground.frati = (epsc-1.0)/(epsc+1.0);
			if ( abs(( ggrid.m_epscf- epsc)/ epsc) >= 1.0e-3 )
			{
				nec_stop("ERROR IN GROUND PARAMETERS -"
					"\n COMPLEX DIELECTRIC CONSTANT FROM FILE IS: %12.5E%+12.5Ej"
					"\n REQUESTED: %12.5E%+12.5Ej",
					real(ggrid.m_epscf), imag(ggrid.m_epscf), 
					real(epsc), imag(epsc) );
			}

			m_output.line("                            FINITE GROUND - SOMMERFELD SOLUTION" );

		} /* if (  ground.type_sommerfeld_norton() ) */

		m_output.endl();
		m_output.nec_printf(
				"                            "
				"RELATIVE DIELECTRIC CONST: %.3f\n"
				"                            "
				"CONDUCTIVITY: %10.3E MHOS/METER\n"
				"                            "
				"COMPLEX DIELECTRIC CONSTANT: %11.4E%+11.4Ej",
				ground.epsr, ground.sig, real(epsc), imag(epsc) );
	}
	else
	{
		m_output.nec_printf("                            PERFECT GROUND" );
	}
}

void nec_context::print_structure_currents(char *pattype, int iptflg, int iptflq,
	 int iptag, int iptagf, int iptagt, int iptaq, int iptaqf, int iptaqt)
{
	int jump;
	nec_float cmag, ph;
	nec_complex curi;
	nec_float fr;
	nec_float etha, ethm, ephm, epha;
	nec_complex eth, eph, ex, ey, ez;

	if ( m_geometry->n != 0)
	{
		if ( iptflg != -1)
		{
			if ( iptflg <= 0)
			{
				m_output.endl(3);
				m_output.line(	"                           -------- CURRENTS AND LOCATION --------");
				m_output.line(	"                                  DISTANCES IN WAVELENGTHS" );
				m_output.endl();
				m_output.line(	"   SEG  TAG    COORDINATES OF SEGM CENTER     SEGM    ------------- CURRENT (AMPS) -------------");
				m_output.line(	"   No:  No:       X         Y         Z      LENGTH     REAL      IMAGINARY    MAGN        PHASE");
			}
			else if ( (iptflg != 3) && (inc <= 1) )
			{
				m_output.endl(3);
				m_output.nec_printf(
						"             -------- RECEIVING PATTERN PARAMETERS --------\n"
						"                      ETA: %7.2f DEGREES\n"
						"                      TYPE: %s\n"
						"                      AXIAL RATIO: %6.3f\n\n"
						"            THETA     PHI      ----- CURRENT ----    SEG\n"
						"            (DEG)    (DEG)     MAGNITUDE    PHASE    No:",
						xpr3, pattype, xpr6 );
			} /* if ( iptflg <= 0) */
		} /* if ( iptflg != -1) */

		structure_power_loss=0.;
		int itmp1=0;
		jump= iptflg+1;

		for (int i = 0; i < m_geometry->n; i++ )
		{
			curi= current_vector[i]* wavelength;
			cmag= abs( curi);
			ph= arg_degrees( curi);

			if ( (nload != 0) && (fabs(real(zarray[i])) >= 1.e-20) )
				structure_power_loss += 0.5*cmag*cmag*real( zarray[i]) * m_geometry->segment_length[i];

			if ( jump == 0)
				continue;

			if ( jump > 0 )
			{
				if ( (iptag != 0) && (m_geometry->segment_tags[i] != iptag) )
					continue;

				itmp1++;
				if ( (itmp1 < iptagf) || (itmp1 > iptagt) )
					continue;

				if ( iptflg != 0)
				{
					if ( iptflg >= 2 )
					{
						fnorm[inc-1]= cmag;
						isave= (i+1);
					}

					if ( iptflg != 3)
					{
						m_output.endl();
						m_output.nec_printf("          %7.2f  %7.2f   %11.4E  %7.2f  %5d",
							xpr1, xpr2, cmag, ph, i+1 );

						continue;
					}
				} /* if ( iptflg != 0) */
			}
			else
			{
				m_output.endl();
				m_output.nec_printf(
					" %5d %4d %9.4f %9.4f %9.4f %9.5f"
					" %11.4E %11.4E %11.4E %8.3f",
					i+1, m_geometry->segment_tags[i],
					m_geometry->x[i], m_geometry->y[i], m_geometry->z[i], m_geometry->segment_length[i],
					real(curi), imag(curi), cmag, ph );

				// added test for plot_card.is_valid()
				if (plot_card.is_valid() && plot_card.currents())
				{
					plot_card.plot_complex(curi);
					plot_card.plot_endl();
				}
			}

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

		if ( iptflq != -1)
		{
			m_output.endl(3);
			m_output.nec_printf(
			    "                                  "
			    "------ CHARGE DENSITIES ------\n"
			    "                                  "
			    "   DISTANCES IN WAVELENGTHS\n\n"
			    "   SEG   TAG    COORDINATES OF SEG CENTER     SEG"
			    "        "
			    "  CHARGE DENSITY (COULOMBS/METER)\n"
			    "   NO:   NO:     X         Y         Z       LENGTH"
			    "   "
			    "  REAL      IMAGINARY     MAGN        PHASE" );

			itmp1 = 0;
			fr = 1.e-6/freq_mhz;

			for(int i = 0; i < m_geometry->n; i++ )
			{
				if ( iptflq != -2 )
				{
					if ( (iptaq != 0) && (m_geometry->segment_tags[i] != iptaq) )
						continue;

					itmp1++;
					if ( (itmp1 < iptaqf) || (itmp1 > iptaqt) )
						continue;

				} /* if ( iptflq == -2) */

				curi = fr * nec_complex(- bii[i], bir[i]);
				cmag = abs( curi);
				ph = arg_degrees( curi);

				m_output.endl();
				m_output.nec_printf(
						" %5d %4d %9.4f %9.4f %9.4f %9.5f"
						" %11.4E %11.4E %11.4E %9.3f",
						i+1, m_geometry->segment_tags[i], m_geometry->x[i], m_geometry->y[i], m_geometry->z[i], m_geometry->segment_length[i],
						real(curi), imag(curi), cmag, ph );

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

		} /* if ( iptflq != -1) */

	} /* if ( n != 0) */

	if ( m_geometry->m != 0)
	{
		m_output.endl(3);
		m_output.nec_printf(
		    "                                      "
		    " --------- SURFACE PATCH CURRENTS ---------\n"
		    "                                                "
		    " DISTANCE IN WAVELENGTHS\n"
		    "                                                "
		    " CURRENT IN AMPS/METER\n\n"
		    "                                 ---------"
		    " SURFACE COMPONENTS --------    "
		    "---------------- RECTANGULAR COMPONENTS ----------------\n"
		    "  PCH   --- PATCH CENTER ---     TANGENT VECTOR 1    "
		    " TANGENT VECTOR 2    ------- X ------    ------- Y ------"
		    "   "
		    " ------- Z ------\n  No:    X       Y       Z       MAG."
		    "       "
		    "PHASE     MAG.       PHASE    REAL   IMAGINARY    REAL  "
		    " IMAGINARY    REAL   IMAGINARY" );

		int j = m_geometry->n-3;
		int itmp1 = -1;

		for(int i = 0; i < m_geometry->m; i++ )
		{
			j += 3;
			itmp1++;
			ASSERT(itmp1 == i);
			
			ex= current_vector[j];
			ey= current_vector[j+1];
			ez= current_vector[j+2];
			eth= ex* m_geometry->t1x[itmp1]+ ey* m_geometry->t1y[itmp1]+ ez* m_geometry->t1z[itmp1];
			eph= ex* m_geometry->t2x[itmp1]+ ey* m_geometry->t2y[itmp1]+ ez* m_geometry->t2z[itmp1];
			ethm= abs( eth);
			etha= arg_degrees( eth);
			ephm= abs( eph);
			epha= arg_degrees( eph);

			m_output.endl();
			m_output.nec_printf(
			      " %4d %7.3f %7.3f %7.3f %11.4E "
			      "%8.2f %11.4E %8.2f"
			      " %9.2E %9.2E %9.2E %9.2E %9.2E %9.2E",
			      i+1, m_geometry->px[itmp1], m_geometry->py[itmp1], m_geometry->pz[itmp1],
			      ethm, etha, ephm, epha, real(ex), imag(ex),
			      real(ey), imag(ey), real(ez), imag(ez));

		  	plot_card.plot_currents(ex,ey,ez);
		} /* for( i=0; i<m; i++ ) */
	} /* if ( m != 0) */
} /* print_structure_currents */



/*!\brief Calculate network data such as the lengths of transmission lines.
	TODO Fix up this horrible mess. As far as I can figure we should simply
	be calculating the x11i[] value for every network with ntyp = 2 or 3 (transmission lines).
	The bug being reported by Remi, could be caused by this network failing to calculate the
	length correctly under some circumstances.

			// This is an extremely strange statement.
			//	ntyp[j]		net_type	ntyp[j]/net_type
			//	1			1				1
			//	2			1				2
			//	3			1				3
			//	1			2				0
			//	2			2				1
			//	3			2				1
*/
void nec_context::calculate_network_data(void)
{
	if ( (network_count == 0) || (inc > 1) )
		return;
		
	int itmp3 = 0;
	int net_type = ntyp[0];

	for (int i = 0; i < 2; i++ )
	{
		if ( net_type == 3)
			net_type = 2;

		for (int j = 0; j < network_count; j++)
		{
			if ( (ntyp[j]/net_type) != 1 )
			{
				itmp3 = ntyp[j]; // can never be zero
			}
			else
			{
				if ( (ntyp[j] >= 2) && (x11i[j] <= 0.0) )
				{
					int idx4 = iseg1[j]-1;
					int idx5 = iseg2[j]-1;
					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);
				}
			}
		}
		if ( itmp3 == 0)	// can only be zero if all the networks are the same type
			return;
		net_type = itmp3;
	}
} /* calculate_network_data */


void nec_context::print_network_data(void)
{
//	int i, j;
	int itmp1, itmp2, itmp3, itmp4, itmp5;
	char *pnet[3] = { "        ", "STRAIGHT", " CROSSED" };
	
	if ( (network_count != 0) && (inc <= 1) )
	{
		m_output.nec_printf( "\n\n\n"
		    "                                            "
		    "---------- NETWORK DATA ----------" );

		itmp3=0;
		itmp1= ntyp[0];

		for(int i = 0; i < 2; i++ )
		{
			if ( itmp1 == 3)
				itmp1=2;

			if ( itmp1 == 2)
			{
				m_output.endl();
				m_output.nec_printf(
						"  -- FROM -  --- TO --      "
						"TRANSMISSION LINE       "
						" --------- SHUNT ADMITTANCES (MHOS) "
						"---------   LINE\n"
						"  TAG   SEG  TAG   SEG    IMPEDANCE      "
						"LENGTH    "
						" ----- END ONE -----      "
						"----- END TWO -----   TYPE\n"
						"  No:   No:  No:   No:         OHMS      "
						"METERS      REAL      IMAGINARY      "
						"REAL      IMAGINARY" );
			}
			else if (itmp1 == 1)
			{
				m_output.endl();
				m_output.nec_printf(
						"  -- FROM -  --- TO --            "
						"--------"
						" ADMITTANCE MATRIX ELEMENTS (MHOS) "
						"---------\n"
						"  TAG   SEG  TAG   SEG   "
						"----- (ONE,ONE) ------  "
						" ----- (ONE,TWO) -----   "
						"----- (TWO,TWO) -------\n"
						"  No:   No:  No:   No:      REAL      "
						"IMAGINARY     "
						" REAL     IMAGINARY       REAL      "
						"IMAGINARY" );
			}
			for (int j = 0; j < network_count; j++)
			{
				itmp2= ntyp[j];

				if ( (itmp2/itmp1) != 1 )
					itmp3 = itmp2;
				else
				{