c_geometry.cpp

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

	{
		m_output->nec_printf( "\n\n\n"
			"                                   "
			" --------- SURFACE PATCH DATA ---------\n"
			"                                            "
			" COORDINATES IN METERS\n\n"
			" PATCH      COORD. OF PATCH CENTER           UNIT NORMAL VECTOR      "
			" PATCH           COMPONENTS OF UNIT TANGENT VECTORS\n"
			"  NO:       X          Y          Z          X        Y        Z      "
			" AREA         X1       Y1       Z1        X2       Y2      Z2" );
	
		for(int i = 0; i < m; i++ )
		{
			nec_float xw1=( t1y[i]* t2z[i]- t1z[i]* t2y[i])* psalp[i];
			nec_float yw1=( t1z[i]* t2x[i]- t1x[i]* t2z[i])* psalp[i];
			nec_float zw1=( t1x[i]* t2y[i]- t1y[i]* t2x[i])* psalp[i];
	
			m_output->nec_printf( "\n"
				" %4d %10.5f %10.5f %10.5f  %8.4f %8.4f %8.4f"
				" %10.5f  %8.4f %8.4f %8.4f  %8.4f %8.4f %8.4f",
				i+1, px[i], py[i], pz[i], xw1, yw1, zw1, pbi[i],
				t1x[i], t1y[i], t1z[i], t2x[i], t2y[i], t2z[i] );
		} /* for( i = 0; i < m; i++ ) */
	} /* if ( m == 0) */

	n_plus_m  = n+m;
	n_plus_2m = n+2*m;
	n_plus_3m = n+3*m;

	x_unscaled.resize(n);
	y_unscaled.resize(n);
	z_unscaled.resize(n);
	si_unscaled.resize(n);
	bi_unscaled.resize(n);
	
	px_unscaled.resize(m);
	py_unscaled.resize(m);
	pz_unscaled.resize(m);
	pbi_unscaled.resize(m);
	
	// Fill the unscaled segments...
	for (int i = 0; i < n; i++ )
	{
		x_unscaled[i]= x[i];
		y_unscaled[i]= y[i];
		z_unscaled[i]= z[i];
		si_unscaled[i]= segment_length[i];
		bi_unscaled[i]= segment_radius[i];
	}
	// Fill the unscaled patches...
	for (int i = 0; i < m; i++ )
	{
		px_unscaled[i]= px[i];
		py_unscaled[i]= py[i];
		pz_unscaled[i]= pz[i];
		pbi_unscaled[i]= pbi[i];
	}
}


/* subroutine wire generates segment geometry */
/* data for a straight wire of segment_count segments. */
void c_geometry::wire( int tag_id, int segment_count, nec_float xw1, nec_float yw1, nec_float zw1,
    nec_float xw2, nec_float yw2, nec_float zw2, nec_float rad,
    nec_float rdel, nec_float rrad )
{
	nec_float xd, yd, zd, delz, rd, fns, radz;
	nec_float xs1, ys1, zs1, xs2, ys2, zs2;
	
	int istart = n;
	n = n + segment_count;
	np= n;
	mp= m;
	m_ipsym=0;
	
	if ( segment_count < 1)
		return;
	
	/* Reallocate tags buffer */
	segment_tags.resize(n + m);
	
	/* Reallocate wire buffers */
	x.resize(n);
	y.resize(n);
	z.resize(n);
	x2.resize(n);
	y2.resize(n);
	z2.resize(n);
	segment_radius.resize(n);
	
	xd= xw2- xw1;
	yd= yw2- yw1;
	zd= zw2- zw1;
	
	if ( fabs(rdel-1.0) >= 1.0e-6)
	{
		delz= sqrt( xd* xd+ yd* yd+ zd* zd);
		xd= xd/ delz;
		yd= yd/ delz;
		zd= zd/ delz;
		delz= delz*(1.- rdel)/(1.- pow(rdel, segment_count) );
		rd= rdel;
	}
	else
	{
		fns= segment_count;
		xd= xd/ fns;
		yd= yd/ fns;
		zd= zd/ fns;
		delz=1.0;
		rd=1.0;
	}
	
	radz= rad;
	xs1= xw1;
	ys1= yw1;
	zs1= zw1;
	
	for (int i = istart; i < n; i++ )
	{
		segment_tags[i]= tag_id;
		xs2= xs1+ xd* delz;
		ys2= ys1+ yd* delz;
		zs2= zs1+ zd* delz;
		x[i]= xs1;
		y[i]= ys1;
		z[i]= zs1;
		x2[i]= xs2;
		y2[i]= ys2;
		z2[i]= zs2;
		ASSERT(0.0 != radz);
		segment_radius[i]= radz;
		delz= delz* rd;
		radz= radz* rrad;
		xs1= xs2;
		ys1= ys2;
		zs1= zs2;
	}
	
	x2[n-1]= xw2;
	y2[n-1]= yw2;
	z2[n-1]= zw2;
}

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

/* subroutine helix generates segment geometry */
/* data for a helix of segment_count segments */
void c_geometry::helix( nec_float s, nec_float hl, nec_float a1, nec_float b1,
    nec_float a2, nec_float b2, nec_float rad, int segment_count, int tag_id )
{
	int ist;
	nec_float turns, zinc, copy, sangle, hdia, turn, pitch, hmaj, hmin;
	
	ist= n;
	n += segment_count;
	np= n;
	mp= m;
	m_ipsym=0;
	
	if ( segment_count < 1)
		return;
	
	turns= fabs( hl/ s);
	zinc= fabs( hl/ segment_count);
	
	segment_tags.resize(n+m); /*????*/
	
	/* Reallocate wire buffers */
	x.resize(n);
	y.resize(n);
	z.resize(n);
	x2.resize(n);
	y2.resize(n);
	z2.resize(n);
	segment_radius.resize(n);
	
	z[ist]=0.;
	for(int i = ist; i < n; i++ )
	{
		segment_radius[i]= rad;
		segment_tags[i]= tag_id;
	
		if ( i != ist )
		z[i]= z[i-1]+ zinc;
	
		z2[i]= z[i]+ zinc;
	
		if ( a2 == a1)
		{
			if ( b1 == 0.)
				b1= a1;
		
			x[i]= a1* cos(2.* pi()* z[i]/ s);
			y[i]= b1* sin(2.* pi()* z[i]/ s);
			x2[i]= a1* cos(2.* pi()* z2[i]/ s);
			y2[i]= b1* sin(2.* pi()* z2[i]/ s);
		}
		else
		{
			if ( b2 == 0.)
				b2= a2;
		
			x[i]=( a1+( a2- a1)* z[i]/ fabs( hl))* cos(2.* pi()* z[i]/ s);
			y[i]=( b1+( b2- b1)* z[i]/ fabs( hl))* sin(2.* pi()* z[i]/ s);
			x2[i]=( a1+( a2- a1)* z2[i]/ fabs( hl))* cos(2.* pi()* z2[i]/ s);
			y2[i]=( b1+( b2- b1)* z2[i]/ fabs( hl))* sin(2.* pi()* z2[i]/ s);	
		} /* if ( a2 == a1) */
	
		if ( hl > 0.)
			continue;
	
		copy= x[i];
		x[i]= y[i];
		y[i]= copy;
		copy= x2[i];
		x2[i]= y2[i];
		y2[i]= copy;	
	} /* for( i = ist; i < n; i++ ) */
	
	if ( a2 != a1)
	{
		sangle= atan( a2/( fabs( hl)+( fabs( hl)* a1)/( a2- a1)));
		m_output->nec_printf(
			"\n       THE CONE ANGLE OF THE SPIRAL IS %10.4f", sangle );
		return;
	}
	
	if ( a1 == b1)
	{
		hdia=2.* a1;
		turn= hdia* pi();
		pitch= atan( s/( pi()* hdia));
		turn= turn/ cos( pitch);
		pitch=180.* pitch/ pi();
	}
	else
	{
		if ( a1 >= b1)
		{
			hmaj=2.* a1;
			hmin=2.* b1;
		}
		else
		{
			hmaj=2.* b1;
			hmin=2.* a1;
		}
	
		hdia= sqrt(( hmaj*hmaj+ hmin*hmin)/2* hmaj);
		turn=2.* pi()* hdia;
		pitch=(180./ pi())* atan( s/( pi()* hdia));
	
	} /* if ( a1 == b1) */
	
	m_output->nec_printf( "\n"
		"       THE PITCH ANGLE IS: %.4f    THE LENGTH OF WIRE/TURN IS: %.4f",
		pitch, turn );
}

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


/* subroutine move moves the structure with respect to its */
/* coordinate system or reproduces structure in new positions. */
/* structure is rotated about x,y,z axes by rox,roy,roz */
/* respectively, then shifted by xs,ys,zs */
void c_geometry::move( nec_float rox, nec_float roy, nec_float roz, nec_float xs,
    nec_float ys, nec_float zs, int its, int nrpt, int itgi )
{
  int nrp, ix, i1, k, ir, i, ii;
  nec_float sps, cps, sth, cth, sph, cph, xx, xy;
  nec_float xz, yx, yy, yz, zx, zy, zz, xi, yi, zi;

  if ( fabs( rox)+ fabs( roy) > 1.0e-10)
    m_ipsym= m_ipsym*3;

  sps= sin( rox);
  cps= cos( rox);
  sth= sin( roy);
  cth= cos( roy);
  sph= sin( roz);
  cph= cos( roz);
  xx= cph* cth;
  xy= cph* sth* sps- sph* cps;
  xz= cph* sth* cps+ sph* sps;
  yx= sph* cth;
  yy= sph* sth* sps+ cph* cps;
  yz= sph* sth* cps- cph* sps;
  zx=- sth;
  zy= cth* sps;
  zz= cth* cps;

  if ( nrpt == 0)
    nrp=1;
  else
    nrp= nrpt;

  ix=1;
  if ( n > 0)
  {
    i1= get_segment_number( its, 1);
    if ( i1 < 1)
      i1= 1;

    ix= i1;
    if ( nrpt == 0)
      k= i1-1;
    else
    {
      k= n;
      /* Reallocate tags buffer */
      segment_tags.resize(n+m + (n+1-i1)*nrpt);
      // mreq = n+m + (n+1-i1)*nrpt;
      // segment_tags.resize(mreq);

      /* Reallocate wire buffers */
      int new_size = (n+(n+1-i1)*nrpt);
      x.resize(new_size);
      y.resize(new_size);
      z.resize(new_size);
      x2.resize(new_size);
      y2.resize(new_size);
      z2.resize(new_size);
      segment_radius.resize(new_size);
    }

    for( ir = 0; ir < nrp; ir++ )
    {
      for( i = i1-1; i < n; i++ )
      {
	xi= x[i];
	yi= y[i];
	zi= z[i];
	x[k]= xi* xx+ yi* xy+ zi* xz+ xs;
	y[k]= xi* yx+ yi* yy+ zi* yz+ ys;
	z[k]= xi* zx+ yi* zy+ zi* zz+ zs;
	xi= x2[i];
	yi= y2[i];
	zi= z2[i];
	x2[k]= xi* xx+ yi* xy+ zi* xz+ xs;
	y2[k]= xi* yx+ yi* yy+ zi* yz+ ys;
	z2[k]= xi* zx+ yi* zy+ zi* zz+ zs;
	segment_radius[k]= segment_radius[i];
	segment_tags[k]= segment_tags[i];
	if ( segment_tags[i] != 0)
	  segment_tags[k]= segment_tags[i]+ itgi;

	k++;

      } /* for( i = i1; i < n; i++ ) */

      i1= n+1;
      n= k;

    } /* for( ir = 0; ir < nrp; ir++ ) */

  } /* if ( n >= n2) */

  if ( m > 0)
  {
    i1 = 0;
    if ( nrpt == 0)
      k= 0;
    else
      k = m;

    /* Reallocate patch buffers */
    int new_size = m * (1+nrpt);
    px.resize(new_size);
    py.resize(new_size);
    pz.resize(new_size);
    t1x.resize(new_size);
    t1y.resize(new_size);
    t1z.resize(new_size);
    t2x.resize(new_size);
    t2y.resize(new_size);
    t2z.resize(new_size);
    pbi.resize(new_size);
    psalp.resize(new_size);

    for( ii = 0; ii < nrp; ii++ )
    {
      for( i = i1; i < m; i++ )
      {
	xi= px[i];
	yi= py[i];
	zi= pz[i];
	px[k]= xi* xx+ yi* xy+ zi* xz+ xs;
	py[k]= xi* yx+ yi* yy+ zi* yz+ ys;
	pz[k]= xi* zx+ yi* zy+ zi* zz+ zs;
	xi= t1x[i];
	yi= t1y[i];
	zi= t1z[i];
	t1x[k]= xi* xx+ yi* xy+ zi* xz;
	t1y[k]= xi* yx+ yi* yy+ zi* yz;
	t1z[k]= xi* zx+ yi* zy+ zi* zz;
	xi= t2x[i];
	yi= t2y[i];
	zi= t2z[i];
	t2x[k]= xi* xx+ yi* xy+ zi* xz;
	t2y[k]= xi* yx+ yi* yy+ zi* yz;
	t2z[k]= xi* zx+ yi* zy+ zi* zz;
	psalp[k]= psalp[i];
	pbi[k]= pbi[i];
	k++;

      } /* for( i = i1; i < m; i++ ) */

      i1= m;
      m = k;

    } /* for( ii = 0; ii < nrp; ii++ ) */

  } /* if ( m >= m2) */

  if ( (nrpt == 0) && (ix == 1) )
    return;

  np= n;
  mp= m;
  m_ipsym=0;

  return;
}

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

/* reflects partial structure along x,y, or z axes or rotates */
/* structure to complete a symmetric structure. */
void c_geometry::reflect( int ix, int iy, int iz, int itx, int nop )
{
	int iti, i, nx, itagi, k;
	nec_float e1, e2, fnop, sam, cs, ss, xk, yk;
	
	np= n;
	mp= m;
	m_ipsym=0;
	iti= itx;
	
	if ( ix >= 0)
	{
		if ( nop == 0)
		return;
	
		m_ipsym=1;
	
		/* reflect along z axis */
		if ( iz != 0)
		{
		m_ipsym=2;
	
		if ( n > 0 )
		{
		/* Reallocate tags buffer */
		segment_tags.resize(2*n + m);
		// segment_tags.resize((2*n+m));
	
		/* Reallocate wire buffers */
		int new_size = 2*n;
		x.resize(new_size);
		y.resize(new_size);
		z.resize(new_size);
		x2.resize(new_size);
		y2.resize(new_size);
		z2.resize(new_size);
		segment_radius.resize(new_size);
	
		for( i = 0; i < n; i++ )
		{
		nx= i+ n;
		e1= z[i];
		e2= z2[i];
	
		if ( (fabs(e1)+fabs(e2) <= 1.0e-5) || (e1*e2 < -1.0e-6) )
		{
			nec_exception* nex = new nec_exception("GEOMETRY DATA ERROR--SEGMENT ");
			nex->append(i+1);
			nex->append("LIES IN PLANE OF SYMMETRY");
			throw nex;
		}
	
		x[nx]= x[i];
		y[nx]= y[i];
		z[nx]=- e1;
		x2[nx]= x2[i];
		y2[nx]= y2[i];
		z2[nx]=- e2;
		itagi= segment_tags[i];
	
		if ( itagi == 0)
			segment_tags[nx]=0;
		if ( itagi != 0)
			segment_tags[nx]= itagi+ iti;
	
		segment_radius[nx]= segment_radius[i];
	
		} /* for( i = 0; i < n; i++ ) */
	
		n= n*2;
		iti= iti*2;
	
		} /* if ( n > 0) */
	
		if ( m > 0 )
		{
		/* Reallocate patch buffers */
		int new_size = 2*m;
		px.resize(new_size);
		py.resize(new_size);
		pz.resize(new_size);
		t1x.resize(new_size);
		t1y.resize(new_size);
		t1z.resize(new_size);
		t2x.resize(new_size);
		t2y.resize(new_size);
		t2z.resize(new_size);
		pbi.resize(new_size);
		psalp.resize(new_size);
	
		for( i = 0; i < m; i++ )
		{
			nx = i+m;
			if ( fabs(pz[i]) <= 1.0e-10)
			{
				nec_exception* nex = new nec_exception("GEOMETRY DATA ERROR--PATCH ");
				nex->append(i+1);
				nex->append("LIES IN PLANE OF SYMMETRY");
				throw nex;
			}
	
			px[nx]= px[i];
			py[nx]= py[i];
			pz[nx]=- pz[i];
			t1x[nx]= t1x[i];
			t1y[nx]= t1y[i];
			t1z[nx]=- t1z[i];
			t2x[nx]= t2x[i];
			t2y[nx]= t2y[i];
			t2z[nx]=- t2z[i];
			psalp[nx]=- psalp[i];
			pbi[nx]= pbi[i];
		}
	
		m= m*2;
	
		} /* if ( m >= m2) */
	
		} /* if ( iz != 0) */
	
		/* reflect along y axis */
		if ( iy != 0)
		{
		if ( n > 0)
		{
		/* Reallocate tags buffer */
		segment_tags.resize(2*n + m);
		// segment_tags.resize((2*n+m));/*????*/
	
		/* Reallocate wire buffers */
		int new_size = 2*n;
		x.resize(new_size);
		y.resize(new_size);
		z.resize(new_size);
		x2.resize(new_size);
		y2.resize(new_size);
		z2.resize(new_size);
		segment_radius.resize(new_size);
	
		for( i = 0; i < n; i++ )
		{
			nx= i+ n;
			e1= y[i];
			e2= y2[i];
		
			if ( (fabs(e1)+fabs(e2) <= 1.0e-5) || (e1*e2 < -1.0e-6) )
			{
				nec_exception* nex = new nec_exception("GEOMETRY DATA ERROR--SEGMENT ");
				nex->append(i+1);
				nex->append("LIES IN PLANE OF SYMMETRY");
				throw nex;
			}
		
			x[nx]= x[i];
			y[nx]=- e1;
			z[nx]= z[i];
			x2[nx]= x2[i];
			y2[nx]=- e2;
			z2[nx]= z2[i];
			itagi= segment_tags[i];
		
			if ( itagi == 0)
				segment_tags[nx]=0;
			if ( itagi != 0)
				segment_tags[nx]= itagi+ iti;
		
			segment_radius[nx]= segment_radius[i];
		
		} /* for( i = n2-1; i < n; i++ ) */
	
		n= n*2;
		iti= iti*2;
	
		} /* if ( n >= n2) */
	
		if ( m > 0 )
		{
			/* Reallocate patch buffers */
			int new_size = 2*m;
			px.resize(new_size);
			py.resize(new_size);
			pz.resize(new_size);
			t1x.resize(new_size);
			t1y.resize(new_size);
			t1z.resize(new_size);
			t2x.resize(new_size);
			t2y.resize(new_size);
			t2z.resize(new_size);
			pbi.resize(new_size);
			psalp.resize(new_size);
		
			for( i = 0; i < m; i++ )
			{
				nx= i+m;
				if ( fabs( py[i]) <= 1.0e-10)
				{
					nec_exception* nex = new nec_exception("GEOMETRY DATA ERROR--PATCH ");