nec_radiation_pattern.cpp

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

	
	int result_counter=0;
	nec_float pint = 0.0;
	nec_float delta_phi_rad = degrees_to_rad(delta_phi);
	nec_float tmp2 = 0.5 * degrees_to_rad(delta_theta);
	phi= m_phi_start- delta_phi;

	for (int kph = 1; kph <= n_phi; kph++ )
	{
		phi += delta_phi;
		nec_float pha = degrees_to_rad(phi);
		thet= m_theta_start- delta_theta;
		
		for (int kth = 1; kth <= n_theta; kth++ )
		{
			thet += delta_theta;
			if ( m_ground.present() && (thet > 90.01) && (_ifar != 1) )
				continue;
		
			nec_float tha = degrees_to_rad(thet);
			
			nec_complex  eth, eph;
			if ( 1 == _ifar)
			{
				bool space_wave_only = (false == m_ground.present());
				nec_complex erd;
				
				m_context->gfld(m_range/_wavelength, pha, thet/_wavelength,
					&eth, &eph, &erd, space_wave_only, _wavelength );
			
				_e_theta[result_counter] = eth;
				_e_phi[result_counter] = eph;
				_e_r[result_counter] = erd;
			}
			else
			{
				m_context->ffld(tha, pha, &eth, &eph, _wavelength);
		
				nec_float ethm2= norm(eth);
				nec_float ethm= sqrt(ethm2);
				nec_float etha= arg_degrees(eth);
				nec_float ephm2= norm(eph);
				nec_float ephm= sqrt( ephm2);
				nec_float epha= arg_degrees( eph);
		
				/* elliptical polarization calc. */
				if ( (ethm2 <= 1.0e-20) && (ephm2 <= 1.0e-20) )
				{
					tilta=0.;
					emajr2=0.;
					eminr2=0.;
					pol_axial_ratio=0.;
					pol_sense_index = 3;
				}
				else
				{
					dfaz= epha- etha;
					if ( epha >= 0.)
						dfaz2= dfaz-360.;
					else
						dfaz2= dfaz+360.;
				
					if ( fabs(dfaz) > fabs(dfaz2) )
						dfaz= dfaz2;
				
					cdfaz= cos(degrees_to_rad(dfaz));
					tstor1= ethm2- ephm2;
					tstor2=2.* ephm* ethm* cdfaz;
					tilta=.5* atan2( tstor2, tstor1);
					stilta= sin( tilta);
					tstor1= tstor1* stilta* stilta;
					tstor2= tstor2* stilta* cos( tilta);
					emajr2=- tstor1+ tstor2+ ethm2;
					eminr2= tstor1- tstor2+ ephm2;
					if ( eminr2 < 0.)
						eminr2=0.;
				
					pol_axial_ratio= sqrt( eminr2/ emajr2);
					tilta= rad_to_degrees(tilta);
					
					if ( pol_axial_ratio <= 1.0e-5)
						pol_sense_index = 0;
					else
					if ( dfaz <= 0.)
						pol_sense_index = 1;
					else
						pol_sense_index = 2;
				
				} /* if ( (ethm2 <= 1.0e-20) && (ephm2 <= 1.0e-20) ) */
			
				/* Gain Normalization Variables */
				nec_float gnmj= db10( gcon* emajr2);
				nec_float gnmn= db10( gcon* eminr2);
				nec_float gnv = db10( gcon* ethm2);
				nec_float gnh = db10( gcon* ephm2);
				nec_float gtot= db10( gcon*(ethm2+ ephm2) );
			
				if (m_rp_normalization > 0)
				{
					nec_float temp_gain;
					
					switch(m_rp_normalization )
					{
					case 1:
						temp_gain = gnmj;
						break;
				
					case 2:
						temp_gain = gnmn;
						break;
				
					case 3:
						temp_gain = gnv;
						break;
				
					case 4:
						temp_gain = gnh;
						break;
				
					case 5:
						temp_gain = gtot;
						break;
						
					default:
						throw new nec_exception("Unknown Gain Normalization Encountered.");
					}
				
					_gain[result_counter] = temp_gain;
				
				} /* if ( m_rp_normalization > 0) */
			
				if ( m_rp_power_average != 0)
				{
					// compute the numerical integral of the  power gain in angular co-ordinates
					tstor1= gcop*( ethm2+ ephm2);
					nec_float tmp3 = tha - tmp2;
					nec_float tmp4 = tha + tmp2;
				
					if ( kth == 1)
						tmp3= tha;
					else if ( kth == n_theta)
						tmp4= tha;
				
					nec_float da = fabs( delta_phi_rad*( cos(tmp3)- cos(tmp4)));
					if ( (kph == 1) || (kph == n_phi) )
						da *=.5;
					pint += tstor1 * da;
				
					if ( m_rp_power_average == 2) // do not print the power gain values (just compute the average)
						continue;
				}
			
				if ( m_rp_output_format != 1)
				{
					_power_gain_vert[result_counter] = gnmj;
					_power_gain_horiz[result_counter] = gnmn;
				}
				else
				{
					_power_gain_vert[result_counter] = gnv;
					_power_gain_horiz[result_counter] = gnh;
				}
			
				ethm *= _wavelength;
				ephm *= _wavelength;
			
				if ( m_range >= 1.0e-20 )
				{
					ethm= ethm* exrm;
					etha= etha+ exra;
					ephm= ephm* exrm;
					epha= epha+ exra;
				}
			
				_e_theta[result_counter] = deg_polar(ethm, etha);
				_e_phi[result_counter] = deg_polar(ephm, epha);
				
				_power_gain_tot[result_counter] = gtot;
				_polarization_axial_ratio[result_counter] = pol_axial_ratio;
				_polarization_tilt[result_counter] = tilta;
				_polarization_sense_index[result_counter] = pol_sense_index;
				
			} /* if ( _ifar == 1) */
			result_counter++;
		} /* for( kth = 1; kth <= n_theta; kth++ ) */
	
	} /* for( kph = 1; kph <= n_phi; kph++ ) */

	if ( m_rp_power_average != 0)
	{
/*		nec_float tmp3 = degrees_to_rad(m_theta_start);
		tmp4 = tmp3 + degrees_to_rad(delta_theta) * (nec_float)(n_theta-1);
		tmp3 = fabs( degrees_to_rad(delta_phi) * (nec_float)( n_phi-1)*( cos( tmp3)- cos( tmp4)));
		pint /= tmp3;
		tmp3 /= pi();
		_average_power_gain = pint;
		_average_power_solid_angle = tmp3; */

		// We now compute the solid angle over which the power is averaged
		nec_float theta_start_rad = degrees_to_rad(m_theta_start);
		nec_float theta_range = degrees_to_rad(delta_theta) * (nec_float)(n_theta-1);
		nec_float phi_range = degrees_to_rad(delta_phi) * (nec_float)(n_phi-1);
		nec_float total_theta = theta_start_rad + theta_range;
		nec_float solid_angle = fabs(phi_range * (cos(theta_start_rad) - cos(total_theta)) );
						
		_average_power_gain = pint / solid_angle;
		_average_power_solid_angle = solid_angle / pi(); // We display it as a multiple of pi()
	}

	_maximum_gain = _gain.max();
	m_analysis_done = true;
}


nec_float nec_radiation_pattern::get_gain_normalization_factor(nec_float gnor)
{
	if ( fabs(gnor) > 1.0e-20)
		return gnor;
		
	if (false == m_analysis_done)
		throw new nec_exception("Internal Error: Radiation Pattern Analysis not done");
		
	return _maximum_gain;
}

void nec_radiation_pattern::write_normalized_gain(ostream& os)
{
	// if  is non-zero then use it as a normaliztion factor
	
	nec_float normalization_factor = get_gain_normalization_factor(m_rp_gnor);
		
	string norm_type;
	switch (m_rp_normalization)
	{
		case 1:
			norm_type = "  MAJOR AXIS"; break;
		case 2:
			norm_type = "  MINOR AXIS"; break;
		case 3:
			norm_type = "    VERTICAL"; break;
		case 4:
			norm_type = "  HORIZONTAL"; break;
		case 5:
			norm_type = "      TOTAL "; break;
		
		default: throw new nec_exception("Unknown Gain Normalization Encountered.");
	}
	
	output_helper oh(os,_result_format);
	
	oh.section_start();
	os << "                              ---------- NORMALIZED GAIN ----------" << endl;
	os << "                                      " << norm_type << " GAIN" << endl;
	os << "                                   NORMALIZATION FACTOR: "; oh.real_out(7,2,normalization_factor,false); os << " db" << endl << endl;
	os << "    ---- ANGLES ----                ---- ANGLES ----                ---- ANGLES ----" << endl;
	os << "    THETA      PHI        GAIN      THETA      PHI        GAIN      THETA      PHI       GAIN" << endl;
	os << "   DEGREES   DEGREES        DB     DEGREES   DEGREES        DB     DEGREES   DEGREES       DB" << endl;

	
	int row_count = 0;
	int n_cols = 3;
	
	int item_count = 0;
	for (int p=0;p<n_phi;p++)
	{
		nec_float phi = m_phi_start + p*delta_phi ;
		for (int t=0;t<n_theta;t++)
		{

			nec_float theta = m_theta_start + t * delta_theta;
			nec_float norm_gain = _gain[item_count++] - normalization_factor;
			
			oh.start_record();
			oh.padding(" ");
			
			oh.real_out(9,2,theta,false);
			oh.separator(); oh.real_out(9,2,phi,false);
			oh.separator(); oh.real_out(9,2,norm_gain,false);
			oh.padding(" ");
			
			if (_result_format == RESULT_FORMAT_NEC)
			{
				if (item_count % n_cols == 0)
				{
					row_count++;
					oh.end_record();
				}
			}
			else
			{
				oh.end_record();
			}
		}
	}
	os << endl;
}