net_solve.cpp

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

/*
	Copyright (C) 2004  Timothy C.A. Molteno
	
	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

	This class is not fully functional yet. It will perform
	the network solution when done!
*/
#include "matrix_algebra.h"
#include "electromag.h"
#include "nec_ground.h"
#include "c_geometry.h"

#include <stdio.h>

#include "nec_results.h"
extern nec_results s_results;

#include "nec_output.h"
extern nec_output_file s_output;
extern nec_output_flags s_output_flags;
extern nec_float wavelength;
extern int nop;
extern complex_array symmetry_array;


/* subroutine netwk solves for structure currents for a given */
/* excitation including the effect of non-radiating networks if */
/* present. */
class c_network
{
public:
	c_network();
	
	void net_solve(
		complex_array& cm,
		nec_complex *cmb,
		nec_complex *cmc,
		nec_complex *cmd,
		int_array& ip,
		complex_array& einc );
	
private:
	
	FILE *output_fp;
	
	c_geometry geometry;
	nec_ground ground;
	/* common  /crnt/
	
	*/
	real_array air, aii; // coefficients of the constant terms in the current interpolation functions for the current vector 
	real_array bir, bii; // coefficients of the sine terms in the current interpolation functions
	real_array cir, cii; // coefficients of the cosine terms in the current interpolation functions
	complex_array current_vector;	// the current vector
	
	/* common  /vsorc/ */
	int_array ivqd, source_segment_array, iqds;
	int nvqd, voltage_source_count, nqds;
	complex_array vqd, vqds, source_voltage_array;
	
	/* common  /netcx/ */
	int masym, neq, npeq, neq2, network_count, ntsol, nprint;
	int_array iseg1, iseg2, ntyp;
	real_array x11r, x11i, x12r;
	real_array x12i, x22r, x22i;
	nec_float input_power, network_power_loss;
	nec_complex zped;
};

void c_network::net_solve( complex_array& cm, nec_complex *cmb,
    nec_complex *cmc, nec_complex *cmd, int_array& ip,
    complex_array& einc )
{	
	/* Network buffers */
	int_array ipnt, nteqa, ntsca;
	complex_array vsrc, rhs, cmn, rhnt, rhnx;
	
	bool jump1, jump2;
	
	int nteq=0, ntsc=0, nseg2, irow2=0;
	int neqz2, neqt, irow1=0, i, nseg1, isc1=0, isc2=0;
	nec_float asmx, asa, y11r, y11i, y12r, y12i, y22r, y22i;
	nec_complex ymit, vlt, cux;
	
	neqz2= neq2;
	if ( neqz2 == 0)
		neqz2=1;
	
	input_power = 0.0;
	network_power_loss = 0.0;
	neqt= neq+ neq2;
	
	int ndimn = (2*network_count + voltage_source_count);
	
	/* Allocate network buffers */
	if ( network_count > 0 )
	{
		rhs.resize( geometry.n_plus_3m ); // this should probably be ndimn!
	
		rhnt.resize( ndimn );
		rhnx.resize( ndimn);
		cmn.resize( ndimn * ndimn );
	
		ntsca.resize( ndimn );
		nteqa.resize( ndimn );
		ipnt.resize( ndimn );
	
		vsrc.resize( voltage_source_count );
	}
	
	if ( ntsol == 0)
	{
		/* compute relative matrix asymmetry */
		if ( masym != 0)
		{
			irow1=0;
			for( i = 0; i < network_count; i++ )
			{
				nseg1= iseg1[i];
				for( isc1 = 0; isc1 < 2; isc1++ )
				{
					if ( irow1 == 0)
					{
						ipnt[irow1]= nseg1;
						nseg1= iseg2[i];
						irow1++;
						continue;
					}
			
					int j = 0;
					for( j = 0; j < irow1; j++ )
						if ( nseg1 == ipnt[j])
							break;
			
					if ( j == irow1 )
					{
						ipnt[irow1]= nseg1;
						irow1++;
					}
			
					nseg1= iseg2[i];
				} /* for( isc1 = 0; isc1 < 2; isc1++ ) */
			} /* for( i = 0; i < network_count; i++ ) */
	
			ASSERT(voltage_source_count >= 0);
			for( i = 0; i < voltage_source_count; i++ )
			{
				nseg1= source_segment_array[i];
				if ( irow1 == 0)
				{
					ipnt[irow1]= nseg1;
					irow1++;
					continue;
				}
			
				int j = 0;
				for( j = 0; j < irow1; j++ )
					if ( nseg1 == ipnt[j])
						break;
			
				if ( j == irow1 )
				{
					ipnt[irow1]= nseg1;
					irow1++;
				}
			} /* for( i = 0; i < voltage_source_count; i++ ) */
		
			if ( irow1 >= 2)
			{
				for( i = 0; i < irow1; i++ )
				{
					isc1 = ipnt[i]-1;
					asmx= geometry.segment_length[isc1];
				
					for (int j = 0; j < neqt; j++ )
						rhs[j] = cplx_00();
				
					rhs[isc1] = cplx_10();
					solves( cm, ip, rhs, neq, 1, geometry.np, geometry.n, geometry.mp, geometry.m, nop, symmetry_array);
					geometry.get_current_coefficients(wavelength, rhs, air, aii, bir, bii, cir, cii, vqds, nqds, iqds);
				
					for (int j = 0; j < irow1; j++ )
					{
						isc1= ipnt[j]-1;
						cmn[j+i*ndimn]= rhs[isc1]/ asmx;
					}
				} /* for( i = 0; i < irow1; i++ ) */
			
				asmx=0.0;
				asa=0.0;
			
				for( i = 1; i < irow1; i++ )
				{
					for (int j = 0; j < i; j++ )
					{
						cux = cmn[i+j*ndimn];
						nec_float pwr= abs(( cux- cmn[j+i*ndimn])/ cux);
						asa += pwr* pwr;
				
						if ( pwr >= asmx)
						{
							asmx= pwr;
							nteq= ipnt[i];
							ntsc= ipnt[j];
						}
					} /* for( j = 0; j < i; j++ ) */	
				} /* for( i = 1; i < irow1; i++ ) */
	
				asa= sqrt( asa*2./ (nec_float)( irow1*( irow1-1)));
				fprintf( output_fp, "\n\n"
					"   MAXIMUM RELATIVE ASYMMETRY OF THE DRIVING POINT ADMITTANCE\n"
					"   MATRIX IS %10.3E FOR SEGMENTS %d AND %d\n"
					"   RMS RELATIVE ASYMMETRY IS %10.3E",
					asmx, nteq, ntsc, asa );
			} /* if ( irow1 >= 2) */
		} /* if ( masym != 0) */
	
		/* solution of network equations */
		if ( network_count != 0)
		{
			// zero the cmn array, and the rhnx array
			cmn.fill(cplx_00());
			rhnx.fill(cplx_00());
			
/*			for( i = 0; i < ndimn; i++ )
			{
				rhnx[i]=cplx_00();
				for (int j = 0; j < ndimn; j++ )
					cmn[j+i*ndimn]=cplx_00();
			}
*/			
			nteq=0;
			ntsc=0;
	
			/*	sort network and source data and
				assign equation numbers to segments */
			for (int j = 0; j < network_count; j++ )
			{
				nseg1= iseg1[j];
				nseg2= iseg2[j];
			
				if ( ntyp[j] <= 1)
				{
					y11r= x11r[j];
					y11i= x11i[j];
					y12r= x12r[j];
					y12i= x12i[j];
					y22r= x22r[j];
					y22i= x22i[j];
				}
				else
				{
					y22r= two_pi() * x11i[j]/ wavelength;
					y12r=0.;
					y12i=1./( x11r[j]* sin( y22r));
					y11r= x12r[j];
					y11i=- y12i* cos( y22r);
					y22r= x22r[j];
					y22i= y11i+ x22i[j];
					y11i= y11i+ x12i[j];
				
					if ( ntyp[j] != 2)
					{
						y12r=- y12r;
						y12i=- y12i;
					}
				} /* if ( ntyp[j] <= 1) */
		
				jump1 = false;
				for( i = 0; i < voltage_source_count; i++ )
				{
					if ( nseg1 == source_segment_array[i])
					{
						isc1 = i;
						jump1 = true;
						break;
					}
				}
		
				jump2 = false;
				if ( ! jump1 )
				{
					isc1=-1;
			
					for( i = 0; i < nteq; i++ )
					{
						if ( nseg1 == nteqa[i])
						{
							irow1 = i;
							jump2 = true;
							break;
						}
					}
			
					if ( ! jump2 )
					{
						irow1= nteq;
						nteqa[nteq]= nseg1;
						nteq++;
					}
				} /* if ( ! jump1 ) */
				else
				{
					for( i = 0; i < ntsc; i++ )
					{
						if ( nseg1 == ntsca[i])
						{
							irow1 = ndimn- (i+1);
							jump2 = true;
							break;