c_evlcom.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
*/
#include "c_evlcom.h"

#include "matrix_algebra.h" // for test()
#include "nec_exception.h"

using namespace std;

/* 
	What on earth is this NM thing?
*/
#define NM	131072
#define MAXH	20
#define CRIT	1.0E-4
#define	PTP	.6283185308
#define NTS	4

/* compute integration parameter xlam=lambda from parameter t. */
void c_evlcom::lambda( nec_float t, nec_complex *xlam, nec_complex *dxlam ) const
{
	*dxlam = m_contour_b - m_contour_a;
	*xlam = m_contour_a + *dxlam*t;
}


/*! \brief gshank integrates the 6 Sommerfeld integrals from start to 
	infinity (until convergence) in lambda.  At the break point, bk, 
	the step increment may be changed from dela to delb.  Shank's 
	algorithm to accelerate convergence of a slowly converging series
	is used. */
void c_evlcom::gshank( nec_complex start, nec_complex dela, complex_array& sum,
	int nans, complex_array& seed, int ibk, nec_complex bk, nec_complex delb )
{
	bool brk = false;
	int ibx, jm;
	static nec_float rbk, amg, den, denm;
	nec_complex a1, a2, as1, as2, del, aa;
	nec_complex q1[6][20], q2[6][20];
	 
	complex_array ans1(6), ans2(6);
	
	rbk=real(bk);
	del=dela;
	if (ibk == 0)
		ibx=1;
	else
		ibx=0;
	
	// I believe that this is a spurious generalization for this routine. Hence the
	// following assert.
	ASSERT(nans == 6);
	
	for (int i = 0; i < nans; i++ )
		ans2[i]=seed[i];
	
	m_contour_b=start;
	for (int intx = 1; intx <= MAXH; intx++ )
	{
		int inx=intx-1;
		m_contour_a = m_contour_b;
		m_contour_b += del;
		
		if ( (ibx == 0) && (real(m_contour_b) >= rbk) )
		{
			/* hit break point.  reset seed and start over. */
			ibx=1;
			m_contour_b=bk;
			del=delb;
			rom1(nans,sum,2);
			if ( ibx != 2 )
			{
				for (int i = 0; i < nans; i++ )
					ans2[i] += sum[i];
				intx = 0;
				continue;
			}
			
			for (int i = 0; i < nans; i++ )
				ans2[i]=ans1[i]+sum[i];
			intx = 0;
			continue;
		} /* if ( (ibx == 0) && (real(m_contour_b) >= rbk) ) */
		
		rom1(nans,sum,2);
		for (int i = 0; i < nans; i++ )
			ans1[i] = ans2[i]+sum[i];
			
		m_contour_a = m_contour_b;
		m_contour_b += del;
		
		if ( (ibx == 0) && (real(m_contour_b) >= rbk) )
		{
			/* hit break point.  reset seed and start over. */
			ibx=2;
			m_contour_b=bk;
			del=delb;
			rom1(nans,sum,2);
			if ( ibx != 2 )
			{
				for (int i = 0; i < nans; i++ )
					ans2[i] += sum[i];
				intx = 0;
				continue;
			}
			
			for (int i = 0; i < nans; i++ )
				ans2[i] = ans1[i]+sum[i];
			intx = 0;
			continue;
		} /* if ( (ibx == 0) && (real(m_contour_b) >= rbk) ) */
		
		rom1(nans,sum,2);
		for (int i = 0; i < nans; i++ )
			ans2[i]=ans1[i]+sum[i];
		
		den=0.;
		for (int i = 0; i < nans; i++ )
		{
			as1=ans1[i];
			as2=ans2[i];
			
			if (intx >= 2)
			{
				for (int j = 1; j < intx; j++ )
				{
					jm=j-1;
					aa=q2[i][jm];
					a1=q1[i][jm]+as1-2.*aa;
				
					if ( (real(a1) != 0.) || (imag(a1) != 0.) )
					{
						a2=aa-q1[i][jm];
						a1=q1[i][jm]-a2*a2/a1;
					}
					else
						a1=q1[i][jm];
				
					a2=aa+as2-2.*as1;
					if ( (real(a2) != 0.) || (imag(a2) != 0.) )
						a2=aa-(as1-aa)*(as1-aa)/a2;
					else
						a2=aa;
				
					q1[i][jm]=as1;
					q2[i][jm]=as2;
					as1=a1;
					as2=a2;	
				}
			}
			
			q1[i][intx-1]=as1;
			q2[i][intx-1]=as2;
			amg=fabs(real(as2))+fabs(imag(as2));
			if (amg > den)
				den=amg;		
		} /* for ( i = 0; i < nans; i++ ) */
		
		denm=1.e-3*den*CRIT;
		jm=intx-3;
		if (jm < 1)
			jm = 1;
		
		for (int j = jm-1; j < intx; j++ )
		{
			brk = false;
			for (int i = 0; i < nans; i++ )
			{
				a1=q2[i][j];
				den=(fabs(real(a1))+fabs(imag(a1)))*CRIT;
				if (den < denm)
					den=denm;
				a1=q1[i][j]-a1;
				amg=fabs(real(a1)+fabs(imag(a1)));
				if (amg > den)
				{
					brk = true;
					break;
				}
			
			} /* for ( i = 0; i < nans; i++ ) */
			
			if ( brk )
				break;	
		} /* for ( j = jm-1; j < intx; j++ ) */
		
		if ( false == brk )
		{
			for (int i = 0; i < nans; i++ )
				sum[i]=.5*(q1[i][inx]+q2[i][inx]);
			return;
		}
		
	} /* for ( intx = 1; intx <= maxh; intx++ ) */
	
	/* No convergence */
	throw new nec_exception("No convergence in gshank() - aborting");
}


/*! \brief rom1 integrates the 6 Sommerfeld integrals from m_contour_a to m_contour_b in lambda.
	The method of variable interval width Romberg integration is used. */
void c_evlcom::rom1( int n, complex_array& sum, int nx )
{
	int ns, nt;
	static nec_float z, ze, s, ep, zend, dz=0.0, dzot=0.0, tr, ti;
	static nec_complex t00, t11, t02;
	
	static complex_array g1(6), g2(6), g3(6), g4(6), g5(6), t01(6), t10(6), t20(6);
	
	ASSERT(n == 6);
	
	z = 0.0;
	ze = 1.0;
	s = 1.0;
	ep = s / (1.0e4 * NM);
	zend=ze-ep;
	
	nec_complex cmplx_zero(0.0,0.0);
	for (int i = 0; i < n; i++ )
		sum[i]=cmplx_zero;
		
	ns=nx;
	nt=0;
	saoa(z,g1);
	
	bool jump = false;
	bool lstep = false;
	
	while( true )
	{
		if ( false == jump )
		{
			dz = s/ns;
			if ( (z+dz) > ze )
			{
				dz=ze-z;
				if ( dz <= ep )
					return;
			}
			
			dzot=dz*.5;
			saoa(z+dzot,g3);
			saoa(z+dz,g5);	
		} /* if ( false == jump ) */
		
		bool nogo = false;
		for (int i = 0; i < n; i++ )
		{
			t00=(g1[i]+g5[i])*dzot;
			t01[i]=(t00+dz*g3[i])*.5;
			t10[i]=(4.*t01[i]-t00)/3.;
			
			/* test convergence of 3 point romberg result */
			test( real(t01[i]), real(t10[i]), &tr, imag(t01[i]), imag(t10[i]), &ti, 0. );
			if ( (tr > CRIT) || (ti > CRIT) )
				nogo = true;
		}
		
		if ( false == nogo )
		{
			for (int i = 0; i < n; i++ )
				sum[i] += t10[i];
			
			nt += 2;
			z += dz;
			if (z > zend)
				return;
			
			for (int i = 0; i < n; i++ )
				g1[i]=g5[i];
			
			if ( (nt >= NTS) && (ns > nx) )
			{
				ns=ns/2;
				nt=1;
			}
			
			jump = false;
			continue;
		} /* if ( false == nogo ) */
		
		saoa(z+dz*.25,g2);
		saoa(z+dz*.75,g4);
		nogo = false;
		for (int i = 0; i < n; i++ )
		{
			t02=(t01[i]+dzot*(g2[i]+g4[i]))*.5;
			t11=(4.*t02-t01[i])/3.;
			t20[i]=(16.*t11-t10[i])/15.;
			
			/* test convergence of 5 point Romberg result */
			test( real(t11), real(t20[i]), &tr, imag(t11), imag(t20[i]), &ti, 0.0 );
			if ( (tr > CRIT) || (ti > CRIT) )
				nogo = true;
		}
		
		if ( false == nogo )
		{
			for (int i = 0; i < n; i++ )
				sum[i] += t20[i];
			
			nt++;
			z += dz;
			if (z > zend)
				return;
			
			for (int i = 0; i < n; i++ )
				g1[i]=g5[i];
			
			if ( (nt >= NTS) && (ns > nx) )
			{
				ns=ns/2;
				nt=1;
			}
			
			jump = false;
			continue;	
		} /* if ( false == nogo ) */
		
		nt=0;
		if (ns < NM)
		{
			ns *= 2;
			dz=s/ns;
			dzot=dz*.5;
			
			for (int i = 0; i < n; i++ )
			{
				g5[i]=g3[i];
				g3[i]=g2[i];
			}
			
			jump = true;
			continue;	
		} /* if (ns < NM) */
		
		if ( false == lstep )
		{
			lstep = true;
			lambda( z, &t00, &t11 );
		}
		
		for (int i = 0; i < n; i++ )
			sum[i] += t20[i];
		
		nt++;
		z += dz;
		if (z > zend)
			return;
		
		for (int i = 0; i < n; i++ )
			g1[i]=g5[i];
		
		if ( (nt >= NTS) && (ns > nx) )
		{
			ns /= 2;
			nt=1;
		}
		
		jump = false;
	} /* while( TRUE ) */
}

/*! \brief saoa computes the integrand for each of the 6 Sommerfeld
	integrals for source and observer above ground. */
void c_evlcom::saoa( nec_float t, complex_array& ans)
{
	static nec_complex xl, dxl, cgam1, cgam2, b0, b0p, com, dgam, den1, den2;
	
	lambda(t, &xl, &dxl);
	if ( m_bessel_flag == true )
	{
		/* Bessel function form */
		bessel(xl*m_rho, &b0, &b0p);
		b0  *=2.;
		b0p *=2.;
		cgam1=sqrt(xl*xl-m_ck1sq);
		cgam2=sqrt(xl*xl-m_ck2sq);
		if (real(cgam1) == 0.0)
			cgam1=nec_complex(0.0,-fabs(imag(cgam1)));
		if (real(cgam2) == 0.)
			cgam2=nec_complex(0.0,-fabs(imag(cgam2)));
	}
	else
	{
		/* Hankel function form */
		hankel(xl*m_rho, &b0, &b0p);
		com=xl-m_ck1;
		cgam1=sqrt(xl+m_ck1)*sqrt(com);
		if (real(com) < 0. && imag(com) >= 0.)
			cgam1=-cgam1;
		com=xl-m_ck2;
		cgam2=sqrt(xl+m_ck2)*sqrt(com);
		if (real(com) < 0. && imag(com) >= 0.)
			cgam2=-cgam2;
	}
	
	if (norm(xl) >= m_tsmag)
	{
		if (imag(xl) >= 0.0)
		{
			nec_float xlr = real(xl);
			if (xlr >= m_ck2)
			{
				if (xlr <= m_ck1r)
					dgam=cgam2-cgam1;
				else
				{
					nec_float sign =1.0;
					dgam=1.0/(xl*xl);
					dgam=sign*((m_ct3*dgam+m_ct2)*dgam+m_ct1)/xl;