nec_context.h

Go to the documentation of this file.

#pragma once
 
/*
  Copyright (C) 2004-2008,2015  Timothy C.A. Molteno
  tim@molteno.net 
  
  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 <memory>
 
#include "common.h"
#include "c_ggrid.h"
#include "math_util.h"
#include "matrix_algebra.h"
#include "electromag.h"
#include "nec_radiation_pattern.h"
#include "nec_results.h"
#include "nec_structure_currents.h"
#include "nec_output.h"
#include "nec_ground.h"
#include "c_plot_card.h"
 
class c_geometry;
 
enum excitation_return {
    FREQ_PRINT_NORMALIZATION = 0,
    FREQ_LOOP_CONTINUE = 1,
    FREQ_LOOP_CARD_CONTINUE = 2,
    FREQ_INNER_LOOP_CONTINUE = 3,
    FREQ_PHASE_COMPLETE = 4
};
 
enum processing_state {
    PROCESSING_MEMORY_ALLOC = 1,
    PROCESSING_STRUCTURE_LOADING = 2,
    PROCESSING_EXCITATION_SETUP = 3,
    PROCESSING_EXCITATION_LOOP = 4,
    PROCESSING_NEAR_FIELD = 5,
    PROCESSING_FAR_FIELD = 6
};
 
class nec_context
{
public:
  nec_context();
  virtual ~nec_context();
  
  // Called after construction...
  void initialize();
  
 
  void calc_prepare();
  
  void reset_processing_to_structure_loading() {
      switch (processing_state) {
          /* If we havent allocated memory, then leave state there. Otherwise set to
           * structure loading
           * */
          case PROCESSING_MEMORY_ALLOC:
              processing_state = PROCESSING_MEMORY_ALLOC;
              break;
          default:
              processing_state = PROCESSING_STRUCTURE_LOADING;
      }
  }
 
  void reset_processing_to_excitation_setup() {
      switch (processing_state) {
          /* If we havent allocated memory, then leave state there. Otherwise set to
           * structure loading
           * */
          case PROCESSING_MEMORY_ALLOC:
              processing_state = PROCESSING_MEMORY_ALLOC;
              break;
          case PROCESSING_STRUCTURE_LOADING:
              processing_state = PROCESSING_STRUCTURE_LOADING;
              break;
          default:
              processing_state = PROCESSING_EXCITATION_SETUP;
      }
  }
 
  inline c_geometry* get_geometry()  {
      return m_geometry.get();
  }
  
  double get_gain(int freq_index, int theta_index, int phi_index)  {
          nec_radiation_pattern* rp = get_radiation_pattern(freq_index);
          if (NULL == rp)       return -999.0;  
          return rp->get_power_gain(theta_index, phi_index);
  }
 
  double get_gain_max(int freq_index = 0)  {
          nec_radiation_pattern* rp = get_radiation_pattern(freq_index);
          if (NULL == rp)       return -999.0;  
          return rp->get_gain_max();
  }
 
  double get_gain_min(int freq_index = 0)  {
          nec_radiation_pattern* rp = get_radiation_pattern(freq_index);
          if (NULL == rp)  return -999.0;  
          return rp->get_gain_min();
  }
 
  double get_gain_mean(int freq_index = 0)  {
          nec_radiation_pattern* rp = get_radiation_pattern(freq_index);
          if (NULL == rp)  return -999.0;  
          return rp->get_gain_mean();
  }
  
  double get_gain_sd(int freq_index = 0)  {
          nec_radiation_pattern* rp = get_radiation_pattern(freq_index);
          if (NULL == rp)  return -999.0;  
          return rp->get_gain_sd();
  }
  
  /********************** RHCP ********************************/
  double get_gain_rhcp_max(int freq_index = 0)  {
          nec_radiation_pattern* rp = get_radiation_pattern(freq_index);
          if (NULL == rp)  return -999.0;  
          return rp->get_gain_rhcp_max();
  }
 
  double get_gain_rhcp_min(int freq_index = 0)  {
          nec_radiation_pattern* rp = get_radiation_pattern(freq_index);
          if (NULL == rp)  return -999.0;  
          return rp->get_gain_rhcp_min();
  }
 
  double get_gain_rhcp_mean(int freq_index = 0)  {
          nec_radiation_pattern* rp = get_radiation_pattern(freq_index);
          if (NULL == rp)  return -999.0;  
          return rp->get_gain_rhcp_mean();
  }
  
  double get_gain_rhcp_sd(int freq_index = 0)  {
          nec_radiation_pattern* rp = get_radiation_pattern(freq_index);
          if (NULL == rp)  return -999.0;  
          return rp->get_gain_rhcp_sd();
  }
  
  /********************** LHCP ********************************/
  double get_gain_lhcp_max(int freq_index = 0)  {
          nec_radiation_pattern* rp = get_radiation_pattern(freq_index);
          if (NULL == rp)  return -999.0;  
          return rp->get_gain_lhcp_max();
  }
 
  double get_gain_lhcp_min(int freq_index = 0)  {
          nec_radiation_pattern* rp = get_radiation_pattern(freq_index);
          if (NULL == rp)  return -999.0;  
          return rp->get_gain_lhcp_min();
  }
 
  double get_gain_lhcp_mean(int freq_index = 0)  {
          nec_radiation_pattern* rp = get_radiation_pattern(freq_index);
          if (NULL == rp)  return -999.0;  
          return rp->get_gain_lhcp_mean();
  }
  
  double get_gain_lhcp_sd(int freq_index = 0)  {
          nec_radiation_pattern* rp = get_radiation_pattern(freq_index);
          if (NULL == rp)  return -999.0;  
          return rp->get_gain_lhcp_sd();
  }
  
  /****************** IMPEDANCE CHARACTERISTICS *********************/
 
  double get_impedance_real(int freq_index = 0, int feed_index = 0)  {
          nec_antenna_input* ipt = get_input_parameters(freq_index);
          if (NULL == ipt) return -999.0;
          vector<nec_complex>& imp(ipt->get_impedance());
          if (feed_index < 0 || (size_t)feed_index >= imp.size()) return -999.0;
          return imp[feed_index].real();
  }
  double get_impedance_imag(int freq_index = 0, int feed_index = 0)  {
          nec_antenna_input* ipt = get_input_parameters(freq_index);
          if (NULL == ipt) return -999.0;  
          vector<nec_complex>& imp(ipt->get_impedance());
          if (feed_index < 0 || (size_t)feed_index >= imp.size()) return -999.0;
          return imp[feed_index].imag();
  }
 
  inline nec_antenna_input* get_input_parameters(int index)  {
          return m_results.get_antenna_input(index);
  }
  
  inline nec_norm_rx_pattern* get_norm_rx_pattern(int index)  {
          return m_results.get_norm_rx_pattern(index);
  }
  
  inline nec_radiation_pattern* get_radiation_pattern(int index)  {
          return m_results.get_radiation_pattern(index);
  }
          
  inline nec_structure_excitation* get_structure_excitation(int index)  {
          return m_results.get_structure_excitation(index);
  }
  
  inline nec_near_field_pattern* get_near_field_pattern(int index)  {
          return m_results.get_near_field_pattern(index);
  }
  
  inline nec_structure_currents* get_structure_currents(int index)  {
          return m_results.get_structure_currents(index);
  }
  
  /* added for the python wrapping : some access functions */
  
  void set_isave(int in_isave)
  {
          isave = in_isave;
  }
  
  int get_inc()
  {
          return inc;
  }
  
  nec_float get_xpr1()
  {
          return xpr1;
  }
  
  nec_float get_xpr2()
  {
          return xpr2;
  }
  
  /* end of functions added for the python wrapping */
  
  inline void set_output(nec_output_file in_output, nec_output_flags in_output_flags)
  {
          m_output = in_output;
          m_output_flags = in_output_flags;
          
          m_output_fp = m_output.get_fp();
  }
  
  inline void set_results_format(enum RESULT_FORMAT result_format)
  {
          m_results.m_result_format = result_format;
  }
  
  inline void set_gain_only(bool flag)
  {
          m_output_flags.set_gain_only(flag);
  }
  
  
  static nec_float benchmark();
  
  void geometry_complete(int gpflag);
 
  
  void medium_parameters(nec_float permittivity, nec_float permeability)  {
    em::constants::permittivity = permittivity;
    em::constants::permeability = permeability;
  }
  
  
  void 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);
 
  void sp_card(int ns,
      nec_float x1, nec_float y1, nec_float z1,
      nec_float x2, nec_float y2, nec_float z2);
 
  void sc_card( int i2,
      nec_float x3, nec_float y3, nec_float z3,
      nec_float x4, nec_float y4, nec_float z4);
 
  void gx_card(int i1, int i2);
 
  void 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 );
 
  void arc( int tag_id, int segment_count, nec_float rada,
      nec_float ang1, nec_float ang2, nec_float rad );
      
  
  void helix(int tag_id, int segment_count, nec_float s, nec_float hl, nec_float a1, nec_float b1,
      nec_float a2, nec_float b2, nec_float rad);
  
  
  void fr_card(int ifrq, int nfrq, nec_float freq_mhz, nec_float del_freq);
  
  void ld_card(int itmp1, int itmp2, int itmp3, int itmp4, nec_float tmp1, nec_float tmp2, nec_float tmp3);
  
  void gn_card(int ground_type, int rad_wire_count, nec_float tmp1, nec_float tmp2, nec_float tmp3, nec_float tmp4, nec_float tmp5, nec_float tmp6);
  
  
  void ex_card(enum excitation_type itmp1, int itmp2, int itmp3, int itmp4, nec_float tmp1, nec_float tmp2, nec_float tmp3,         
                    nec_float tmp4, nec_float tmp5, nec_float tmp6);
  
 
 
  void ex_card(enum excitation_type itmp1, int itmp2, int itmp3, int itmp4, int itmp5,
                  nec_float tmp1, nec_float tmp2, nec_float tmp3, nec_float tmp4, nec_float tmp5, nec_float tmp6)
  {
    int itmp45 = 10*itmp4 + itmp5;
    return this->ex_card( itmp1, itmp2, itmp3, itmp45, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6 );  
  }
  
  
  
  void tl_card(int itmp1, int itmp2, int itmp3, int itmp4, nec_float tmp1, nec_float tmp2, nec_float tmp3, nec_float tmp4, nec_float tmp5, nec_float tmp6);
  
  void nt_card(int itmp1, int itmp2, int itmp3, int itmp4, nec_float tmp1, nec_float tmp2, nec_float tmp3, nec_float tmp4, nec_float tmp5, nec_float tmp6);
  
  void xq_card(int itmp1);
  
  void gd_card(nec_float tmp1, nec_float tmp2, nec_float tmp3, nec_float tmp4);
  
  void rp_card(int calc_mode,
    int n_theta, int n_phi,
    int output_format, int normalization, int D, int A,
    nec_float theta0, nec_float phi0, nec_float delta_theta, nec_float delta_phi,
    nec_float radial_distance, nec_float gain_norm);
  
  void pt_card(int itmp1, int itmp2, int itmp3, int itmp4);
  
  
  void pq_card(int itmp1, int itmp2, int itmp3, int itmp4);
  
  
  
  void kh_card(nec_float tmp1);
  
  
  void ne_card(int itmp1, int itmp2, int itmp3, int itmp4, nec_float tmp1, nec_float tmp2, nec_float tmp3, nec_float tmp4, nec_float tmp5, nec_float tmp6);
  
  void nh_card(int itmp1, int itmp2, int itmp3, int itmp4, nec_float tmp1, nec_float tmp2, nec_float tmp3, nec_float tmp4, nec_float tmp5, nec_float tmp6);
  
  void set_extended_thin_wire_kernel(bool ekflag);
  
  
  void cp_card(int itmp1, int itmp2, int itmp3, int itmp4);
  
  
  void pl_card(const char* ploutput_filename, int itmp1, int itmp2, int itmp3, int itmp4);
  
  
  /*!****************************************************
  *** normal exit of nec2++ when all jobs complete ok ***
  ******************************************************/
  void all_jobs_completed()
  {
  }
  
  void write_results(ostream& os)
  {
    m_results.write(os);
  }
  
  
  
  void simulate(bool far_field_flag = false);
 
  nec_output_file m_output;
  nec_ground ground;
  std::unique_ptr<c_geometry> m_geometry;
  c_plot_card plot_card;
  
  c_ground_wave ground_wave;
 
  
  int iptflq;
  int iptaq, iptaqf, iptaqt;
  
  
  int iptflg;
  int iptag, iptagf, iptagt;
  
  
  int iflow;
  int ifrq, nfrq;
  nec_float delfrq;
  
  // strcture loading
  int_array ldtyp, ldtag, ldtagf, ldtagt;
  real_array zlr, zli, zlc;
  
  // normalized receiving pattern
  real_matrix fnorm;
  
  int nthi, nphi;
  nec_float thetis, phiss;
  
  
  nec_results m_results;
  
  
  nec_output_flags m_output_flags;
  
  
  nec_float _wavelength;
  
  /* common  /cmb/ */
  complex_array cm;  // primary interaction matrix
  
  /* common  /matpar/ */
  int icase, npblk;
  int64_t nlast;
  int nbbx, npbx, nlbx, nbbl, npbl, nlbl;
  
  /* common  /save/ */
  int_array ip;
  nec_float freq_mhz;
  
  /* common  /crnt/ */
  real_array air, aii; 
  real_array bir, bii; 
  real_array cir, cii; 
  complex_array current_vector;  
  
  int ifar; 
    
  /* common  /zload/ */
  int nload;
  complex_array zarray;
  
  /* common  /yparm/ */
  int ncoup, icoup;
  int_array nctag, ncseg;
  complex_array y11a, y12a;
  
  /* common  /vsorc/ */
  int_array ivqd, source_segment_array, iqds;
  int nvqd, voltage_source_count, nqds;
  complex_array vqd, vqds, source_voltage_array;
  
  /* common  /netcx/ */
  int masym, neq2, network_count, ntsol, nprint;
  int64_t neq, npeq;
  int_array iseg1, iseg2, ntyp;
  real_array x11r, x11i, x12r;
  real_array x12i, x22r, x22i;
  nec_float input_power, network_power_loss;
  nec_complex zped;
  
  /* common  /fpat/ */
  enum excitation_type m_excitation_type;
  
  int m_rp_output_format;
  int m_rp_normalization;
  
  int m_near, nfeh, nrx, nry, nrz, nth, nph, ipd, iavp;
  nec_float thets, phis, dth, dph, rfld, gnor; 
  nec_float xpr6, structure_power_loss, xnr, ynr, znr, dxnr, dynr, dznr;
  
  
  /* common  /dataj/ */
  int ind1, indd1, ind2, indd2;
  bool m_use_exk; /* Was iexk */
  
  nec_float m_s, m_b, xj, yj, zj, cabj, sabj, salpj;
  nec_float rkh; /* matrix integration limit */
  nec_float t1xj, t1yj, t1zj, t2xj, t2yj, t2zj;
  nec_complex  exk, eyk, ezk, exs, eys, ezs, exc, eyc, ezc;
  
  /* common  /smat/ */
  int nop; /* My addition */
  complex_array symmetry_array;
  
  /* common  /incom/ */
  int isnor;
  nec_float xo, yo, zo, sn, xsn, ysn;
  
  /* common  /tmi/ */
  int ija; /* changed to ija to avoid conflict */
  nec_float zpk, rkb2;
  
  /*common  /tmh/ */
  nec_float zpka, rhks;
 
  
  // some auxiliary functions to be made private once
  // the radiation pattern calculation is done entirely
  // inside this class...
  void gfld(nec_float rho, nec_float phi, nec_float rz,
    nec_complex *eth, nec_complex *epi,
    nec_complex *erd, bool space_only, nec_float in_wavelength  );
  
  void ffld(nec_float thet, nec_float phi,
    nec_complex *eth, nec_complex *eph, nec_float in_wavelength  );
  
 
private:
 
  bool step_warning_issued;
 
  void ne_nh_card(int in_nfeh, int itmp1, int itmp2, int itmp3, int itmp4, nec_float tmp1, nec_float tmp2, nec_float tmp3, nec_float tmp4, nec_float tmp5, nec_float tmp6);
  
  
  void print_freq_int_krnl(
    nec_float f, 
    nec_float lambda, 
    nec_float int_dist, 
    bool using_extended_kernel);
    
  void   antenna_env(void);
  
  /*no more used */
  void  print_structure_currents(char *pattype, int iptflg, int iptflq, int iptag, int iptagf, int iptagt, int iptaq, 
    int iptaqf, int iptaqt);
    
  void  calculate_network_data(void);
  void  print_network_data(void);
  void  print_norm_rx_pattern();
  void  print_input_impedance();
  void  print_power_budget(void);
  nec_float compute_structure_power_loss(void) const;
  void  structure_segment_loading();
    
  
  enum excitation_return
    excitation_loop(enum processing_state in_freq_loop_state, int mhz);
 
  excitation_return excitation_process_inner(int mhz);
  excitation_return excitation_compute_near_field(int mhz);
  void excitation_compute_far_field();
 
  void  setup_excitation();
  
  /* pointers to output files */
  FILE *m_output_fp;
  
  int inc, isave;
  enum processing_state processing_state;
  int nthic, nphic;
  int iped;
  
  nec_float impedance_norm_factor; // was zpnorm
  
  nec_float xpr1, xpr2, xpr3, xpr4, xpr5, xpr7;
  
  nec_structure_currents* structure_currents;
  
  void load();
 
  void cmset(int64_t nrow, complex_array& in_cm, nec_float rkhx);
  void compute_matrix_ss(int j1, int j2, int im1, int im2,
      complex_array& in_cm, int64_t nrow, int itrp);
  void cmsw(int j1, int j2, int i1, int i2, complex_array& in_cm,
      complex_array& cw, int64_t ncw, int64_t nrow, int itrp);
  void cmws( int j, int i1, int i2, complex_array& in_cm,
                        int64_t nr, complex_array& cw, int64_t nw, int itrp );
        
  void cmww(int j, int i1, int i2, complex_array& in_cm, int64_t nr,
      complex_array& cw, int64_t nw, int itrp);
  void couple(complex_array& cur, nec_float wlam);
 
  void efld(nec_float xi, nec_float yi, nec_float zi, nec_float ai, bool not_on_source_segment);
  nec_float efld_compute(nec_float xij, nec_float yij, nec_float ai, nec_float salpr,
      nec_float zij, bool ijx, complex_array& egnd);
  void eksc(nec_float s, nec_float z, nec_float rh, nec_float xk, int ij,
      nec_complex *ezs, nec_complex *ers, nec_complex *ezc,
      nec_complex *erc, nec_complex *ezk, nec_complex *erk);
  void ekscx(nec_float bx, nec_float s, nec_float z, nec_float rhx, nec_float xk,
      int ij, int inx1, int inx2, nec_complex *ezs,
      nec_complex *ers, nec_complex *ezc, nec_complex *erc,
      nec_complex *ezk, nec_complex *erk);
          void 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);
        
          void etmns_voltage_source(complex_array& e);
          void etmns_linear_wave(nec_float cth, nec_float sth, nec_float cph, nec_float sph,
              nec_float pxl, nec_float pyl, nec_float pzl, nec_float wx, nec_float wy, nec_float wz,
              nec_float qx, nec_float qy, nec_float qz, nec_complex rrv, nec_complex rrh,
              nec_float incident_amplitude, complex_array& e);
          void etmns_circular_wave(nec_float cth, nec_float sth, nec_float cph, nec_float sph,
              nec_float pxl, nec_float pyl, nec_float pzl, nec_float wx, nec_float wy, nec_float wz,
              nec_float qx, nec_float qy, nec_float qz, nec_complex rrv, nec_complex rrh,
              nec_float p6, nec_float incident_amplitude, enum excitation_type excite_type, complex_array& e);
          void etmns_current_source(nec_float p1, nec_float p2, nec_float p3, nec_float p4,
              nec_float p5, nec_float p6, complex_array& e);
          int etmns_patch_base(int patch_index) const;
          void etmns_zero_element(complex_array& e, int i);
        
          void fblock( int nrow, int ncol, int64_t imax, int ipsym );
 
  void gf(nec_float zk, nec_float *co, nec_float *si);
  void gh(nec_float zk, nec_float *hr, nec_float *hi);
  void gx(nec_float zz, nec_float rh, nec_float xk,
      nec_complex *gz, nec_complex *gzp);
  void gxx(nec_float zz, nec_float rh, nec_float a, nec_float a2, nec_float xk,
      int ira, nec_complex *g1, nec_complex *g1p, nec_complex *g2,
      nec_complex *g2p, nec_complex *g3, nec_complex *gzp);
  void hfk(nec_float el1, nec_float el2, nec_float rhk,
      nec_float zpkx, nec_float *sgr, nec_float *sgi);
  void hintg(nec_float xi, nec_float yi, nec_float zi);
  void hsfld(nec_float xi, nec_float yi, nec_float zi, nec_float ai);
  void hsflx(nec_float s, nec_float rh, nec_float zpx, nec_complex *hpk,
      nec_complex *hps, nec_complex *hpc);
 
  void intx(nec_float el1, nec_float el2, nec_float b, int ij,
      nec_float *sgr, nec_float *sgi);
 
  void nefld(nec_float xob, nec_float yob, nec_float zob, nec_complex *ex,
      nec_complex *ey, nec_complex *ez);
  struct netwk_state {
    int_array  ipnt, nteqa, ntsca;
    complex_array vsrc, rhs, cmn, rhnt, rhnx;
    int nteq  = 0;
    int ntsc  = 0;
    int ndimn = 0;
    int neqz2 = 0;
    int neqt  = 0;
  };
  void netwk(complex_array& in_cm, int_array& in_ip, complex_array& einc);
  void netwk_build_matrix(complex_array& in_cm, int_array& in_ip, netwk_state& st);
  void netwk_solve_and_check(complex_array& in_cm, int_array& in_ip, netwk_state& st);
  void netwk_compute_currents(complex_array& in_cm, int_array& in_ip, complex_array& einc, netwk_state& st);
  void netwk_compute_inputs(complex_array& einc, netwk_state& st);
  void nfpat(void);
  void nhfld(nec_float xob, nec_float yob, nec_float zob, nec_complex *hx,
      nec_complex *hy, nec_complex *hz);
  void pcint(nec_float xi, nec_float yi, nec_float zi, nec_float cabi,
      nec_float sabi, nec_float salpi, complex_array& e);
  void impedance_print(int in1, int in2, int in3, nec_float fl1, nec_float fl2,
      nec_float fl3, nec_float fl4, nec_float fl5, nec_float fl6, const char *ia);
  void qdsrc(int is, nec_complex v, complex_array& e);
 
  
  void rom2(nec_float a, nec_float b, complex_array& sum, nec_float dmin);
  void sflds(const nec_float t, complex_array& e);
  void solgf(nec_complex *a, nec_complex *b, nec_complex *c,
      nec_complex *d, nec_complex *xy, int *ip, int np, int n1,
      int n, int mp, int m1, int m, int n1c, int n2c, int n2cz);
  void unere(nec_float xob, nec_float yob, nec_float zob, bool ground_reflection);
  nec_complex zint(nec_float sigl, nec_float rolam);
  
  void init_voltage_sources();
 
}; /* nec_context */