nec2cpp.cpp

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

/* Translation to C++ by Tim Molteno

	Based on the C port by N. Kyriazis
 	Including some pieces from additional work by
		Jeroen Vreeken <[e-mail]>
	
	Fixed a few bugs in the process.
	
	Using the std vector library in preparation
	for moving to ATLAS for the matrix and vector
	operations.

		Debian Build Instructions
		
	apt-get install atlas3-base atlas3-headers atlas3-base-dev
	apt-get install refblas3-dev lapack3-dev lapack3-doc
	
	For more information on using LAPACK for doing 
	efficient computation, see
		http://seehuhn.de/comp/linear.html
*/

/* Original disclaimer that came with the FORTRAN code */

/******* Translated to the C language by N. Kyriazis  20 Aug 2003 *******/
/*									*/
/* Program NEC(input,tape5=input,output,tape11,tape12,tape13,tape14,	*/
/* tape15,tape16,tape20,tape21)						*/
/*									*/
/* Numerical Electromagnetics Code (NEC2)  developed at Lawrence	*/
/* Livermore lab., Livermore, CA.  (contact G. Burke at 415-422-8414	*/
/* for problems with the NEC code. For problems with the vax implem- 	*/
/* entation, contact J. Breakall at 415-422-8196 or E. Domning at 415 	*/
/* 422-5936) 								*/
/* file created 4/11/80. 						*/
/*									*/
/*                ***********Notice********** 				*/
/* This computer code material was prepared as an account of work 	*/
/* sponsored by the United States government.  Neither the United 	*/
/* States nor the United States Department Of Energy, nor any of 	*/
/* their employees, nor any of their contractors, subcontractors, 	*/
/* or their employees, makes any warranty, express or implied, or	*/
/* assumes any legal liability or responsibility for the accuracy, 	*/
/* completeness or usefulness of any information, apparatus, product 	*/
/* or process disclosed, or represents that its use would not infringe 	*/
/* privately-owned rights. 						*/
/*									*/
/************************************************************************/


#include "nec2cpp.h"
#include "nec_exception.h"
#include <signal.h>
#include <vector>
#include <string>

using namespace std;

#include "nec_context.h"

#ifndef _WIN32
/* Signal handler */
static void sig_handler( int signal );
#endif


/*-------------------------------------------------------------------*/

int nec_main( int argc, char **argv, nec_output_file& s_output );

/*	New main() function

	This places an exception handler around the old main loop to 
	allow errors to be nicely caught!
*/
int main( int argc, char **argv )
{
	nec_output_file s_output;
	
	
	try
	{
		// This is an nec_3vector testharness
		{
			nec_3vector x(1,2,3);
			nec_3vector y(3,5,1);
			
			nec_3vector a = x + y;
			ASSERT(a == nec_3vector(4,7,4));
		
			a = x * y;
			ASSERT(a == nec_3vector(-13,8,-1));
			a = y * x;
			ASSERT(a == nec_3vector(13,-8,1));
			
			a = x + 3.0;
			ASSERT(a == nec_3vector(4,5,6));
			
			nec_float ip = x.dot(y);
			
			ASSERT(ip == 16);
		}
	
		nec_main(argc, argv, s_output);
	}
	catch (const char* message)
	{
		nec_error_mode nem(s_output);
		s_output.line("NEC++ Runtime Error: ");
		s_output.line(message);
		exit(1);
	}
	catch (nec_exception* nex)
	{
		nec_error_mode nem(s_output);
		s_output.line("NEC++ Runtime Error: ");
		s_output.line(nex->get_message().c_str());
		exit(1);
	}
}

#include "c_geometry.h"
void benchmark();

void benchmark()
{
	cout << "The nec2++ benchmark." << endl;
	cout << "nec2++ version " nec_version << endl << endl;
	
	nec_float bench = nec_context::benchmark();
	cout << "Your computer's score is: " << bench << " NEC's" << endl;
}

int readmn(FILE* input_fp, FILE* output_fp, char *gm, int *i1, int *i2, int *i3, int *i4, nec_float *f1,
	nec_float *f2, nec_float *f3, nec_float *f4, nec_float *f5, nec_float *f6);


#include "XGetopt.h"

int nec_main( int argc, char **argv, nec_output_file& s_output )
{
	nec_output_flags s_output_flags;
	FILE *input_fp=NULL;
	FILE *output_fp=NULL;
	
	string input_filename, output_filename;
	
	char ain[3], line_buf[81];
	
	/* input card mnemonic list */
	/* "XT" stands for "exit", added for testing */
	#define CMD_NUM  20
	char *atst[CMD_NUM] =
	{
		"FR", "LD", "GN", "EX", "NT", "TL", \
		"XQ", "GD", "RP", "NX", "PT", "KH", \
		"NE", "NH", "PQ", "EK", "CP", "PL", \
		"EN", "WG"
	};
	
	int itmp3, itmp2, itmp4;
	
	int ain_num;    /* ain mnemonic as a number */
	
	nec_float tmp1, tmp2, tmp3, tmp4, tmp5, tmp6;
	nec_float ex_timer;
	
	/* getopt() variables */
	extern char *optarg;
	int option;
	
#ifndef _WIN32
	/*** signal handler related code ***/
	/* new and old actions for sigaction() */
	struct sigaction sa_new, sa_old;
	
	
	/* initialize new actions */
	sa_new.sa_handler = sig_handler;
	sigemptyset( &sa_new.sa_mask );
	sa_new.sa_flags = 0;
	
	/* register function to handle signals */
	sigaction( SIGINT,  &sa_new, &sa_old );
	sigaction( SIGSEGV, &sa_new, 0 );
	sigaction( SIGFPE,  &sa_new, 0 );
	sigaction( SIGTERM, &sa_new, 0 );
	sigaction( SIGABRT, &sa_new, 0 );
#endif
	
	/*** command line arguments handler ***/
	if ( argc == 1 )
	{
		usage();
		exit(-1);
	}
	
	bool results_to_stdout = false;
	bool results_in_csv_format = false;
	
	/* process command line options */
	while( (option = XGetopt(argc, argv, "i:o:hvscgb") ) != -1 )
	{
		switch( option )
		{
		case 'i' : /* specify input file name */
			input_filename = optarg;
			break;
		
		case 'o' : /* specify output file name */
			output_filename = optarg;
			break;
		
		case 'g': /* return only the maximum gain to stdout */
			s_output_flags.set_gain_only(true);
			break;

		case 's': /* return output to stdout */
			results_to_stdout = true;
			break;

		case 'c': /* use CVS result data */
			results_in_csv_format = true;
			break;
	
		case 'h' : /* print usage and exit */
			usage();
			exit(0);
		
		case 'v' : /* print nec2++ version */
#ifdef _MSC_VER
			cout << ( "nec2++ " nec_version ) << " compiler: " << _MSC_VER << endl;
#else
			cout << ( "nec2++ " nec_version ) << (" compiler: " __VERSION__) << endl;
#endif
			exit(0);

		case 'b' : /* Run benchmark */
			benchmark();
			exit(0);
		
		default: /* print usage and exit */
			usage();
			exit(-1);
		
		}	
	}

	/*** open input file ***/
	if ( (input_fp = fopen(input_filename.c_str(), "r")) == NULL )
	{
		string mesg = "nec2++: "  + input_filename;
		
		perror( mesg.c_str() );
		exit(-1);
	}

	/* make an output file name if not */
	/* specified by user on invocation */
	if ( output_filename  == "" )
	{
		/* strip the input file name extension if there is one */
		output_filename = input_filename.substr(0, input_filename.find(".",0)) + ".out";
	}
	
	/* open output file */
	if ( (output_fp = fopen(output_filename.c_str(), "w")) == NULL )
	{
		string mesg = "nec2++: "  + output_filename;
		
		perror( mesg.c_str() );
		exit(-1);
	}
	
	s_output.set_file(output_fp);
    
	secnds( &ex_timer );

	// allocate a new nec_context;
	nec_context s_context;
	s_context.set_output(s_output, s_output_flags);
	s_context.set_results_stdout(results_to_stdout);
	if (results_in_csv_format)
		s_context.set_results_format(RESULT_FORMAT_CSV);
	s_context.initialize();
    
  /* main execution loop, exits at various points */
  /* depending on error conditions or end of jobs */
	while( true )
	{
	
		s_output.end_section();
		s_output.set_indent(31);
		s_output.line(" __________________________________________");
		s_output.line("|                                          |");
		s_output.line("| NUMERICAL ELECTROMAGNETICS CODE (nec2++) |");
		s_output.line("| Translated to 'C++' in Double Precision  |");
		s_output.line("|        Version " nec_version "        |");
		s_output.line("|__________________________________________|");
	
		/* read a line from input file */
		if ( load_line(line_buf, input_fp) == EOF )
			throw new nec_exception("Error reading input file.");
		
		/* separate card's id mnemonic */
		strncpy( ain, line_buf, 2 );
		ain[2] = '\0';
	
		/* If its an "XT" card, exit (used for debugging) */
		if ( strcmp(ain, "XT") == 0 )
		{
			nec_error_mode em(s_output);
			s_output.end_section();
			s_output.line("nec2++: Exiting after an \"XT\" command in main()");
			exit(0);
		}
	
		/* if its a "cm" or "ce" card start reading comments */
		if ( (strcmp(ain, "CM") == 0) ||
			(strcmp(ain, "CE") == 0) )
		{
			s_output.end_section();
			s_output.set_indent(31);
			s_output.line("---------------- COMMENTS ----------------");
			s_output.line(&line_buf[2]);
			while( strcmp(ain, "CM") == 0 )
			{
				/* read a line from input file */
				if ( load_line(line_buf, input_fp) == EOF )
					throw new nec_exception("Error reading input file (comments not terminated?)");
			
				/* separate card's id mnemonic */
				strncpy( ain, line_buf, 2 );
				ain[2] = '\0';
			
				/* write comment to output file */
				s_output.line(&line_buf[2]);
			}
		
			/* no "ce" card at end of comments */
			if ( strcmp(ain, "CE") != 0 )
			{
				throw new nec_exception("ERROR: INCORRECT LABEL FOR A COMMENT CARD");
			}
		}
		else
		{
			rewind( input_fp );
		}
	
		/* initializations etc from original fortran code */
		int data_card_count=0;
		
		/* set up geometry data in subroutine parse_geometry */
		c_geometry* geo = s_context.get_geometry();
		geo->parse_geometry(&s_context, input_fp);
		
		s_context.calc_prepare();
		s_output.end_section();
		
		/*
			Main input section, exits at various points
			depending on error conditions or end of job.
			This is called the card input loop.
		*/
				
		bool next_job = false; /* start next job (next structure) flag */
		while ( ! next_job )
		{
			int itmp1;
			
			/* main input section - standard read statement - jumps */
			/* to appropriate section for specific parameter set up */
			int parameter_count = readmn(input_fp, output_fp, ain, &itmp1, &itmp2, &itmp3, &itmp4,
				&tmp1, &tmp2, &tmp3, &tmp4, &tmp5, &tmp6 );
			
			/* If its an "XT" card, exit */
			if ( strcmp(ain, "XT" ) == 0 )
			{
				nec_error_mode em(s_output);
				s_output.endl();
				s_output.line("nec2++: Exiting after an \"XT\" command in main()" );
				exit(0);
			}
			
			data_card_count++;
			fprintf( output_fp,
				"\n*****  DATA CARD N0. %3d "
				"%s %3d %5d %5d %5d %12.5E %12.5E %12.5E %12.5E %12.5E %12.5E",
				data_card_count, ain, itmp1, itmp2, itmp3, itmp4,
				tmp1, tmp2, tmp3, tmp4, tmp5, tmp6 );
			
			/* identify card id mnemonic (except "ce" and "cm") */
			for( ain_num = 0; ain_num < CMD_NUM; ain_num++ )
				if ( strncmp( ain, atst[ain_num], 2) == 0 )
					break;