nec_radiation_pattern.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_radiation_pattern.h"
#include "nec_context.h"
#include "c_geometry.h"
#include "nec_exception.h"

int nec_radiation_pattern::get_index(int theta_index, int phi_index) const
{
	if (theta_index >= n_theta) 
		throw new nec_exception("nec_radiation_pattern: Theta index too large");
	if (phi_index >= n_phi) 
		throw new nec_exception("nec_radiation_pattern: Phi index too large");
		
	return phi_index*n_theta + theta_index;
}

nec_radiation_pattern::nec_radiation_pattern(int in_n_theta, int in_n_phi,
	nec_float in_theta_start, nec_float in_phi_start,
	nec_float in_delta_theta, nec_float in_delta_phi,
	nec_float in_range,
	nec_ground& in_ground,
	int in_ifar, nec_float in_wavelength,
	nec_float pinr, nec_float pnlr,
	int in_rp_output_format, int in_rp_normalization, int in_rp_ipd, int in_rp_power_average,
	nec_float in_gnor,
	c_plot_card& in_plot_card)
		: m_ground(in_ground), m_plot_card(in_plot_card)
{
	n_theta = in_n_theta;
	n_phi = in_n_phi;
	
	m_theta_start = in_theta_start;
	m_phi_start = in_phi_start;
	
	delta_theta = in_delta_theta;
	delta_phi = in_delta_phi;
	
	m_range = in_range; // was rfld
	
	m_rp_output_format = in_rp_output_format;
	m_rp_normalization = in_rp_normalization;
	m_rp_ipd = in_rp_ipd;
	m_rp_power_average = in_rp_power_average; // was iavp
	
	m_rp_gnor = in_gnor;
	
	int n_angles = n_theta * n_phi;
	
	_gain.resize(n_angles);
	_power_gain_vert.resize(n_angles);
	_power_gain_horiz.resize(n_angles);
	_power_gain_tot.resize(n_angles);
	_polarization_axial_ratio.resize(n_angles);
	_polarization_tilt.resize(n_angles);
	_polarization_sense_index.resize(n_angles);
			
	_e_theta.resize(n_angles);
	_e_phi.resize(n_angles);
	_e_r.resize(n_angles);
	
	_ifar = in_ifar;
	_wavelength = in_wavelength;
	_pinr = pinr;
	_pnlr = pnlr;
	
	m_analysis_done = false;
	_maximum_gain = -999.0;
}


/*! \brief Write the analyzed data to a file
*/
void nec_radiation_pattern::write_to_file_aux(ostream& os)
{
	if (false == m_analysis_done)
		throw new nec_exception("Internal Error: Radiation Pattern Analysis not done");
	
	static char  *hpol[4] = { "LINEAR", "RIGHT ", "LEFT  ", " " };
	static char  *gain_type[2] = { "----- POWER GAINS ----- ", "--- DIRECTIVE GAINS ---" };
	static char  *igax[4] = { " MAJOR", " MINOR", " VERTC", " HORIZ" };
	
	int i;
	
	output_helper oh(os,_result_format);
		
	if ( _ifar >= 2) 
	{
		oh.section_start();
		os << "                                 ------ FAR FIELD GROUND PARAMETERS ------" << endl << endl;
		
		if ( _ifar > 3)
		{
			os << endl;
			os << "                                        RADIAL WIRE GROUND SCREEN" << endl;
			os << "                                        "; oh.int_out(5, m_ground.radial_wire_count); os << " WIRES" << endl;
			os << "                                        WIRE LENGTH= "; oh.real_out(8,2, m_ground.radial_wire_length,false); os << " METERS" << endl;
			os << "                                        WIRE RADIUS= "; oh.real_out(10,3, m_ground.radial_wire_radius); os << " METERS" << endl;
		} /* if ( _ifar > 3) */
		
		if ( _ifar != 4 ) 
		{
			std::string hclif;
			if ( (_ifar == 2) || (_ifar == 5) )
				hclif = "LINEAR";
			if ( (_ifar == 3) || (_ifar == 6) )
				hclif= "CIRCLE";
		
		
			os << endl;
			os << "                                        " << hclif << " CLIFF" << endl;
			os << "                                        EDGE DISTANCE= "; oh.real_out(9,2,m_ground.cliff_edge_distance,false); os << " METERS" << endl;
			os << "                                        HEIGHT= "; oh.real_out(8,2,m_ground.cliff_height,false); os << " METERS" << endl;
			os << "                                        SECOND MEDIUM -" << endl;
			os << "                                        RELATIVE DIELECTRIC CONST.= "; oh.real_out(7,3,m_ground.epsr2, false); os << endl;
			os << "                                        CONDUCTIVITY= "; oh.real_out(10,3,m_ground.sig2,false); os << " MHOS" << endl;
		} /* if ( _ifar != 4 ) */
		
	} /* if ( _ifar >= 2) */

	if ( _ifar == 1)
	{
		oh.section_start();
		os << "                             ------- RADIATED FIELDS NEAR GROUND --------" << endl << endl;
		os << "    ------- LOCATION -------     --- E(THETA) ---     ---- E(PHI) ----    --- E(RADIAL) ---" << endl;
		os << "      RHO    PHI        Z           MAG    PHASE         MAG    PHASE        MAG     PHASE" << endl;
		os << "    METERS DEGREES    METERS      VOLTS/M DEGREES      VOLTS/M DEGREES     VOLTS/M  DEGREES" << endl;
	}
	else // _ifar != 1
	{
		oh.section_start();
		os << "                             ---------- RADIATION PATTERNS -----------" << endl << endl;
		
		if ( m_range >= 1.0e-20)
		{
			nec_float exrm = 1.0 / m_range;
			nec_float exra = m_range/ _wavelength;
			exra = -360.0*(exra - floor(exra));
		
			os << "                             RANGE: "; oh.real_out(13,6,m_range); os << " METERS" << endl;
			os << "                             EXP(-JKR)/R: "; oh.real_out(12,5,exrm); os << " AT PHASE: "; oh.real_out(7,2,exra,false); os << " DEGREES" << endl;
		}
		
		int itmp1 = 2 * m_rp_output_format;
		int itmp2 = itmp1+1;
		
		os << " ---- ANGLES -----     "; oh.string_out(23,gain_type[m_rp_ipd]); os << "      ---- POLARIZATION ----   ---- E(THETA) ----    ----- E(PHI) ------" << endl;
		os << "  THETA      PHI      "; oh.string_out(6,igax[itmp1]); os << "  "; oh.string_out(6,igax[itmp2]); os << "   TOTAL       AXIAL      TILT  SENSE   MAGNITUDE    PHASE    MAGNITUDE     PHASE" << endl;
		os << " DEGREES   DEGREES        DB       DB       DB       RATIO   DEGREES            VOLTS/M   DEGREES     VOLTS/M   DEGREES" << endl;
		
	} /* if ( _ifar == 1) */


	i=0;
	nec_float phi = m_phi_start- delta_phi;

	for (int kph = 1; kph <= n_phi; kph++ )
	{
		phi += delta_phi;
		nec_float 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;
		
			/* elliptical polarization */
			if ( _ifar == 1)
			{
				nec_complex e_theta = _e_theta[i];
				nec_complex e_phi = _e_phi[i];
				nec_complex e_r = _e_r[i];
				
				oh.start_record();
				oh.padding(" ");
				oh.real_out(9,2,m_range,false);
				oh.separator(); oh.real_out(7,2,phi,false);
				oh.separator(); oh.real_out(9,2,thet,false);
				oh.separator(); oh.real_out(11,4,abs(e_theta));
				oh.separator(); oh.real_out(7,2,arg_degrees(e_theta),false); 
				oh.separator(); oh.real_out(11,4,abs(e_phi));
				oh.separator(); oh.real_out(7,2,arg_degrees(e_phi),false);
				oh.separator(); oh.real_out(11,4,abs(e_r));
				oh.separator(); oh.real_out(7,2,arg_degrees(e_r),false);
				oh.end_record();
			}
			else
			{
				nec_complex e_theta = _e_theta[i];
				nec_complex e_phi = _e_phi[i];
				
				char* pol_sense = hpol[_polarization_sense_index[i]];
				
				oh.start_record();
				oh.padding(" ");
				oh.real_out(7,2,thet,false); oh.separator(); oh.real_out(9,2,phi,false); oh.separator();
				oh.padding(" ");
				oh.real_out(8,2,_power_gain_vert[i],false);
				oh.separator(); oh.real_out(8,2,_power_gain_horiz[i],false);
				oh.separator(); oh.real_out(8,2,_power_gain_tot[i],false);
				oh.separator(); oh.real_out(11,4,_polarization_axial_ratio[i],false);
				oh.separator(); oh.real_out(9,2,_polarization_tilt[i],false);
				oh.separator(); oh.string_out(6,pol_sense);
				oh.separator(); oh.real_out(11,4,abs(e_theta));
				oh.separator(); oh.real_out(9,2,arg_degrees(e_theta),false);
				oh.separator(); oh.real_out(11,4,abs(e_phi));
				oh.separator(); oh.real_out(9,2,arg_degrees(e_phi),false);
				oh.end_record();
				
				m_plot_card.plot_patterns(thet, phi,
					e_theta, e_phi,
					_power_gain_vert[i], _power_gain_horiz[i], _power_gain_tot[i]);
				
			} /* if ( _ifar != 1) */
		
			i++;
		} /* for( kth = 1; kth <= n_theta; kth++ ) */
	
	} /* for( kph = 1; kph <= n_phi; kph++ ) */

	if ( m_rp_power_average != 0)
	{		
		oh.section_start();
		os << "  AVERAGE POWER GAIN: "; oh.real_out(11,4,_average_power_gain); 
		os << " - SOLID ANGLE USED IN AVERAGING: ("; oh.real_out(7,4,_average_power_solid_angle,false);
		os << ")*PI STERADIANS" << endl;
	}

	if ( m_rp_normalization != 0)
		write_normalized_gain(os);
}


/*! \brief Generate the data for the radiation pattern
*/
void nec_radiation_pattern::analyze(nec_context* m_context)
{
	static nec_float s_impedance = 376.73; // impedance of free space
	
	if (m_analysis_done)
		return;
		
	int pol_sense_index;
	nec_float exrm=0., exra=0., prad, gcon, gcop;
	nec_float phi, thet;
	nec_float tilta, emajr2, eminr2, pol_axial_ratio;
	nec_float dfaz, dfaz2, cdfaz, tstor1=0., tstor2, stilta;
		

	if ((m_context->m_excitation_type == EXCITATION_VOLTAGE) ||
		(m_context->m_excitation_type == EXCITATION_VOLTAGE_DISC) )
	{
		gcop= _wavelength * _wavelength * 2.0 * pi()/(s_impedance * _pinr);
		prad= _pinr- m_context->structure_power_loss- _pnlr;
		gcon= gcop;
		if ( m_context->ipd != 0)
			gcon= gcon* _pinr/ prad;
	}
	else 
	if ( m_context->m_excitation_type == EXCITATION_CURRENT)
	{
		_pinr=394.51* m_context->xpr6* m_context->xpr6* _wavelength* _wavelength;
		gcop= _wavelength* _wavelength*2.* pi()/(s_impedance * _pinr);
		prad= _pinr- m_context->structure_power_loss- _pnlr;
		gcon= gcop;
		if ( m_context->ipd != 0)
			gcon= gcon* _pinr/ prad;
	}
	else
	{
		prad=0.;
		gcon=4.* pi()/(1.0+ m_context->xpr6* m_context->xpr6);
		gcop= gcon;
	}

	if ( m_range >= 1.0e-20)
	{
		exrm=1./ m_range;
		exra= m_range/ _wavelength;
		exra=-360.*( exra- floor( exra));
	}