nec_context.cpp

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

				stop(-1);
			}
			... rest of code here
		}
		// done with loads.
		if ( iwarn == true )
			m_output.nec_printf(
				"\n  NOTE, SOME OF THE ABOVE SEGMENTS "
				"HAVE BEEN LOADED TWICE - IMPEDANCES ADDED" );
	
		if ( nop != 1)
		{
			for( i = 0; i < np; i++ )
			{
				zt= zarray[i];
				l1= i;
			
				for( l2 = 1; l2 < nop; l2++ )
				{
					l1 += np;
					zarray[l1]= zt;
				}
			}
		}
		
		
	*/
	/* cycle over loading cards */
	while( true )
	{
		istepx = istep;
		istep++;
	
		if ( istep > nload)
		{
			if ( iwarn == true )
				m_output.line("NOTE, SOME OF THE ABOVE SEGMENTS HAVE BEEN LOADED TWICE - IMPEDANCES ADDED" );
		
			if ( nop == 1)
				return;
		
			for(int i = 0; i < np; i++ )
			{
				zt= zarray[i];
				l1= i;
			
				for( l2 = 1; l2 < nop; l2++ )
				{
					l1 += np;
					zarray[l1]= zt;
				}
			}
			return;
		
		} /* if ( istep > nload) */
	
		if ( ldtyp[istepx] > 5 )
		{
			nec_stop("IMPROPER LOAD TYPE CHOSEN,"
				" REQUESTED TYPE IS %d", ldtyp[istepx] );
		}
	
		/* search segments for proper itags */
		ldtags = ldtag[istepx];
		int jump = ldtyp[istepx]+1;
		ichk=0;
		l1= 1;
		l2= n;
	
		if ( ldtags == 0)
		{
			if ( (ldtagf[istepx] != 0) || (ldtagt[istepx] != 0) )
			{
				l1= ldtagf[istepx];
				l2= ldtagt[istepx];
			}
		}
	
		for(int i = l1-1; i < l2; i++ )
		{
			if ( ldtags != 0)
			{
				if ( ldtags != m_geometry->segment_tags[i])
					continue;
			
				if ( ldtagf[istepx] != 0)
				{
					ichk++;
					if ( (ichk < ldtagf[istepx]) || (ichk > ldtagt[istepx]) )
						continue;
				}
				else
					ichk=1;
			} /* if ( ldtags != 0) */
			else
				ichk=1;
		
			ASSERT(0.0 != m_geometry->segment_length[i]);
			ASSERT(0.0 != m_geometry->segment_radius[i]);
			
			/* 	Calculation of lamda*imped. per unit length,
				jump to appropriate section for loading type */
			switch( jump )
			{
				case 1:
					zt= zlr[istepx]/ m_geometry->segment_length[i]+ tpcj* zli[istepx]/( m_geometry->segment_length[i]*wavelength);
					if ( fabs( zlc[istepx]) > 1.0e-20)
						zt += wavelength/( tpcj* m_geometry->segment_length[i]* zlc[istepx]);
					break;
			
				case 2:
					zt= tpcj* m_geometry->segment_length[i]* zlc[istepx]/ wavelength;
					if ( fabs( zli[istepx]) > 1.0e-20)
						zt += m_geometry->segment_length[i]* wavelength/( tpcj* zli[istepx]);
					if ( fabs( zlr[istepx]) > 1.0e-20)
						zt += m_geometry->segment_length[i]/ zlr[istepx];
					zt=1./ zt;
					break;
			
				case 3:
					zt= zlr[istepx]* wavelength+ tpcj* zli[istepx];
					if ( fabs( zlc[istepx]) > 1.0e-20)
						zt += 1./( tpcj* m_geometry->segment_length[i]* m_geometry->segment_length[i]* zlc[istepx]);
					break;
			
				case 4:
					zt= tpcj* m_geometry->segment_length[i]* m_geometry->segment_length[i]* zlc[istepx];
					if ( fabs( zli[istepx]) > 1.0e-20)
						zt += 1./( tpcj* zli[istepx]);
					if ( fabs( zlr[istepx]) > 1.0e-20)
						zt += 1./( zlr[istepx]* wavelength);
					zt=1./ zt;
					break;
			
				case 5:
					zt= nec_complex( zlr[istepx], zli[istepx])/ m_geometry->segment_length[i];
					break;
			
				case 6:
				{
					zt= zint( zlr[istepx]* wavelength, m_geometry->segment_radius[i]);
				}
			
			} /* switch( jump ) */
		
			if (( fabs( real( zarray[i]))+ fabs( imag( zarray[i]))) > 1.0e-20)
				iwarn=true;
			zarray[i] += zt;
		
		} /* for( i = l1-1; i < l2; i++ ) */
		
		if ( ichk == 0 )
		{
			nec_exception* nex = new nec_exception("LOADING DATA CARD ERROR, NO SEGMENT HAS AN ITAG = ");
			nex->append(ldtags);
			throw nex;
		}
	
		/* printing the segment loading data, jump to proper print */
		switch( jump )
		{
			case 1:
				impedance_print( ldtags, ldtagf[istepx], ldtagt[istepx], zlr[istepx],
					zli[istepx], zlc[istepx],0.,0.,0.," SERIES ");
				break;
		
			case 2:
				impedance_print( ldtags, ldtagf[istepx], ldtagt[istepx], zlr[istepx],
					zli[istepx], zlc[istepx],0.,0.,0.,"PARALLEL");
				break;
		
			case 3:
				impedance_print( ldtags, ldtagf[istepx], ldtagt[istepx], zlr[istepx],
					zli[istepx], zlc[istepx],0.,0.,0., "SERIES (PER METER)");
				break;
		
			case 4:
				impedance_print( ldtags, ldtagf[istepx], ldtagt[istepx], zlr[istepx],
					zli[istepx], zlc[istepx],0.,0.,0.,"PARALLEL (PER METER)");
				break;
		
			case 5:
				impedance_print( ldtags, ldtagf[istepx], ldtagt[istepx],0.,0.,0.,
					zlr[istepx], zli[istepx],0.,"FIXED IMPEDANCE ");
				break;
		
			case 6:
				impedance_print( ldtags, ldtagf[istepx], ldtagt[istepx],
					0.,0.,0.,0.,0., zlr[istepx],"  WIRE  ");
		} /* switch( jump ) */
	
	} /* while( true ) */
}

/* cmset sets up the complex structure matrix in the array in_cm */
void nec_context::cmset( int nrow, complex_array& in_cm, nec_float rkhx)
{
	int mp2, iout, it, j, i1, i2, in2;
	int im1, im2, ist, ij, jss, jm1, jm2, jst, k, kk;
	complex_array scm;
	
	int np = m_geometry->np;
	int mp = m_geometry->mp;
	
	mp2 = 2 * mp;
//	 = np + mp2;
//	neq= m_geometry->n_plus_2m; // n + 2* m;
	
	rkh= rkhx;
	iout=2* npblk* nrow;
	it= nlast;
			
	in_cm.fill(0,it*nrow,cplx_00());
	
	i1= 1;
	i2= it;
	in2= i2;
	
	if ( in2 > np)
		in2= np;
	
	im1= i1- np;
	im2= i2- np;
	
	if ( im1 < 1)
		im1=1;
	
	ist=1;
	if ( i1 <= np)
		ist= np- i1+2;
	
	/* wire source loop */
	int n = m_geometry->n;
		
	for( j = 1; j <= n; j++ )
	{
		m_geometry->trio(j);
		
		for (int i = 0; i < m_geometry->jsno; i++ )
		{
			ij = m_geometry->jco[i];
			m_geometry->jco[i] = ((ij-1)/ np)* mp2+ij;
		}
	
		if ( i1 <= in2)
			cmww( j, i1, in2, in_cm, nrow, in_cm, nrow,1);
	
		if ( im1 <= im2)
		{
			complex_array temp = in_cm.sub_array((ist-1)*nrow);
			cmws( j, im1, im2, temp, nrow, in_cm, 1);
		}
		/* matrix elements modified by loading */
		if ( nload == 0)
			continue;
	
		if ( j > np)
			continue;
	
		int ipr = j;
		if ( (ipr < 1) || (ipr > it) )
			continue;
	
		nec_complex zaj= zarray[j-1];
	
		for (int i = 0; i < m_geometry->jsno; i++ )
		{
			jss = m_geometry->jco[i];
			in_cm[(jss-1)+(ipr-1)*nrow] -= ( m_geometry->ax[i]+ m_geometry->cx[i])* zaj;
		}
	} /* for( j = 1; j <= n; j++ ) */
	
	int m = m_geometry->m;
	if ( m != 0)
	{
		/* matrix elements for patch current sources */
		jm1=1- mp;
		jm2=0;
		jst=1- mp2;
	
		for (int i = 0; i < nop; i++ )
		{
			jm1 += mp;
			jm2 += mp;
			jst += npeq;
		
			if ( i1 <= in2)
			{
				complex_array temp = in_cm.sub_array((jst-1));
				cmsw( jm1, jm2, i1, in2, temp, in_cm, 0, nrow, 1);
			}
		
			if ( im1 <= im2)
			{
				complex_array temp = in_cm.sub_array((jst-1)+(ist-1)*nrow);
				compute_matrix_ss( jm1, jm2, im1, im2, temp, nrow, 1);
			}
		}
	
	} /* if ( m != 0) */
	
	if ( icase == 1)
		return;
	
	/* Allocate to scratch memory */
	scm.resize(m_geometry->n_plus_2m);
	
	/* combine elements for symmetry modes */
	for (int i = 0; i < it; i++ )
	{
		int row_offset = i*nrow;
		
		for( j = 0; j < npeq; j++ )
		{
			for( k = 0; k < nop; k++ )
			{
				int ka = j+ k*npeq;
				scm[k] = in_cm[row_offset + ka];
			}
		
			in_cm[row_offset + j] = scm.sum(0,nop);
			
			for( k = 1; k < nop; k++ )
			{
				int ka = j+ k*npeq;
				nec_complex deter = scm[0];
			
				for( kk = 1; kk < nop; kk++ )
				{
					deter += scm[kk] * symmetry_array[k+kk*nop];
				}
				in_cm[row_offset + ka] = deter;
			}
		}
	}
	
	scm.resize(0);
}

/* compute_matrix_ss computes matrix elements for surface-surface interactions. */

void nec_context::compute_matrix_ss( int j1, int j2, int im1, int im2,
    complex_array& in_cm, int nrow, int itrp )
{
	int i1, i2, icomp, ii1, i, il, ii2, jj1, jl, jj2;
	nec_float t1xi, t1yi, t1zi, t2xi, t2yi, t2zi, xi, yi, zi;
	nec_complex g11, g12, g21, g22;
	
	i1=( im1+1)/2;
	i2=( im2+1)/2;
	icomp= i1*2-3;
	ii1=-2;
	if ( icomp+2 < im1)
		ii1=-3;
	
	/* loop over observation patches */
	il = -1;
	for( i = i1; i <= i2; i++ )
	{
		il++;
		icomp += 2;
		ii1 += 2;
		ii2 = ii1+1;
	
		t1xi= m_geometry->t1x[il]* m_geometry->psalp[il];
		t1yi= m_geometry->t1y[il]* m_geometry->psalp[il];
		t1zi= m_geometry->t1z[il]* m_geometry->psalp[il];
		t2xi= m_geometry->t2x[il]* m_geometry->psalp[il];
		t2yi= m_geometry->t2y[il]* m_geometry->psalp[il];
		t2zi= m_geometry->t2z[il]* m_geometry->psalp[il];
		xi= m_geometry->px[il];
		yi= m_geometry->py[il];
		zi= m_geometry->pz[il];
	
		/* loop over source patches */
		jj1=-2;
		for(int j = j1; j <= j2; j++ )
		{
			jl=j-1;
			jj1 += 2;
			jj2 = jj1+1;
		
			m_s= m_geometry->pbi[jl];
			xj= m_geometry->px[jl];
			yj= m_geometry->py[jl];
			zj= m_geometry->pz[jl];
			t1xj= m_geometry->t1x[jl];
			t1yj= m_geometry->t1y[jl];
			t1zj= m_geometry->t1z[jl];
			t2xj= m_geometry->t2x[jl];
			t2yj= m_geometry->t2y[jl];
			t2zj= m_geometry->t2z[jl];
		
			hintg( xi, yi, zi);
		
			g11=-( t2xi* exk+ t2yi* eyk+ t2zi* ezk);
			g12=-( t2xi* exs+ t2yi* eys+ t2zi* ezs);
			g21=-( t1xi* exk+ t1yi* eyk+ t1zi* ezk);
			g22=-( t1xi* exs+ t1yi* eys+ t1zi* ezs);
		
			if ( i == j )
			{
				g11 -= .5;
				g22 += .5;
			}
		
			/* normal fill */
			if ( itrp == 0)
			{
				if ( icomp >= im1 )
				{
					in_cm[ii1+jj1*nrow]= g11;
					in_cm[ii1+jj2*nrow]= g12;
				}
			
				if ( icomp >= im2 )
					continue;
			
				in_cm[ii2+jj1*nrow]= g21;
				in_cm[ii2+jj2*nrow]= g22;
				continue;	
			} /* if ( itrp == 0) */
		
			/* transposed fill */
			if ( icomp >= im1 )
			{
				in_cm[jj1+ii1*nrow]= g11;
				in_cm[jj2+ii1*nrow]= g12;
			}
		
			if ( icomp >= im2 )
				continue;
		
			in_cm[jj1+ii2*nrow]= g21;
			in_cm[jj2+ii2*nrow]= g22;
		
		} /* for( j = j1; j <= j2; j++ ) */
	} /* for( i = i1; i <= i2; i++ ) */
}
	

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

/* computes matrix elements for e along wires due to patch current */
void nec_context::cmsw( int j1, int j2, int i1, int i2, complex_array& in_cm,
    complex_array& cw, int ncw, int nrow, int itrp )
{
  int neqs, k, icgo, i, ipch, jl, j, js, il;
  int jsnox; /* -1 offset to "m_geometry->jsno" for array indexing */
  nec_float xi, yi, zi, cabi, sabi, salpi, fsign=1., pyl, pxl;
  complex_array emel;
  
  emel.resize(9);

  neqs= m_geometry->n_plus_2m;
  jsnox = m_geometry->jsno-1;

  if ( itrp >= 0)
  {
    k=-1;
    icgo=0;

    /* observation loop */
    for( i = i1-1; i < i2; i++ )
    {
      k++;
      xi= m_geometry->x[i];
      yi= m_geometry->y[i];
      zi= m_geometry->z[i];
      cabi= m_geometry->cab[i];
      sabi= m_geometry->sab[i];
      salpi= m_geometry->salp[i];
      ipch=0;

      if ( m_geometry->icon1[i] >= PCHCON)
      {
	ipch= m_geometry->icon1[i]-PCHCON;
	fsign=-1.;
      }

      if ( m_geometry->icon2[i] >= PCHCON)
      {
	ipch= m_geometry->icon2[i]-PCHCON;
	fsign=1.;
      }

      /* source loop */
      jl = -1;
      for( j = j1; j <= j2; j++ )
      {
	jl += 2;
	js = j-1;
	t1xj= m_geometry->t1x[js];
	t1yj= m_geometry->t1y[js];
	t1zj= m_geometry->t1z[js];
	t2xj= m_geometry->t2x[js];
	t2yj= m_geometry->t2y[js];
	t2zj= m_geometry->t2z[js];
	xj= m_geometry->px[js];
	yj= m_geometry->py[js];
	zj= m_geometry->pz[js];
	m_s = m_geometry->pbi[js];

	/* ground loop */
	#warning weird problem to test here. replaced loop variable called ip
	for (int ipgnd = 1; ipgnd <= ground.ksymp; ipgnd++ )
	{

	  if ( ((ipch == j) || (icgo != 0)) && (ipgnd != 2) )
	  {
	    if ( icgo <= 0 )
	    {
	      ASSERT(ipgnd !