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

void nec_context::rom2 ( nec_float  a,
nec_float  b,
complex_array sum,
nec_float  dmin 
) [private]

For the Sommerfeld ground option, rom2 integrates over the source segment to obtain the total field due to ground. The method of variable interval width Romberg integration is used. There are 9 field components - the x, y, and z components due to constant, sine, and cosine current distributions.

Parameters:
a lower limit of integral
b upper limit of integral
dmin Set in EFLD to 1% of the magnitude of the electric field

Definition at line 6187 of file nec_context.cpp.

References safe_array< T >::fill(), m_output, c_geometry::n_plus_m, and safe_array< T >::size().

Referenced by efld().

{
      ASSERT(sum.size() == 9);
      
      bool flag = true;
      static bool step_warning_issued = false;
      
      int nts = 4, nx = 1, n = 9;
      nec_float dz=0., dzot=0.0;
      nec_float rx = 1.0e-4;
      
      complex_array g1(9), g2(9), g3(9), g4(9), g5(9);
      complex_array t01(9), t10(9), t20(9);
      
      nec_float _z = a;
      nec_float ze = b;
      nec_float _s = b - a;
      
      if ( _s < 0.0)
      {
            throw new nec_exception("ERROR - B LESS THAN A IN ROM2");
      }
      
      nec_float ep = _s/(1.0e4 * m_geometry->n_plus_m);
      nec_float zend = ze - ep;
      
      sum.fill(0,n,cplx_00());
      
      int ns = nx;
      int nt = 0;
      sflds( _z, g1);

      while ( true )
      {
            if ( flag )
            {
                  dz = _s / ns;
                  if ( _z + dz >= ze)
                  {
                        dz = ze- _z;
                        if ( dz <= ep)
                              return;
                  }
            
                  dzot= dz*.5;
                  sflds( _z + dzot, g3);
                  sflds( _z + dz, g5);
            }
      
            // Evaluate 3-point Romberg result and test convergence. 
            for (int i = 0; i < n; i++ )
            {
                  nec_complex t00 = (g1[i]+ g5[i]) * dzot;
                  t01[i] = (t00 + dz * g3[i]) * 0.5;
                  t10[i] = (4.0 * t01[i] - t00)/3.;
            }
            
            nec_float tmag1 = sqrt( norm(t01[0]) + norm(t01[1]) + norm(t01[2]) );
            nec_float tmag2 = sqrt( norm(t10[0]) + norm(t10[1]) + norm(t10[2]) );

            nec_float tr = test_simple( tmag1, tmag2, dmin);
            
            if ( tr <= rx)
            {
                  for (int i = 0; i < n; i++ )
                        sum[i] += t10[i];
                  nt += 2;
            
                  // Now we do some housekeeping before looping back
                  _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;
                  }
                  
                  flag = true;
                  continue;
            
            } /* if ( tr <= rx) */
      
            sflds( _z+ dz*.25, g2);
            sflds( _z+ dz*.75, g4);
            
            tmag1 = 0.0;
            tmag2 = 0.0;
      
            /* Evaluate 5 point Romberg result and test convergence. */
            for (int i = 0; i < n; i++ )
            {
                  nec_complex t02 = (t01[i]+ dzot*( g2[i]+ g4[i]))*0.5;
                  nec_complex t11 = (4.0 * t02- t01[i] )/3.0;
                  t20[i] = (16.* t11- t10[i])/15.0;
                  
                  if (i <= 2)
                  {
                        tmag1 += norm(t11);
                        tmag2 += norm(t20[i]);
                  }
            }
      
            tmag1 = sqrt(tmag1);
            tmag2 = sqrt(tmag2);
            tr = test_simple( tmag1, tmag2, dmin);

            // check this with the FORTRAN flow, Neoklis (probably correctly) feels that it
            // should be if (tr > rx)
            // TODO test this code here on Sommerfeld Norton Conditions
            if ( tr > rx)
            {
                  nt=0;
                  if ( ns < m_geometry->n_plus_m )
                  {
                        ns = ns*2;
                        dz = _s/ ns;
                        dzot = dz*.5;
                  
                        for (int i = 0; i < n; i++ )
                        {
                              g5[i] = g3[i];
                              g3[i] = g2[i];
                        }
                  
                        flag=false;
                        continue;
                  }
            
                  nec_error_mode em(m_output);
                  m_output.string("ROM2 -- STEP SIZE LIMITED AT Z = ");
                  m_output.real_out(12,5,_z);
                  m_output.endl();
                  
                  if (false == step_warning_issued)
                  {
                        m_output.line("About the above warning:");
                        m_output.line("Probably caused by a wire too close to the ground in the Somerfeld/");
                        m_output.line("Norton ground method.  Execution continues but results may be inaccurate.");
                        step_warning_issued = true;
                  }
            } /* if ( tr <= rx) */
      
            for (int i = 0; i < n; i++ )
                  sum[i] += t20[i];
            nt = nt+1;
      
            // Now we do the same housekeeping before looping back
            _z = _z + dz;
            if ( _z > zend)
                  return;
      
            for (int i = 0; i < n; i++ )
                  g1[i] = g5[i];
      
            flag = true;
            
            if ( (nt >= nts) && (ns > nx) )
            {
                  /* Double step size */
                  ns = ns/2;
                  nt = 1;
            }
      } /* while( true ) */
}


Generated by  Doxygen 1.6.0   Back to index