Logo Search packages:      
Sourcecode: necpp version File versions  Download package

void nec_context::etmns ( nec_float  p1,
nec_float  p2,
nec_float  p3,
nec_float  p4,
nec_float  p5,
nec_float  p6,
nec_float  incident_amplitude,
enum excitation_type  excite_type,
complex_array e 
) [private]

Fills the array e with the negative of the electric field tangent to the segments and the tangential magnetic field on the surfaces.

Parameters:
e is the right hand side of the matrix equation.

Definition at line 3638 of file nec_context.cpp.

References c_geometry::cab, safe_array< T >::fill(), c_geometry::m, c_geometry::n, c_geometry::n_plus_m, c_geometry::patch_angle(), c_geometry::psalp, c_geometry::px, c_geometry::py, c_geometry::pz, c_geometry::sab, c_geometry::salp, c_geometry::segment_length, c_geometry::t1x, c_geometry::t1y, c_geometry::t1z, c_geometry::t2x, c_geometry::t2y, c_geometry::t2z, c_geometry::x, c_geometry::y, and c_geometry::z.

{
      nec_float cth, sth, cph, sph, cet, set, pxl, pyl, pzl, wx;
      nec_float wy, wz, qx, qy, qz, ds, dsh, rs, r;
      nec_complex er, et, ezh, erh, rrv, rrh, tt1, tt2;
      
      int n = m_geometry->n;
      int m = m_geometry->m;
      
//    int neq= n+2*m;
      ASSERT(neq == n+2*m);
      
      nqds=0;
      
      ASSERT(excite_type >= 0);
            
      /* applied field of voltage sources for transmitting case */
      if ( (excite_type == EXCITATION_VOLTAGE) || (excite_type == EXCITATION_VOLTAGE_DISC) )
      {
            e.fill(0,neq,cplx_00());
      
            for (int i = 0; i < voltage_source_count; i++ )
            {
                  int source_index = source_segment_array[i]-1;
                  e[source_index] = -source_voltage_array[i]/(m_geometry->segment_length[source_index]* wavelength);
            }
      
            if ( nvqd == 0)
                  return;
      
            for (int i = 0; i < nvqd; i++ )
            {
                  qdsrc( ivqd[i], vqd[i], e);
            }
            return;
      } /* if ( (excite_type == EXCITATION_VOLTAGE) || (excite_type == EXCITATION_VOLTAGE_DISC) ) */

      /* incident plane wave, linearly polarized. */
      if ((excite_type == EXCITATION_LINEAR) ||
            (excite_type == EXCITATION_CIRC_RIGHT) ||
            (excite_type == EXCITATION_CIRC_LEFT))
      {
            cth= cos( p1);
            sth= sin( p1);
            cph= cos( p2);
            sph= sin( p2);
            cet= cos( p3);
            set= sin( p3);
            pxl= cth* cph* cet- sph* set;
            pyl= cth* sph* cet+ cph* set;
            pzl=- sth* cet;
            wx=- sth* cph;
            wy=- sth* sph;
            wz=- cth;
            qx= wy* pzl- wz* pyl;
            qy= wz* pxl- wx* pzl;
            qz= wx* pyl- wy* pxl;

            if (ground.present())
            {
                  if (ground.type_perfect())
                  {
                        rrv=-cplx_10();
                        rrh=-cplx_10();
                  }
                  else
                  {
                        rrv= sqrt(1.0 - ground.get_zrati_sqr() * sth*sth);
                        rrh= ground.zrati* cth;
                        rrh=( rrh- rrv)/( rrh+ rrv);
                        rrv= ground.zrati* rrv;
                        rrv=-( cth- rrv)/( cth+ rrv);
                  }
            }
      
            if ( excite_type == EXCITATION_LINEAR)
            {
                  if ( n != 0)
                  {
                        for (int i = 0; i < n; i++ )
                        {
                              nec_float arg=- two_pi() *( wx* m_geometry->x[i]+ wy* m_geometry->y[i]+ wz* m_geometry->z[i]);
                              nec_complex e_amplitude(cos(arg), sin(arg));
                              
                              e[i] = -(pxl*m_geometry->cab[i]+ pyl*m_geometry->sab[i]+ pzl*m_geometry->salp[i]) * e_amplitude * incident_amplitude;
                        }
                  
                        if (ground.present())
                        {
                              tt1=( pyl* cph- pxl* sph)*( rrh- rrv);
                              nec_complex cx= rrv* pxl- tt1* sph;
                              nec_complex cy= rrv* pyl+ tt1* cph;
                              nec_complex cz=- rrv* pzl;
                        
                              for (int i = 0; i < n; i++ )
                              {
                                    nec_float arg=- two_pi()*( wx* m_geometry->x[i]+ wy* m_geometry->y[i]- wz* m_geometry->z[i]);
                                    nec_complex e_amplitude(cos(arg), sin(arg));
                                    e[i] -= ( cx* m_geometry->cab[i]+ cy* m_geometry->sab[i]+ cz* m_geometry->salp[i])* e_amplitude * incident_amplitude;
                              }
                        }
                  } /* if ( n != 0) */
      
                  if ( m == 0) // if no surface patches we're done!
                        return;
      
                  {
                        int i= -1;
                        int i1= n-2;
                        for (int is = 0; is < m; is++ )
                        {
                              i++;
                              ASSERT(is == i);
                              i1 += 2;
                              int i2 = i1+1;
                              nec_float arg = -m_geometry->patch_angle(i,wx,wy,wz);
                              ASSERT( arg == (- two_pi()*( wx* m_geometry->px[i]+ wy* m_geometry->py[i]+ wz* m_geometry->pz[i])));
                              
                              nec_complex e_amplitude = nec_complex(cos(arg), sin(arg)) * m_geometry->psalp[i] * em::inverse_impedance();
                              e[i2]=( qx* m_geometry->t1x[i]+ qy* m_geometry->t1y[i]+ qz* m_geometry->t1z[i])* e_amplitude * incident_amplitude;
                              e[i1]=( qx* m_geometry->t2x[i]+ qy* m_geometry->t2y[i]+ qz* m_geometry->t2z[i])* e_amplitude * incident_amplitude;
                        }
                  }
                  
                  if (ground.present())
                  {
                        tt1=( qy* cph- qx* sph)*( rrv- rrh);
                        nec_complex cx=-( rrh* qx- tt1* sph);
                        nec_complex cy=-( rrh* qy+ tt1* cph);
                        nec_complex cz= rrh* qz;
                        
                        int i= -1;
                        int i1= n-2;
                        for (int is = 0; is < m; is++ )
                        {
                              i++;
                              i1 += 2;
                              int i2 = i1+1;
                              nec_float arg = -m_geometry->patch_angle(i,wx,wy,wz);
                              ASSERT((-two_pi()*( wx* m_geometry->px[i]+ wy* m_geometry->py[i]- wz* m_geometry->pz[i])) == arg);
                              
                              nec_complex e_amplitude = cplx_exp(arg) * m_geometry->psalp[i] * em::inverse_impedance();
                              
                              e[i2]= e[i2]+( cx* m_geometry->t1x[i]+ cy* m_geometry->t1y[i]+ cz* m_geometry->t1z[i])* e_amplitude * incident_amplitude;
                              e[i1]= e[i1]+( cx* m_geometry->t2x[i]+ cy* m_geometry->t2y[i]+ cz* m_geometry->t2z[i])* e_amplitude * incident_amplitude;
                        }
                  }
                  return;
            } /* if ( excite_type == EXCITATION_LINEAR) */

            /* incident plane wave, elliptic polarization. */
            tt1=-(cplx_01())* p6;
            if ( excite_type == EXCITATION_CIRC_LEFT)
                  tt1=- tt1;
      
            if ( n != 0)
            {
                  nec_complex cx= pxl+ tt1* qx;
                  nec_complex cy= pyl+ tt1* qy;
                  nec_complex cz= pzl+ tt1* qz;
            
                  for (int i = 0; i < n; i++ )
                  {
                        nec_float arg=- two_pi()*( wx* m_geometry->x[i]+ wy* m_geometry->y[i]+ wz* m_geometry->z[i]);
                        nec_complex e_amplitude(cos(arg), sin(arg));
                        e[i] = -(cx* m_geometry->cab[i] + cy* m_geometry->sab[i] + cz*m_geometry->salp[i])* e_amplitude * incident_amplitude;
                  }
                  
                  
                  if (ground.present())
                  {
                        tt2=( cy* cph- cx* sph)*( rrh- rrv);
                        nec_complex ccx= rrv* cx- tt2* sph;
                        nec_complex ccy= rrv* cy+ tt2* cph;
                        nec_complex ccz=- rrv* cz;
                  
                        for (int i = 0; i < n; i++ )
                        {
                              nec_float arg=- two_pi()*( wx* m_geometry->x[i]+ wy* m_geometry->y[i]- wz* m_geometry->z[i]);
                              nec_complex e_amplitude(cos(arg), sin(arg));
                              e[i] -= (ccx* m_geometry->cab[i]+ ccy* m_geometry->sab[i]+ ccz* m_geometry->salp[i]) * e_amplitude * incident_amplitude;
                        }
                  }
            } /* if ( n != 0) */

            if ( m == 0)
                  return;

            nec_complex cx= qx- tt1* pxl;
            nec_complex cy= qy- tt1* pyl;
            nec_complex cz= qz- tt1* pzl;
      
            {
                  int i= -1;
                  int i1= n-2;
                  for (int is = 0; is < m; is++ )
                  {
                        i++;
                        i1 += 2;
                        int i2 = i1+1;
                        
                        nec_float arg = -m_geometry->patch_angle(i,wx,wy,wz);
                        ASSERT(arg == -two_pi()*( wx* m_geometry->px[i]+ wy* m_geometry->py[i]+ wz* m_geometry->pz[i]));

                        nec_complex e_amplitude = nec_complex(cos(arg), sin(arg)) * m_geometry->psalp[i] * em::inverse_impedance();
                        
                        e[i2]=( cx* m_geometry->t1x[i]+ cy* m_geometry->t1y[i]+ cz* m_geometry->t1z[i])* e_amplitude * incident_amplitude;
                        e[i1]=( cx* m_geometry->t2x[i]+ cy* m_geometry->t2y[i]+ cz* m_geometry->t2z[i])* e_amplitude * incident_amplitude;
                  }
            }
            
            if (ground.present())
            {
                  tt1=( cy* cph- cx* sph)*( rrv- rrh);
                  cx=-( rrh* cx- tt1* sph);
                  cy=-( rrh* cy+ tt1* cph);
                  cz= rrh* cz;

                  int i= -1;
                  int i1= n-2;
                  for (int is=0; is < m; is++ )
                  {
                        i++;
                        i1 += 2;
                        int i2 = i1+1;
                        nec_float arg = -m_geometry->patch_angle(i,wx,wy,wz);
                        ASSERT(arg == -two_pi()*( wx* m_geometry->px[i]+ wy* m_geometry->py[i]- wz* m_geometry->pz[i]));
                        
                        nec_complex e_amplitude = nec_complex(cos(arg), sin(arg)) * m_geometry->psalp[i] * em::inverse_impedance();
                        
                        e[i2] += ( cx* m_geometry->t1x[i]+ cy* m_geometry->t1y[i]+ cz* m_geometry->t1z[i])* e_amplitude * incident_amplitude;
                        e[i1] += ( cx* m_geometry->t2x[i]+ cy* m_geometry->t2y[i]+ cz* m_geometry->t2z[i])* e_amplitude * incident_amplitude;
                  }
            }
      
            return;
      
      } /* if ( excite_type <= 3) */

      ASSERT(excite_type == EXCITATION_CURRENT);
      /* incident field of an elementary current source. */

      nec_float cosp4 = cos(p4);
      
      wx = cosp4 * cos(p5);
      wy = cosp4 * sin(p5);
      wz = sin(p4);
      ds = p6*em::impedance_over_2pi();
      dsh= p6/(2.0 * two_pi());
            
      // TODO: Split the following loop up into two loops i=1:n and i=1:m
      // to loop over the patches and the segments seperately
      
      for (int i = 0; i < m_geometry->n_plus_m; i++ )
      {
            if ( i < n )
            {
                  pxl = m_geometry->x[i] - p1;
                  pyl = m_geometry->y[i] - p2;
                  pzl = m_geometry->z[i] - p3;
            }
            else  // patches
            {
                  int patch_index = i - n;
                  pxl = m_geometry->px[patch_index] - p1;
                  pyl = m_geometry->py[patch_index] - p2;
                  pzl = m_geometry->pz[patch_index] - p3;
            }
            
            rs = norm2(pxl,pyl,pzl);
            if ( rs < 1.0e-30)
                  continue;
      
            r = sqrt(rs);
            pxl = pxl/r;
            pyl = pyl/r;
            pzl = pzl/r;
            cth = pxl*wx + pyl*wy + pzl*wz;
            sth = sqrt(1.0 - cth*cth);
            qx = pxl - wx*cth;
            qy = pyl - wy*cth;
            qz = pzl - wz*cth;
      
            nec_float arg = norm(qx,qy,qz);
            if ( arg >= 1.e-30)
            {
                  qx = qx / arg;
                  qy = qy / arg;
                  qz = qz / arg;
            }
            else
            {
                  qx = 1.0;
                  qy = 0.0;
                  qz = 0.0;
            } /* if ( arg >= 1.e-30) */
      
            arg = two_pi() * r;
            tt1 = nec_complex(cos(arg), -sin(arg));
      
            if ( i < n )
            {
                  tt2 = nec_complex(1.0,-1.0/arg) / rs;
                  er = ds* tt1* tt2* cth;
                  et = 0.5 * ds * tt1 *((cplx_01()) * two_pi()/r + tt2)* sth;
                  ezh= er*cth - et*sth;
                  erh= er*sth + et*cth;
                  nec_complex cx = ezh*wx + erh*qx;
                  nec_complex cy = ezh*wy + erh*qy;
                  nec_complex cz = ezh*wz + erh*qz;
                  e[i] = -(cx* m_geometry->cab[i]+ cy* m_geometry->sab[i]+ cz* m_geometry->salp[i]);
            }
            else // patches
            {
                  int patch_index = i - n;
                  int i1 = i + patch_index*2; // was i1 += 2; and starts at n-2
                  int i2 = i1+1;
                  
                  pxl = wy*qz - wz*qy; // cross product here...
                  pyl = wz*qx - wx*qz;
                  pzl = wx*qy - wy*qx;
                  
                  tt2= dsh* tt1* nec_complex(1.0/r, two_pi())/ r* sth* m_geometry->psalp[patch_index];
                  nec_complex cx= tt2* pxl;
                  nec_complex cy= tt2* pyl;
                  nec_complex cz= tt2* pzl;
                  e[i2] = cx*m_geometry->t1x[patch_index] + cy*m_geometry->t1y[patch_index] + cz*m_geometry->t1z[patch_index];
                  e[i1] = cx*m_geometry->t2x[patch_index] + cy*m_geometry->t2y[patch_index] + cz*m_geometry->t2z[patch_index];
            } /* if ( i < n) */
      } /* for( i = 0; i < m_geometry->n_plus_m; i++ ) */
}


Generated by  Doxygen 1.6.0   Back to index