c_geometry.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 "c_geometry.h"

#include "nec_context.h"
#include "nec_exception.h"

c_geometry::c_geometry()
{
	n = 0;
	np = 0; 	// n is the number of segments
	
	m = 0;
	mp = 0;		// m is the number of patches?
	
	m_ipsym = 0;
	
	n_plus_m = 0;
	n_plus_2m = 0;
	n_plus_3m = 0;
	
	jsno = 0;
	nscon = 0;
	maxcon = 0;
	
	m_context = NULL;
	m_output = NULL;
}

void c_geometry::set_context(nec_context* in_context)
{
	m_context = in_context;
	m_output = &m_context->m_output;
}

/*! \brief Get a segment number for a specified tag.
	\param in_tag The tag
	\param in_m The mth segment with the specified tag will be returned.
	\return The segment number of the mth segment having the
	tag number in_tag.  if in_tag=0 segment number m is returned.
*/
int c_geometry::get_segment_number( int in_tag, int in_m)
{
	ASSERT(in_tag >= 0);
	ASSERT(in_m >= 0);
	
	
	if (in_m <= 0)
	{
		throw new nec_exception("CHECK DATA, PARAMETER SPECIFYING SEGMENT POSITION IN A GROUP OF EQUAL TAGS MUST NOT BE ZERO" );
	}
	
	if ( 0 == in_tag)
	{
		return( in_m );
	}
	
	int tag_seg_count=0;
	for (int i = 0; i < n; i++ )
	{
		if ( segment_tags[i] == in_tag )
		{
			tag_seg_count++;
			if ( tag_seg_count == in_m)
			{
				return( i+1 );
			}
		}
	}
	
	throw new nec_exception("NO SEGMENT HAS AN ITAG OF ", in_tag);
	
	return 0;
}


/*
	parse_geometry is the main routine for input of geometry data.

	ANTLR Grammar for geometry input
	
		geometry
			:	(geometry_element {transform}*)+
				geometry_end
			;
		
		geometry_element
			:	wire
			|	arc
			|	helix
			|	patch
			|	multi_patch
			;
		
		transform
			:	scale
			|	reflect
			|	rotate
			|	move
			;
		
		wire = "GW" int int real real real real real real {taper}
		taper = "GC" xxx
		arc = "GA"
		helix = "GH"
		
		patch = "SP"
		multi_patch = "SM" {patch_data}+
		patch_data = "SC"
		
		scale = "GS"
		reflect = "GX"
		rotate = "GR"
		move = "GM"
		
		geometry_end = "GE" .+
		
		int:
			[0-9]+;
			
		real:
			int:{.int};
		
*/

#include "c_plot_card.h"
	

void c_geometry::parse_geometry(nec_context* in_context, FILE* input_fp )
{
	char gm[3];
	char ifx[2] = {'*', 'X'}, ify[2]={'*','Y'}, ifz[2]={'*','Z'};
	char ipt[4] = { 'P', 'R', 'T', 'Q' };
	
	/* input card mnemonic list */
	/* "XT" stands for "exit", added for testing */
	#define GM_NUM  12
	char *atst[GM_NUM] =
	{
		"GW", "GX", "GR", "GS", "GE", "GM", \
		"SP", "SM", "GA", "SC", "GH", "GF"
	};
	
	bool print_structure_spec = true;
	
	int nwire, isct, i1, i2, iy, iz;
	int ix; 
	int card_int_1, card_int_2; /* The two integer parameters from the geometry card */
	nec_float rad, xs1, xs2, ys1, ys2, zs1, zs2, x4=0, y4=0, z4=0;
	nec_float x3=0, y3=0, z3=0, xw1, xw2, yw1, yw2, zw1, zw2;
	nec_float dummy;
	
	m_ipsym=0;
	nwire=0;
	n=0;
	np=0;
	m=0;
	mp=0;
	isct=0;
	
	/* read geometry data card and branch to */
	/* section for operation requested */
	do
	{
		read_geometry_card(input_fp, gm, &card_int_1, &card_int_2, &xw1, &yw1, &zw1, &xw2, &yw2, &zw2, &rad);
	
		/* identify card id mnemonic */
		int gm_num;
		for( gm_num = 0; gm_num < GM_NUM; gm_num++ )
			if ( strncmp( gm, atst[gm_num], 2) == 0 )
				break;

		if ( print_structure_spec )
		{
			m_output->end_section();
			m_output->set_indent(32);
			m_output->line("-------- STRUCTURE SPECIFICATION --------");
			m_output->line("COORDINATES MUST BE INPUT IN" );
			m_output->line("METERS OR BE SCALED TO METERS" );
			m_output->line("BEFORE STRUCTURE INPUT IS ENDED" );
			m_output->set_indent(0);
		
			m_output->line("  WIRE                                                                                 SEG FIRST  LAST  TAG");
			m_output->line("   No:        X1         Y1         Z1         X2         Y2         Z2       RADIUS   No:   SEG   SEG  No:" );
		
			print_structure_spec = false;
		}

		if ( gm_num != 10 )
			isct=0;

		switch( gm_num )
		{
		case 0:
		/* "gw" card, generate segment data for straight wire.
			GW		STRAIGHT WIRE, ENDS 1,2
				card_int_1- TAG NO.
				card_int_2- NO. SEGMENTS
				xw1- X1
				F2- Y1
				F3- Z1
				F4- X2
				F5- Y2
				F6- Z2
				F7- WIRE RAD., 0=USE GC FOR TAPERED WIRE
		*/
		{
			int wire_segment_count = card_int_2;
			int wire_tag = card_int_1;
			
			nwire++;
		
			// output some wire diagnostics.
			m_output->nec_printf( "\n"
				" %5d  %10.4f %10.4f %10.4f %10.4f"
				" %10.4f %10.4f %10.4f %5d %5d %5d %4d",
				nwire, xw1, yw1, zw1, xw2, yw2, zw2, rad, wire_segment_count, n+1, n + wire_segment_count, wire_tag );
		
			if ( rad != 0)	// rad == 0 implies a tapered wire
			{
				xs1 = 1.0;
				ys1 = 1.0;
			}
			else
			{
				read_geometry_card(input_fp, gm, &ix, &iy, &xs1, &ys1, &zs1,
					&dummy, &dummy, &dummy, &dummy);
			
				if ( strcmp(gm, "GC" ) != 0 )
				{
					throw new nec_exception("GEOMETRY DATA CARD ERROR" );
				}
			
				m_output->nec_printf(
					"\n  ABOVE WIRE IS TAPERED.  SEGMENT LENGTH RATIO: %9.5f\n"
					"                                 "
					"RADIUS FROM: %9.5f TO: %9.5f", xs1, ys1, zs1 );
			
				if ( (ys1 == 0) || (zs1 == 0) )
				{
					throw new nec_exception("GEOMETRY DATA CARD ERROR" );
				}
			
				rad= ys1;
				ys1= pow( (zs1/ys1), (1./(wire_segment_count-1.)) );
			}
		
			wire(wire_tag, wire_segment_count, xw1, yw1, zw1, xw2, yw2, zw2, rad, xs1, ys1);
		}
		continue;

		/* reflect structure along x,y, or z */
		/* axes or rotate to form cylinder.  */
		case 1: /* "gx" card */
			iy= card_int_2/10;
			iz= card_int_2- iy*10;
			ix= iy/10;
			iy= iy- ix*10;
		
			if ( ix != 0)
			ix=1;
			if ( iy != 0)
			iy=1;
			if ( iz != 0)
			iz=1;
		
			m_output->nec_printf(
				"\n  STRUCTURE REFLECTED ALONG THE AXES %c %c %c"
				" - TAGS INCREMENTED BY %d\n",
				ifx[ix], ify[iy], ifz[iz], card_int_1 );
		
			reflect( ix, iy, iz, card_int_1, card_int_2);
		
			continue;
	
		case 2: /* "gr" card */
		
			m_output->nec_printf(
				"\n  STRUCTURE ROTATED ABOUT Z-AXIS %d TIMES"
				" - LABELS INCREMENTED BY %d\n", card_int_2, card_int_1 );
		
			ix=-1;
			iz = 0;
			reflect( ix, iy, iz, card_int_1, card_int_2);
		
			continue;
	
		case 3: /* "gs" card, scale structure dimensions by factor xw1. */
		
			if ( n > 0)
			{
				for(int i = 0; i < n; i++ )
				{
					x[i]= x[i]* xw1;
					y[i]= y[i]* xw1;
					z[i]= z[i]* xw1;
					x2[i]= x2[i]* xw1;
					y2[i]= y2[i]* xw1;
					z2[i]= z2[i]* xw1;
					segment_radius[i]= segment_radius[i]* xw1;
				}
			} /* if ( n >= n2) */
		
			if ( m > 0)
			{
				yw1= xw1* xw1;
				for(int i = 0; i < m; i++ )
				{
					px[i]= px[i]* xw1;
					py[i]= py[i]* xw1;
					pz[i]= pz[i]* xw1;
					pbi[i]= pbi[i]* yw1;
				}
			} /* if ( m >= m2) */
		
			m_output->nec_printf(
				"\n     STRUCTURE SCALED BY FACTOR: %10.5f", xw1 );
		
			continue;

		case 4: /* "ge" card, terminate structure geometry input. */

			geometry_complete(in_context, card_int_1, card_int_2);
			
			return;

		/* "gm" card, move structure or reproduce */
		/* original structure in new positions.   */
		case 5:
		
			m_output->nec_printf(
				"\n     THE STRUCTURE HAS BEEN MOVED, MOVE DATA CARD IS:\n"
				"   %3d %5d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f",
				card_int_1, card_int_2, xw1, yw1, zw1, xw2, yw2, zw2, rad );
		
			xw1= degrees_to_rad(xw1);
			yw1= degrees_to_rad(yw1);
			zw1= degrees_to_rad(zw1);
		
			move( xw1, yw1, zw1, xw2, yw2, zw2, (int)( rad+.5), card_int_2, card_int_1);
			continue;
		
		case 6: /* "sp" card, generate single new patch */
		
			i1= m+1;
			card_int_2++;
		
			if ( card_int_1 != 0)
			{
				throw new nec_exception("PATCH DATA ERROR" );
			}
		
			m_output->nec_printf( "\n"
				" %5d%c %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f",
				i1, ipt[card_int_2-1], xw1, yw1, zw1, xw2, yw2, zw2 );
		
			if ( (card_int_2 == 2) || (card_int_2 == 4) )
				isct=1;
		
			if ( card_int_2 > 1)
			{
				// read another geometry card for the rest of the patch data
				read_geometry_card(input_fp, gm, &ix, &iy, &x3, &y3, &z3, &x4, &y4, &z4, &dummy);
			
				if ( (card_int_2 == 2) || (card_int_1 > 0) )
				{
					x4= xw1+ x3- xw2;
					y4= yw1+ y3- yw2;
					z4= zw1+ z3- zw2;
				}
			
				m_output->nec_printf( "\n"
					"      %11.5f %11.5f %11.5f %11.5f %11.5f %11.5f",
					x3, y3, z3, x4, y4, z4 );
			
				if ( strcmp(gm, "SC") != 0 )
				{
					throw new nec_exception("PATCH DATA ERROR" );
				}
			} /* if ( card_int_2 > 1) */
			else
			{
				xw2= degrees_to_rad(xw2);
				yw2= degrees_to_rad(yw2);
			}
		
			patch( card_int_1, card_int_2, xw1, yw1, zw1, xw2, yw2, zw2, x3, y3, z3, x4, y4, z4);
			continue;
	
		case 7: /* "sm" card, generate multiple-patch surface */
		
			i1= m+1;
			m_output->nec_printf( "\n"
				" %5d%c %10.5f %11.5f %11.5f %11.5f %11.5f %11.5f"
				"     SURFACE - %d BY %d PATCHES",
				i1, ipt[1], xw1, yw1, zw1, xw2, yw2, zw2, card_int_1, card_int_2 );
		
			if ( (card_int_1 < 1) || (card_int_2 < 1) )
			{
				throw new nec_exception("PATCH DATA ERROR" );
			}
		
			read_geometry_card(input_fp, gm, &ix, &iy, &x3, &y3, &z3, &x4, &y4, &z4, &dummy);
		
			if ( (card_int_2 == 2) || (card_int_1 > 0) )
			{
				x4= xw1+ x3- xw2;
				y4= yw1+ y3- yw2;
				z4= zw1+ z3- zw2;
			}
		
			m_output->nec_printf( "\n"
				"      %11.5f %11.5f %11.5f %11.5f %11.5f %11.5f",
				x3, y3, z3, x4, y4, z4 );
		
			if ( strcmp(gm, "SC" ) != 0 )
			{
				throw new nec_exception("PATCH DATA ERROR" );
			}
		
			patch( card_int_1, card_int_2, xw1, yw1, zw1, xw2, yw2, zw2, x3, y3, z3, x4, y4, z4);
			continue;
		
		case 8: /* "ga" card, generate segment data for wire arc */
		
			nwire++;
			i1= n+1;
			i2= n+ card_int_2;
		
			m_output->nec_printf( "\n"
				" %5d  ARC RADIUS: %9.5f  FROM: %8.3f TO: %8.3f DEGREES"
				"       %11.5f %5d %5d %5d %4d",
				nwire, xw1, yw1, zw1, xw2, card_int_2, i1, i2, card_int_1 );
		
			arc( card_int_1, card_int_2, xw1, yw1, zw1, xw2);
			continue;
	
		case 9: /* "sc" card */
		
			if ( isct == 0)
			{
				throw new nec_exception("PATCH DATA ERROR" );
			}
		
			i1= m+1;
			card_int_2++;
		
			if ( (card_int_1 != 0) || ((card_int_2 != 2) && (card_int_2 != 4)) )
			{
				throw new nec_exception("PATCH DATA ERROR" );
			}
		
			xs1= x4;
			ys1= y4;
			zs1= z4;
			xs2= x3;
			ys2= y3;
			zs2= z3;
			x3= xw1;
			y3= yw1;
			z3= zw1;
		
			if ( card_int_2 == 4)
			{
				x4= xw2;
				y4= yw2;
				z4= zw2;
			}
		
			xw1= xs1;
			yw1= ys1;
			zw1= zs1;
			xw2= xs2;
			yw2= ys2;
			zw2= zs2;
		
			if ( card_int_2 != 4)
			{
				x4= xw1+ x3- xw2;
				y4= yw1+ y3- yw2;
				z4= zw1+ z3- zw2;
			}
		
			m_output->nec_printf( "\n"
				" %5d%c %10.5f %11.5f %11.5f %11.5f %11.5f %11.5f",
				i1, ipt[card_int_2-1], xw1, yw1, zw1, xw2, yw2, zw2 );
		
			m_output->nec_printf( "\n"
				"      %11.5f %11.5f %11.5f  %11.5f %11.5f %11.5f",
				x3, y3, z3, x4, y4, z4 );
		
			patch( card_int_1, card_int_2, xw1, yw1, zw1, xw2, yw2, zw2, x3, y3, z3, x4, y4, z4);
		
			continue;
	
		case 10: /* "gh" card, generate helix */
		
			nwire++;
			i1= n+1;
			i2= n+ card_int_2;
		
			m_output->nec_printf( "\n"
				" %5d HELIX STRUCTURE - SPACING OF TURNS: %8.3f AXIAL"
				" LENGTH: %8.3f  %8.3f %5d %5d %5d %4d\n      "
				" RADIUS X1:%8.3f Y1:%8.3f X2:%8.3f Y2:%8.3f ",
				nwire, xw1, yw1, rad, card_int_2, i1, i2, card_int_1, zw1, xw2, yw2, zw2 );
		
			helix( xw1, yw1, zw1, xw2, yw2, zw2, rad, card_int_2, card_int_1);
		
			continue;

		case 11: /* "gf" card, not supported */
			throw new nec_exception("NGF solution option not supported");
	
		default: /* error message */
		
			m_output->nec_printf( "\n  GEOMETRY DATA CARD ERROR" );
			m_output->nec_printf( "\n"
				" %2s %3d %5d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f",
				gm, card_int_1, card_int_2, xw1, yw1, zw1, xw2, yw2, zw2, rad );
		
			throw new nec_exception("GEOMETRY DATA CARD ERROR");
    } /* switch( gm_num ) */

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

/**
	We have finished with the geometry description, now connect 
	things up.
*/
void c_geometry::geometry_complete(nec_context* in_context, int card_int_1, int card_int_2)
{
	if (0 == np + mp)
		throw new nec_exception("Geometry has no wires or patches.");
		
	/* TCAM: The following does not make sense for the semantics of the plot
	card. I have left this in, I hope someone will tell me why it is required.
	*/
	if (card_int_2 != 0)
		in_context->plot_card.set_plot_real_imag_currents();
		
	// now proceed and complete the geometry setup...
	
	connect_segments( card_int_1);

	if ( n != 0)
	{
		/* Allocate wire buffers */
		segment_length.resize(n);
		sab.resize(n);
		cab.resize(n);
		salp.resize(n);
	
		m_output->nec_printf( "\n\n\n"
			"                              "
			" ---------- SEGMENTATION DATA ----------\n"
			"                                       "
			" COORDINATES IN METERS\n"
			"                           "
			" I+ AND I- INDICATE THE SEGMENTS BEFORE AND AFTER I\n" );
	
		m_output->nec_printf( "\n"
			"   SEG    COORDINATES OF SEGM CENTER     SEGM    ORIENTATION"
			" ANGLES    WIRE    CONNECTION DATA   TAG\n"
			"   No:       X         Y         Z      LENGTH     ALPHA     "
			" BETA    RADIUS    I-     I    I+   NO:" );
	
		for(int i = 0; i < n; i++ )
		{
			nec_float xw1= x2[i]- x[i];
			nec_float yw1= y2[i]- y[i];
			nec_float zw1= z2[i]- z[i];
			x[i]=( x[i]+ x2[i])*.5;
			y[i]=( y[i]+ y2[i])*.5;
			z[i]=( z[i]+ z2[i])*.5;
			
			nec_float xw2= xw1* xw1+ yw1* yw1+ zw1* zw1;
			nec_float yw2= sqrt( xw2);
			yw2=( xw2/ yw2+ yw2)*.5;
			segment_length[i]= yw2;
			cab[i]= xw1/ yw2;
			sab[i]= yw1/ yw2;
			xw2= zw1/ yw2;
	
			if ( xw2 > 1.)
				xw2=1.;
			if ( xw2 < -1.)
				xw2=-1.;
	
			salp[i]= xw2;
			xw2= rad_to_degrees(asin( xw2));
			yw2= rad_to_degrees(atan2( yw1, xw1));
	
			m_output->nec_printf( "\n"
				" %5d %9.4f %9.4f %9.4f %9.4f"
				" %9.4f %9.4f %9.4f %5d %5d %5d %5d",
				i+1, x[i], y[i], z[i], segment_length[i], xw2, yw2,
				segment_radius[i], icon1[i], i+1, icon2[i], segment_tags[i] );
	
			in_context->plot_card.plot_segments(i,x,y,z,segment_length,xw2,yw2,segment_radius,icon1,icon2);
	
			if ( (segment_length[i] <= 1.e-20) || (segment_radius[i] <= 0.) )
			{
				throw new nec_exception("SEGMENT DATA ERROR" );
			}
		} /* for( i = 0; i < n; i++ ) */
	} /* if ( n != 0) */

	if ( m != 0)