PyPO User Manual
Propagation.h
Go to the documentation of this file.
1 #include <iostream>
2 #include <vector>
3 #include <complex>
4 #include <array>
5 #include <cmath>
6 #include <thread>
7 #include <iomanip>
8 #include <algorithm>
9 #include <new>
10 
11 #include "Utils.h"
12 #include "Structs.h"
13 #include "InterfaceReflector.h"
14 
15 #ifndef __Propagation_h
16 #define __Propagation_h
17 
18 /*! \file Propagation.h
19  \brief Functions for PO calculations on CPU.
20 
21  Provides functions for performing PO calculations on the CPU.
22 */
23 
24 /**
25  * Main class for running PO calculations.
26  *
27  * Contains functions that run PO calculations. All calculations can be performed in a parallel way.
28  *
29  * @see Utils
30  * @see Structs
31  * @see InterfaceReflector
32  */
33 template <class T, class U, class V, class W>
34 class Propagation
35 {
36  T k; /**<Wavenumber of radiation, double/float.*/
37  int numThreads; /**<Number of computing threads to employ.*/
38  int gs; /**<Number of cells on source surface.*/
39  int gt; /**<Number of cells on target surface.*/
40 
41  int step; /**<Number of threads per block.*/
42 
43  T t_direction; /**<Time direction (experimental!).*/
44 
45  T EPS; /**<Relative electric permittivity of source medium, double/float.*/
46  T prefactor;
47  T C_L; /**<Speed of light in vacuum, in mm / s.*/
48  T MU_0; /**<Magnetic permeability.*/
49  T EPS_VAC; /**<Vacuum electric permittivity.*/
50  T ZETA; /**<Impedance of surrounding medium.*/
51  T ZETA_INV; /**<Conductance of surrounding medium.*/
52  T PIf; /**<Template pi (redundant?).*/
53  T Three; /**<Template 3*/
54 
55 
56  std::complex<T> j; /**<Complex unit.*/
57  std::complex<T> z0; /**<Complex zero.*/
58 
59  std::array<std::array<T, 3>, 3> eye; /**<Real-valued 3x3 unit matrix.*/
60 
61  void joinThreads();
62 
63  void _debugArray(T *arr, int idx);
64  void _debugArray(std::array<T, 3> arr);
65  void _debugArray(std::array<std::complex<T>, 3> arr);
66 
67 public:
68 
69  std::vector<std::thread> threadPool; /**<Vector of thread objects.*/
70 
71  Propagation(T omega, int numThreads, int gs, int gt, T epsilon, T t_direction, bool verbose = false);
72 
73  // Make T precision utility kit
74  Utils<T> ut; /**<Utils object for vector operations.*/
75 
76  // Functions for propagating fields between two surfaces
77  void propagateBeam_JM_PEC(int start, int stop,
78  V *cs, V *ct,
79  W *currents,
80  int gmode, int ncx, int ncy,
81  U *res);
82 
83  void propagateBeam_JM(int start, int stop,
84  V *cs, V *ct,
85  W *currents,
86  int gmode, int ncx, int ncy,
87  U *res);
88 
89  void propagateBeam_EH(int start, int stop,
90  V *cs, V *ct,
91  W *currents,
92  int gmode, int ncx, int ncy,
93  U *res);
94 
95  void propagateBeam_JMEH_PEC(int start, int stop,
96  V *cs, V *ct,
97  W *currents,
98  int gmode, int ncx, int ncy,
99  U *res);
100 
101  void propagateBeam_JMEH(int start, int stop,
102  V *cs, V *ct,
103  W *currents,
104  int gmode, int ncx, int ncy,
105  U *res);
106 
107  void propagateBeam_EHP(int start, int stop,
108  V *cs, V *ct,
109  W *currents,
110  int gmode, int ncx, int ncy,
111  U *res);
112 
113  std::array<std::array<std::complex<T>, 3>, 2> fieldAtPoint(V *cs, W *currents,
114  int gmode, int ncx, int ncy,
115  const std::array<T, 3> &point_target);
116 
117  std::array<std::array<std::complex<T>, 3>, 2> fieldAtPoint_Old(V *cs, W *currents,
118  const std::array<T, 3> &point_target);
119 
120 
121  void parallelProp_JM(V *cs, V *ct,
122  W *currents,
123  int gmode, int ncx, int ncy,
124  U *res);
125 
126  void parallelProp_EH(V *cs, V *ct,
127  W *currents,
128  int gmode, int ncx, int ncy,
129  U *res);
130 
131  void parallelProp_JMEH(V *cs, V *ct,
132  W *currents,
133  int gmode, int ncx, int ncy,
134  U *res);
135 
136  void parallelProp_EHP(V *cs, V *ct,
137  W *currents,
138  int gmode, int ncx, int ncy,
139  U *res);
140 
141  // Functions for calculating angular far-field from reflector directly - no phase term
142  void propagateToFarField(int start, int stop,
143  V *cs, V *ct,
144  W *currents,
145  int gmode, int ncx, int ncy,
146  U *res);
147 
148  std::array<std::array<std::complex<T>, 3>, 2> farfieldAtPoint(V *cs, W *currents,
149  int gmode, int ncx, int ncy,
150  const std::array<T, 3> &point_target);
151 
152  std::array<std::array<std::complex<T>, 3>, 2> farfieldAtPoint_Old(V *cs, W *currents,
153  const std::array<T, 3> &point_target);
154 
155  void parallelFarField(V *cs, V *ct,
156  W *currents,
157  int gmode, int ncx, int ncy,
158  U *res);
159 
160  // Scalar propagation
161  void propagateScalarBeam(int start, int stop,
162  V *cs, V *ct,
163  W *field,
164  int gmode, int ncx, int ncy,
165  U *res);
166 
167 
168  std::complex<T> fieldScalarAtPoint(V *cs, W *field,
169  int gmode, int ncx, int ncy,
170  const std::array<T, 3> &point_target);
171 
172  void parallelPropScalar(V *cs, V *ct,
173  W *field,
174  int gmode, int ncx, int ncy,
175  U *res);
176 
177 };
178 
179 /**
180  * Constructor.
181  *
182  * Set important parameters internally, given input.
183  *
184  * @param omega Angular frequency of radiation, double/float.
185  * @param numThreads Number of computing threads to employ.
186  * @param gs Number of cells on source surface.
187  * @param gt Number of cells on target grid.
188  * @param epsilon Relative electric permittivity of source surface.
189  * @param t_direction Time direction. This changes the sign in the Green's function used to propagate the field.
190  * t_direction is -1 for forward propagation, +1 for backward
191  * @param verbose Whether to print internal state info.
192  */
193 template <class T, class U, class V, class W>
194 Propagation<T, U, V, W>::Propagation(T k, int numThreads, int gs, int gt, T epsilon, T t_direction, bool verbose)
195 {
196  this->Three = 3.0f;
197  this->PIf = 3.14159265359f;
198  this->C_L = 2.99792458e8f; // m s^-1
199  this->MU_0 = 1.2566370614e-6f; // kg m s^-2 A^-2
200  this->EPS_VAC = 1 / (MU_0 * C_L*C_L);
201  this->EPS = epsilon * EPS_VAC; // epsilon is relative permeability
202  this->ZETA = sqrt(MU_0 / EPS);
203  this->ZETA_INV = 1 / ZETA;
204 
205  //printf("EPS %.16g\n", ZETA); // %s is format specifier
206  //printf("MU_0 %.16g\n", ZETA); // %s is format specifier
207  //printf("ZETA %.16g\n", ZETA); // %s is format specifier
208 
209  std::complex<T> j(0., 1.);
210  std::complex<T> z0(0., 0.);
211  this->j = j;
212  this->z0 = z0;
213 
214  this->k = k;
215  //printf("k %.16g\n", k); // %s is format specifier
216 
217  this->prefactor = k*k / (4.*PIf);
218  //printf("prefactor %.16g\n", prefactor); // %s is format specifier
219 
220  this->numThreads = numThreads;
221  this->gs = gs;
222  this->gt = gt;
223 
224  this->step = ceil(gt / numThreads);
225 
226  threadPool.resize(numThreads);
227 
228  this->t_direction = t_direction;
229 
230  this->eye[0].fill(0);
231  this->eye[1].fill(0);
232  this->eye[2].fill(0);
233 
234  this->eye[0][0] = 1;
235  this->eye[1][1] = 1;
236  this->eye[2][2] = 1;
237 
238  if (verbose)
239  {
240  printf("***--------- PO info ---------***\n");
241  printf("--- Source : %d cells\n", gs);
242  printf("--- Target : %d cells\n", gt);
243  printf("--- Wavenumber : %.3f / mm\n", k);
244  printf("--- Threads : %d\n", numThreads);
245  printf("--- Device : CPU\n");
246  printf("***--------- PO info ---------***\n");
247  printf("\n");
248  }
249 }
250 
251 
252 /**
253  * Calculate JM on target - PEC version.
254  *
255  * Calculate the J, M currents on a perfectly conducting target surface.
256  *
257  * @param start Index of first loop iteration in parallel block.
258  * @param stop Index of last loop iteration in parallel block.
259  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
260  * @param ct Pointer to reflcontainer or reflcontainerf object containing target grids.
261  * @param currents Pointer to c2Bundle or c2Bundlef object containing currents on source.
262  * @param res Pointer to c2Bundle or c2Bundlef object, to be filled with calculation results.
263  *
264  * @see reflcontainer
265  * @see reflcontainerf
266  * @see c2Bundle
267  * @see c2Bundlef
268  */
269 template <class T, class U, class V, class W>
270 void Propagation<T,U, V, W>::propagateBeam_JM_PEC(int start, int stop,
271  V *cs, V *ct,
272  W *currents,
273  int gmode, int ncx, int ncy,
274  U *res)
275 {
276  // Arrays of Ts
277  std::array<T, 3> point; // Point on target
278  std::array<T, 3> norms; // Normal vector at point
279 
280  // Arrays of complex Ts to hold return values
281  std::array<std::complex<T>, 3> jt;
282  std::array<std::complex<T>, 3> mt;
283 
284  // Return containers
285  std::array<std::array<std::complex<T>, 3>, 2> beam_e_h; // Container for storing fieldAtPoint return
286 
287  int jc = 0; // Counter
288 
289  for(int i=start; i<stop; i++)
290  {
291 
292  point[0] = ct->x[i];
293  point[1] = ct->y[i];
294  point[2] = ct->z[i];
295 
296  norms[0] = ct->nx[i];
297  norms[1] = ct->ny[i];
298  norms[2] = ct->nz[i];
299 
300  // Calculate total incoming E and H field at point on target
301  beam_e_h = fieldAtPoint(cs, currents,
302  gmode, ncx, ncy,
303  point);
304 
305  ut.ext(norms, beam_e_h[1], jt);
306 
307  res->r1x[i] = 2*jt[0].real();
308  res->r1y[i] = 2*jt[1].real();
309  res->r1z[i] = 2*jt[2].real();
310 
311  res->i1x[i] = 2*jt[0].imag();
312  res->i1y[i] = 2*jt[1].imag();
313  res->i1z[i] = 2*jt[2].imag();
314 
315  res->r2x[i] = 0;
316  res->r2y[i] = 0;
317  res->r2z[i] = 0;
318 
319  res->i2x[i] = 0;
320  res->i2y[i] = 0;
321  res->i2z[i] = 0;
322  }
323 }
324 
325 /**
326  * Calculate JM on target.
327  *
328  * Calculate the J, M currents on a target surface.
329  *
330  * @param start Index of first loop iteration in parallel block.
331  * @param stop Index of last loop iteration in parallel block.
332  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
333  * @param ct Pointer to reflcontainer or reflcontainerf object containing target grids.
334  * @param currents Pointer to c2Bundle or c2Bundlef object containing currents on source.
335  * @param res Pointer to c2Bundle or c2Bundlef object, to be filled with calculation results.
336  *
337  * @see reflcontainer
338  * @see reflcontainerf
339  * @see c2Bundle
340  * @see c2Bundlef
341  */
342 template <class T, class U, class V, class W>
343 void Propagation<T,U, V, W>::propagateBeam_JM(int start, int stop,
344  V *cs, V *ct,
345  W *currents,
346  int gmode, int ncx, int ncy,
347  U *res)
348 {
349  // Scalars (T & complex T)
350  std::complex<T> e_dot_p_r_perp; // E-field - perpendicular reflected POI polarization vector dot product
351  std::complex<T> e_dot_p_r_parr; // E-field - parallel reflected POI polarization vector dot product
352 
353  // Arrays of Ts
354  std::array<T, 3> S_i_norm; // Normalized incoming Poynting vector
355  std::array<T, 3> p_i_perp; // Perpendicular incoming POI polarization vector
356  std::array<T, 3> p_i_parr; // Parallel incoming POI polarization vector
357  std::array<T, 3> S_r_norm; // Normalized reflected Poynting vector
358  std::array<T, 3> p_r_perp; // Perpendicular reflected POI polarization vector
359  std::array<T, 3> p_r_parr; // Parallel reflected POI polarization vector
360  std::array<T, 3> S_out_n; // Container for Poynting-normal ext products
361  std::array<T, 3> point; // Point on target
362  std::array<T, 3> norms; // Normal vector at point
363  std::array<T, 3> e_out_h_r; // Real part of E-field - H-field ext product
364 
365  // Arrays of complex Ts
366  std::array<std::complex<T>, 3> e_r; // Reflected E-field
367  std::array<std::complex<T>, 3> h_r; // Reflected H-field
368  std::array<std::complex<T>, 3> n_out_e_i_r; // Electric current
369  std::array<std::complex<T>, 3> temp1; // Temporary container 1 for intermediate irrelevant values
370  std::array<std::complex<T>, 3> temp2; // Temporary container 2
371 
372  std::array<std::complex<T>, 3> jt;
373  std::array<std::complex<T>, 3> mt;
374 
375  // Return containers
376  std::array<std::array<std::complex<T>, 3>, 2> beam_e_h; // Container for storing fieldAtPoint return
377 
378  int jc = 0; // Counter
379 
380  for(int i=start; i<stop; i++)
381  {
382 
383  point[0] = ct->x[i];
384  point[1] = ct->y[i];
385  point[2] = ct->z[i];
386 
387  norms[0] = ct->nx[i];
388  norms[1] = ct->ny[i];
389  norms[2] = ct->nz[i];
390 
391  // Calculate total incoming E and H field at point on target
392  beam_e_h = fieldAtPoint(cs, currents,
393  gmode, ncx, ncy,
394  point);
395 
396  // Calculate normalized incoming poynting vector.
397  ut.conj(beam_e_h[1], temp1); // h_conj
398  ut.ext(beam_e_h[0], temp1, temp2); // e_out_h
399 
400  for (int n=0; n<3; n++)
401  {
402  e_out_h_r[n] = temp2[n].real(); // e_out_h_r
403  }
404 
405  //std::cout << this->Et_container.size() << std::endl;
406 
407  ut.normalize(e_out_h_r, S_i_norm); // S_i_norm
408 
409  // Calculate incoming polarization vectors
410  ut.ext(S_i_norm, norms, S_out_n); // S_i_out_n
411  ut.normalize(S_out_n, p_i_perp); // p_i_perp
412  ut.ext(p_i_perp, S_i_norm, p_i_parr); // p_i_parr
413 
414  // Now calculate reflected poynting vector.
415  ut.snell(S_i_norm, norms, S_r_norm); // S_r_norm
416 
417  // Calculate normalised reflected polarization vectors
418  ut.ext(S_r_norm, norms, S_out_n); // S_r_out_n
419  ut.normalize(S_out_n, p_r_perp); // p_r_perp
420  ut.ext(S_r_norm, p_r_perp, p_r_parr); // p_r_parr
421 
422  // Now, calculate reflected field from target
423  ut.dot(beam_e_h[0], p_r_perp, e_dot_p_r_perp); // e_dot_p_r_perp
424  ut.dot(beam_e_h[0], p_r_parr, e_dot_p_r_parr); // e_dot_p_r_parr
425 
426  // Calculate reflected field from reflection matrix
427  for(int n=0; n<3; n++)
428  {
429  e_r[n] = -e_dot_p_r_perp * p_i_perp[n] - e_dot_p_r_parr * p_i_parr[n];
430 
431  //this->Et_container[k][i] = e_r[k];
432  }
433 
434  ut.ext(S_r_norm, e_r, temp1); // h_r_temp
435  ut.s_mult(temp1, ZETA_INV, h_r); // h_r
436 
437  for(int n=0; n<3; n++)
438  {
439  temp1[n] = e_r[n] + beam_e_h[0][n]; // e_i_r
440  temp2[n] = h_r[n] + beam_e_h[1][n]; // h_i_r
441  }
442 
443  ut.ext(norms, temp2, jt);
444  ut.ext(norms, temp1, n_out_e_i_r);
445  ut.s_mult(n_out_e_i_r, -1., mt);
446 
447  res->r1x[i] = jt[0].real();
448  res->r1y[i] = jt[1].real();
449  res->r1z[i] = jt[2].real();
450 
451  res->i1x[i] = jt[0].imag();
452  res->i1y[i] = jt[1].imag();
453  res->i1z[i] = jt[2].imag();
454 
455  res->r2x[i] = mt[0].real();
456  res->r2y[i] = mt[1].real();
457  res->r2z[i] = mt[2].real();
458 
459  res->i2x[i] = mt[0].imag();
460  res->i2y[i] = mt[1].imag();
461  res->i2z[i] = mt[2].imag();
462  }
463 }
464 
465 /**
466  * Calculate JMEH on target - PEC version.
467  *
468  * Calculate the J, M currents and E, H fields on a perfectly conducting target surface.
469  *
470  * @param start Index of first loop iteration in parallel block.
471  * @param stop Index of last loop iteration in parallel block.
472  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
473  * @param ct Pointer to reflcontainer or reflcontainerf object containing target grids.
474  * @param currents Pointer to c2Bundle or c2Bundlef object containing currents on source.
475  * @param res Pointer to c2Bundle or c2Bundlef object, to be filled with calculation results.
476  *
477  * @see reflcontainer
478  * @see reflcontainerf
479  * @see c2Bundle
480  * @see c2Bundlef
481  */
482 template <class T, class U, class V, class W>
483 void Propagation<T,U, V, W>::propagateBeam_JMEH_PEC(int start, int stop,
484  V *cs, V *ct,
485  W *currents,
486  int gmode, int ncx, int ncy,
487  U *res)
488 {
489  // Arrays of Ts
490  std::array<T, 3> point; // Point on target
491  std::array<T, 3> norms; // Normal vector at point
492 
493  // Arrays of complex Ts to hold return values
494  std::array<std::complex<T>, 3> jt;
495  std::array<std::complex<T>, 3> mt;
496 
497  // Return containers
498  std::array<std::array<std::complex<T>, 3>, 2> beam_e_h; // Container for storing fieldAtPoint return
499 
500  int jc = 0; // Counter
501 
502  for(int i=start; i<stop; i++)
503  {
504 
505  point[0] = ct->x[i];
506  point[1] = ct->y[i];
507  point[2] = ct->z[i];
508 
509  norms[0] = ct->nx[i];
510  norms[1] = ct->ny[i];
511  norms[2] = ct->nz[i];
512 
513  // Calculate total incoming E and H field at point on target
514  beam_e_h = fieldAtPoint(cs, currents,
515  gmode, ncx, ncy,
516  point);
517 
518  res->r3x[i] = beam_e_h[0][0].real();
519  res->r3y[i] = beam_e_h[0][1].real();
520  res->r3z[i] = beam_e_h[0][2].real();
521 
522  res->i3x[i] = beam_e_h[0][0].imag();
523  res->i3y[i] = beam_e_h[0][1].imag();
524  res->i3z[i] = beam_e_h[0][2].imag();
525 
526  res->r4x[i] = beam_e_h[1][0].real();
527  res->r4y[i] = beam_e_h[1][1].real();
528  res->r4z[i] = beam_e_h[1][2].real();
529 
530  res->i4x[i] = beam_e_h[1][0].imag();
531  res->i4y[i] = beam_e_h[1][1].imag();
532  res->i4z[i] = beam_e_h[1][2].imag();
533 
534  // J_e = 2 n ^ H_inc
535  ut.ext(norms, beam_e_h[1], jt);
536 
537  res->r1x[i] = 2*jt[0].real();
538  res->r1y[i] = 2*jt[1].real();
539  res->r1z[i] = 2*jt[2].real();
540 
541  res->i1x[i] = 2*jt[0].imag();
542  res->i1y[i] = 2*jt[1].imag();
543  res->i1z[i] = 2*jt[2].imag();
544 
545  // J_m = 0
546 
547  res->r2x[i] = 0;
548  res->r2y[i] = 0;
549  res->r2z[i] = 0;
550 
551  res->i2x[i] = 0;
552  res->i2y[i] = 0;
553  res->i2z[i] = 0;
554  }
555 }
556 
557 /**
558  * Calculate JMEH on target.
559  *
560  * Calculate the E, H fields and J, M currents on a target surface.
561  *
562  * @param start Index of first loop iteration in parallel block.
563  * @param stop Index of last loop iteration in parallel block.
564  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
565  * @param ct Pointer to reflcontainer or reflcontainerf object containing target grids.
566  * @param currents Pointer to c2Bundle or c2Bundlef object containing currents on source.
567  * @param res Pointer to c2Bundle or c2Bundlef object, to be filled with calculation results.
568  *
569  * @see reflcontainer
570  * @see reflcontainerf
571  * @see c2Bundle
572  * @see c2Bundlef
573  */
574 template <class T, class U, class V, class W>
575 void Propagation<T, U, V, W>::propagateBeam_EH(int start, int stop,
576  V *cs, V *ct,
577  W *currents,
578  int gmode, int ncx, int ncy,
579  U *res)
580 {
581  // Arrays of Ts
582  std::array<T, 3> point; // Point on target
583 
584  // Return containers
585  std::array<std::array<std::complex<T>, 3>, 2> beam_e_h; // Container for storing fieldAtPoint return
586 
587  int jc = 0; // Counter
588 
589  for(int i=start; i<stop; i++)
590  {
591 
592  point[0] = ct->x[i];
593  point[1] = ct->y[i];
594  point[2] = ct->z[i];
595 
596  // Calculate total incoming E and H field at point on target
597  beam_e_h = fieldAtPoint(cs, currents,
598  gmode, ncx, ncy,
599  point);
600 
601  res->r1x[i] = beam_e_h[0][0].real();
602  res->r1y[i] = beam_e_h[0][1].real();
603  res->r1z[i] = beam_e_h[0][2].real();
604 
605  res->i1x[i] = beam_e_h[0][0].imag();
606  res->i1y[i] = beam_e_h[0][1].imag();
607  res->i1z[i] = beam_e_h[0][2].imag();
608 
609  res->r2x[i] = beam_e_h[1][0].real();
610  res->r2y[i] = beam_e_h[1][1].real();
611  res->r2z[i] = beam_e_h[1][2].real();
612 
613  res->i2x[i] = beam_e_h[1][0].imag();
614  res->i2y[i] = beam_e_h[1][1].imag();
615  res->i2z[i] = beam_e_h[1][2].imag();
616  }
617 }
618 
619 /**
620  * Calculate JM and EH on target.
621  *
622  * Calculate the J, M currents and E, H fields on a target surface.
623  *
624  * @param start Index of first loop iteration in parallel block.
625  * @param stop Index of last loop iteration in parallel block.
626  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
627  * @param ct Pointer to reflcontainer or reflcontainerf object containing target grids.
628  * @param currents Pointer to c2Bundle or c2Bundlef object containing currents on source.
629  * @param res Pointer to c4Bundle or c4Bundlef object, to be filled with calculation results.
630  *
631  * @see reflcontainer
632  * @see reflcontainerf
633  * @see c2Bundle
634  * @see c2Bundlef
635  * @see c4Bundle
636  * @see c4Bundlef
637  */
638 template <class T, class U, class V, class W>
639 void Propagation<T, U, V, W>::propagateBeam_JMEH(int start, int stop,
640  V *cs, V *ct,
641  W *currents,
642  int gmode, int ncx, int ncy,
643  U *res)
644 {
645  // Scalars (T & complex T)
646  std::complex<T> e_dot_p_r_perp; // E-field - perpendicular reflected POI polarization vector dot product
647  std::complex<T> e_dot_p_r_parr; // E-field - parallel reflected POI polarization vector dot product
648 
649  // Arrays of Ts
650  std::array<T, 3> S_i_norm; // Normalized incoming Poynting vector
651  std::array<T, 3> p_i_perp; // Perpendicular incoming POI polarization vector
652  std::array<T, 3> p_i_parr; // Parallel incoming POI polarization vector
653  std::array<T, 3> S_r_norm; // Normalized reflected Poynting vector
654  std::array<T, 3> p_r_perp; // Perpendicular reflected POI polarization vector
655  std::array<T, 3> p_r_parr; // Parallel reflected POI polarization vector
656  std::array<T, 3> S_out_n; // Container for Poynting-normal ext products
657  std::array<T, 3> point; // Point on target
658  std::array<T, 3> norms; // Normal vector at point
659  std::array<T, 3> e_out_h_r; // Real part of E-field - H-field ext product
660 
661  // Arrays of complex Ts
662  std::array<std::complex<T>, 3> e_r; // Reflected E-field
663  std::array<std::complex<T>, 3> h_r; // Reflected H-field
664  std::array<std::complex<T>, 3> n_out_e_i_r; // Electric current
665  std::array<std::complex<T>, 3> temp1; // Temporary container 1 for intermediate irrelevant values
666  std::array<std::complex<T>, 3> temp2; // Temporary container 2
667 
668  std::array<std::complex<T>, 3> jt;
669  std::array<std::complex<T>, 3> mt;
670 
671  // Return containers
672  std::array<std::array<std::complex<T>, 3>, 2> beam_e_h; // Container for storing fieldAtPoint return
673 
674  int jc = 0; // Counter
675 
676  for(int i=start; i<stop; i++)
677  {
678 
679  point[0] = ct->x[i];
680  point[1] = ct->y[i];
681  point[2] = ct->z[i];
682 
683  norms[0] = ct->nx[i];
684  norms[1] = ct->ny[i];
685  norms[2] = ct->nz[i];
686 
687  // Calculate total incoming E and H field at point on target
688  beam_e_h = fieldAtPoint(cs, currents,
689  gmode, ncx, ncy,
690  point);
691 
692  //res->r1x[i] = 0.;//beam_e_h[0][0].real();
693 
694  res->r3x[i] = beam_e_h[0][0].real();
695  res->r3y[i] = beam_e_h[0][1].real();
696  res->r3z[i] = beam_e_h[0][2].real();
697 
698  res->i3x[i] = beam_e_h[0][0].imag();
699  res->i3y[i] = beam_e_h[0][1].imag();
700  res->i3z[i] = beam_e_h[0][2].imag();
701 
702  res->r4x[i] = beam_e_h[1][0].real();
703  res->r4y[i] = beam_e_h[1][1].real();
704  res->r4z[i] = beam_e_h[1][2].real();
705 
706  res->i4x[i] = beam_e_h[1][0].imag();
707  res->i4y[i] = beam_e_h[1][1].imag();
708  res->i4z[i] = beam_e_h[1][2].imag();
709 
710  // Calculate normalised incoming poynting vector.
711  ut.conj(beam_e_h[1], temp1); // h_conj
712 
713  ut.ext(beam_e_h[0], temp1, temp2); // e_out_h
714 
715  for (int n=0; n<3; n++)
716  {
717  e_out_h_r[n] = temp2[n].real(); // e_out_h_r
718  }
719 
720  //std::cout << this->Et_container.size() << std::endl;
721 
722  ut.normalize(e_out_h_r, S_i_norm); // S_i_norm
723 
724  // Calculate incoming polarization vectors
725  ut.ext(S_i_norm, norms, S_out_n); // S_i_out_n
726 
727  ut.normalize(S_out_n, p_i_perp); // p_i_perp
728  ut.ext(p_i_perp, S_i_norm, p_i_parr); // p_i_parr
729 
730  // Now calculate reflected poynting vector.
731  ut.snell(S_i_norm, norms, S_r_norm); // S_r_norm
732 
733  // Calculate normalised reflected polarization vectors
734  ut.ext(S_r_norm, norms, S_out_n); // S_r_out_n
735 
736  ut.normalize(S_out_n, p_r_perp); // p_r_perp
737  ut.ext(S_r_norm, p_r_perp, p_r_parr); // p_r_parr
738 
739  // Now, calculate reflected field from target
740  ut.dot(beam_e_h[0], p_r_perp, e_dot_p_r_perp); // e_dot_p_r_perp
741  ut.dot(beam_e_h[0], p_r_parr, e_dot_p_r_parr); // e_dot_p_r_parr
742 
743  //res->r1x[i] = beam_e_h[0][0].real();
744 
745  // Calculate reflected field from reflection matrix
746  for(int n=0; n<3; n++)
747  {
748  e_r[n] = -e_dot_p_r_perp * p_i_perp[n] - e_dot_p_r_parr * p_i_parr[n];
749 
750  //this->Et_container[k][i] = e_r[k];
751  }
752 
753  ut.ext(S_r_norm, e_r, temp1); // h_r_temp
754  ut.s_mult(temp1, ZETA_INV, h_r); // h_r
755 
756  for(int n=0; n<3; n++)
757  {
758  temp1[n] = e_r[n] + beam_e_h[0][n]; // e_i_r
759  temp2[n] = h_r[n] + beam_e_h[1][n]; // h_i_r
760  }
761 
762  ut.ext(norms, temp2, jt);
763  ut.ext(norms, temp1, n_out_e_i_r);
764  ut.s_mult(n_out_e_i_r, -1., mt);
765 
766 
767 
768  res->r1x[i] = jt[0].real();
769  res->r1y[i] = jt[1].real();
770  res->r1z[i] = jt[2].real();
771 
772  res->i1x[i] = jt[0].imag();
773  res->i1y[i] = jt[1].imag();
774  res->i1z[i] = jt[2].imag();
775 
776  res->r2x[i] = mt[0].real();
777  res->r2y[i] = mt[1].real();
778  res->r2z[i] = mt[2].real();
779 
780  res->i2x[i] = mt[0].imag();
781  res->i2y[i] = mt[1].imag();
782  res->i2z[i] = mt[2].imag();
783  }
784 }
785 
786 /**
787  * Calculate reflected EH and P on target.
788  *
789  * Calculate the reflected E, H fields and P, the reflected Poynting vector field, on a target surface.
790  *
791  * @param start Index of first loop iteration in parallel block.
792  * @param stop Index of last loop iteration in parallel block.
793  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
794  * @param ct Pointer to reflcontainer or reflcontainerf object containing target grids.
795  * @param currents Pointer to c2Bundle or c2Bundlef object containing currents on source.
796  * @param res Pointer to c2rBundle or c2rBundlef object, to be filled with calculation results.
797  *
798  * @see reflcontainer
799  * @see reflcontainerf
800  * @see c2Bundle
801  * @see c2Bundlef
802  * @see c2rBundle
803  * @see c2rBundlef
804  */
805 template <class T, class U, class V, class W>
806 void Propagation<T, U, V, W>::propagateBeam_EHP(int start, int stop,
807  V *cs, V *ct,
808  W *currents,
809  int gmode, int ncx, int ncy,
810  U *res)
811 {
812  // Scalars (T & complex T)
813  std::complex<T> e_dot_p_r_perp; // E-field - perpendicular reflected POI polarization vector dot product
814  std::complex<T> e_dot_p_r_parr; // E-field - parallel reflected POI polarization vector dot product
815 
816  // Arrays of Ts
817  std::array<T, 3> S_i_norm; // Normalized incoming Poynting vector
818  std::array<T, 3> p_i_perp; // Perpendicular incoming POI polarization vector
819  std::array<T, 3> p_i_parr; // Parallel incoming POI polarization vector
820  std::array<T, 3> S_r_norm; // Normalized reflected Poynting vector
821  std::array<T, 3> p_r_perp; // Perpendicular reflected POI polarization vector
822  std::array<T, 3> p_r_parr; // Parallel reflected POI polarization vector
823  std::array<T, 3> S_out_n; // Container for Poynting-normal ext products
824  std::array<T, 3> point; // Point on target
825  std::array<T, 3> norms; // Normal vector at point
826  std::array<T, 3> e_out_h_r; // Real part of E-field - H-field ext product
827 
828  // Arrays of complex Ts
829  std::array<std::complex<T>, 3> e_r; // Reflected E-field
830  std::array<std::complex<T>, 3> h_r; // Reflected H-field
831  std::array<std::complex<T>, 3> n_out_e_i_r; // Electric current
832  std::array<std::complex<T>, 3> temp1; // Temporary container 1 for intermediate irrelevant values
833  std::array<std::complex<T>, 3> temp2; // Temporary container 2
834 
835  // Return containers
836  std::array<std::array<std::complex<T>, 3>, 2> beam_e_h; // Container for storing fieldAtPoint return
837 
838  int jc = 0; // Counter
839 
840  for(int i=start; i<stop; i++)
841  {
842 
843  point[0] = ct->x[i];
844  point[1] = ct->y[i];
845  point[2] = ct->z[i];
846 
847  norms[0] = ct->nx[i];
848  norms[1] = ct->ny[i];
849  norms[2] = ct->nz[i];
850 
851  // Calculate total incoming E and H field at point on target
852  beam_e_h = fieldAtPoint(cs, currents, gmode, ncx, ncy,
853  point);
854 
855  // Calculate normalised incoming poynting vector.
856  ut.conj(beam_e_h[1], temp1); // h_conj
857  ut.ext(beam_e_h[0], temp1, temp2); // e_out_h
858 
859  for (int n=0; n<3; n++)
860  {
861  e_out_h_r[n] = temp2[n].real(); // e_out_h_r
862  }
863 
864  //std::cout << this->Et_container.size() << std::endl;
865 
866  ut.normalize(e_out_h_r, S_i_norm); // S_i_norm
867 
868  // Calculate incoming polarization vectors
869  ut.ext(S_i_norm, norms, S_out_n); // S_i_out_n
870  ut.normalize(S_out_n, p_i_perp); // p_i_perp
871  ut.ext(p_i_perp, S_i_norm, p_i_parr); // p_i_parr
872 
873  // Now calculate reflected poynting vector.
874  ut.snell(S_i_norm, norms, S_r_norm); // S_r_norm
875 
876  // Calculate normalised reflected polarization vectors
877  ut.ext(S_r_norm, norms, S_out_n); // S_r_out_n
878  ut.normalize(S_out_n, p_r_perp); // p_r_perp
879  ut.ext(S_r_norm, p_r_perp, p_r_parr); // p_r_parr
880 
881  // Now, calculate reflected field from target
882  ut.dot(beam_e_h[0], p_r_perp, e_dot_p_r_perp); // e_dot_p_r_perp
883  ut.dot(beam_e_h[0], p_r_parr, e_dot_p_r_parr); // e_dot_p_r_parr
884 
885  // Calculate reflected field from reflection matrix
886  for(int n=0; n<3; n++)
887  {
888  e_r[n] = -e_dot_p_r_perp * p_i_perp[n] - e_dot_p_r_parr * p_i_parr[n];
889  }
890 
891  ut.ext(S_r_norm, e_r, temp1); // h_r_temp
892  ut.s_mult(temp1, ZETA_INV, h_r); // h_r
893 
894  res->r1x[i] = e_r[0].real();
895  res->r1y[i] = e_r[1].real();
896  res->r1z[i] = e_r[2].real();
897 
898  res->i1x[i] = e_r[0].imag();
899  res->i1y[i] = e_r[1].imag();
900  res->i1z[i] = e_r[2].imag();
901 
902  res->r2x[i] = h_r[0].real();
903  res->r2y[i] = h_r[1].real();
904  res->r2z[i] = h_r[2].real();
905 
906  res->i2x[i] = h_r[0].imag();
907  res->i2y[i] = h_r[1].imag();
908  res->i2z[i] = h_r[2].imag();
909 
910  res->r3x[i] = S_r_norm[0];
911  res->r3y[i] = S_r_norm[1];
912  res->r3z[i] = S_r_norm[2];
913  }
914 }
915 
916 /**
917  * Calculate scalar field on target.
918  *
919  * Calculate the scalar fields on a target surface.
920  *
921  * @param start Index of first loop iteration in parallel block.
922  * @param stop Index of last loop iteration in parallel block.
923  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
924  * @param ct Pointer to reflcontainer or reflcontainerf object containing target grids.
925  * @param field Pointer to arrC1 or arrC1f object containing field on source.
926  * @param res Pointer to arrC1 or arrC1f object, to be filled with calculation results.
927  *
928  * @see reflcontainer
929  * @see reflcontainerf
930  * @see c2Bundle
931  * @see c2Bundlef
932  */
933 template <class T, class U, class V, class W>
934 void Propagation<T, U, V, W>::propagateScalarBeam(int start, int stop,
935  V *cs, V *ct,
936  W *field,
937  int gmode, int ncx, int ncy,
938  U *res)
939 {
940  std::array<T, 3> point_target;
941  std::complex<T> ets;
942 
943  for(int i=start; i<stop; i++)
944  {
945  point_target[0] = ct->x[i];
946  point_target[1] = ct->y[i];
947  point_target[2] = ct->z[i];
948 
949  ets = fieldScalarAtPoint(cs, field,
950  gmode, ncx, ncy,
951  point_target);
952 
953  res->x[i] = ets.real();
954  res->y[i] = ets.imag();
955  }
956 }
957 
958 /**
959  * Calculate field on target.
960  *
961  * Calculate integrated E and H fields on a target point.
962  *
963  * We use the GRASP formulae for the E and H fields in sqrt{Watts}, calculated
964  * from electric and magnetic currents in sqrt{Watts}.
965  *
966  * The GRASP manual cites Collin "Antenna Theory" (1969) for the derivation.
967  *
968  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
969  * @param currents Pointer to c2Bundle or c2Bundlef object containing currents on source.
970  * @param point_target Array of 3 double/float containing target point co-ordinate.
971  *
972  * @see reflcontainer
973  * @see reflcontainerf
974  * @see c2Bundle
975  * @see c2Bundlef
976  */
977 template <class T, class U, class V, class W>
978 std::array<std::array<std::complex<T>, 3>, 2> Propagation<T, U, V, W>::fieldAtPoint(V *cs,
979  W *currents,
980  int gmode, int ncx, int ncy,
981  const std::array<T, 3> &point_target)
982 {
983  int i;
984  // Scalars (T & complex T)
985  T R; // Distance between source and target points
986  T R_inv; // 1 / R
987  T kR; // The product of the wavevector and the the distance
988  T kR_inv; // The inverse of the product of the wavevector and the the distance
989  T norm_dot_R_hat; // Source normal dotted with wavevector direction
990  T half = 0.5;
991 
992  std::complex<T> Green; // Container for Green's function
993  std::complex<T> js_dot_R; // Magnitude of electric currents on to R_hat
994  std::complex<T> ms_dot_R; // Magnitude of magnetic currents on to R_hat
995  std::complex<T> kR_inv_sum1; // Container for the first complex sum of 1/kR powers
996  std::complex<T> kR_inv_sum2; // Container for the second complex sum of 1/kR powers
997  std::complex<T> kR_inv_sum3; // Container for the third complex sum of 1/kR powers
998 
999 
1000  // Arrays of Ts
1001  std::array<T, 3> source_point; // Container for xyz co-ordinates
1002  std::array<T, 3> source_norm; // Container for xyz normals.
1003  std::array<T, 3> R_vec; // Distance vector between source and target points
1004  std::array<T, 3> R_hat; // Unit vector between source and target points
1005  std::array<T, 3> k_arr; // Wavevector
1006 
1007  // Arrays of complex Ts
1008  std::array<std::complex<T>, 3> e_field; // Electric field surface integral at point
1009  std::array<std::complex<T>, 3> h_field; // Magnetic field surface integral at point
1010  std::array<std::complex<T>, 3> ye_field; // Electric field surface integral at point
1011  std::array<std::complex<T>, 3> yh_field; // Magnetic field surface integral at point
1012  std::array<std::complex<T>, 3> js; // Electric current at source point
1013  std::array<std::complex<T>, 3> ms; // Magnetic current at source point
1014  std::array<std::complex<T>, 3> k_out_ms; // Outer product between k and ms
1015  std::array<std::complex<T>, 3> k_out_js; // Outer product between k and js
1016  std::array<std::complex<T>, 3> js_dot_R_R; // Electric currents projected on to R_hat
1017  std::array<std::complex<T>, 3> ms_cross_R; // Cross product of magnetic currents and R_hat
1018  std::array<std::complex<T>, 3> ms_dot_R_R; // Magnetic currents projected on to R_hat
1019  std::array<std::complex<T>, 3> js_cross_R; // Cross product of electric currents and R_hat
1020 
1021 
1022  // Return container
1023  std::array<std::array<std::complex<T>, 3>, 2> e_h_field; // Return container containing e and h-fields
1024 
1025  e_field.fill(z0);
1026  h_field.fill(z0);
1027 
1028  for(int xu=0; xu<ncx; xu++)
1029  {
1030  for (int n=0; n<3; n++)
1031  {
1032  ye_field[n] = z0; // Intermediate electric field due to integral over y/v
1033  yh_field[n] = z0; // Intermediate magnetic field due to integral over y/v
1034  }
1035 
1036  for(int yv=0; yv<ncy; yv++)
1037  {
1038  i = xu*ncy + yv;
1039 
1040  source_point[0] = cs->x[i];
1041  source_point[1] = cs->y[i];
1042  source_point[2] = cs->z[i];
1043 
1044  source_norm[0] = cs->nx[i];
1045  source_norm[1] = cs->ny[i];
1046  source_norm[2] = cs->nz[i];
1047 
1048  js[0] = {currents->r1x[i], currents->i1x[i]};
1049  js[1] = {currents->r1y[i], currents->i1y[i]};
1050  js[2] = {currents->r1z[i], currents->i1z[i]};
1051 
1052  ms[0] = {currents->r2x[i], currents->i2x[i]};
1053  ms[1] = {currents->r2y[i], currents->i2y[i]};
1054  ms[2] = {currents->r2z[i], currents->i2z[i]};
1055 
1056  ut.diff(point_target, source_point, R_vec);
1057  ut.abs(R_vec, R);
1058 
1059  //printf("R %.16g\n", R); // %s is format specifier
1060 
1061  R_inv = 1 / R;
1062 
1063  //printf("1/R %.16g\n", R_inv); // %s is format specifier
1064 
1065  kR = k*R;
1066 
1067  //printf("kR %.16g\n", kR); // %s is format specifier
1068 
1069  kR_inv = 1/kR;
1070 
1071  //printf("1/kR %.16g\n", kR_inv); // %s is format specifier
1072 
1073  ut.s_mult(R_vec, R_inv, R_hat);
1074 
1075  //printf("R_hat: (%.16g, %.16g, %.16g)\n", R_hat[0], R_hat[1], R_hat[2]); // %s is format specifier
1076 
1077  // Check if source point illuminates target point or not.
1078  ut.dot(source_norm, R_hat, norm_dot_R_hat);
1079  //printf("source_norm: (%.16g, %.16g, %.16g)\n", source_norm[0], source_norm[1], source_norm[2]); // %s is format specifier
1080  //printf("n . R_hat: %.16g\n", norm_dot_R_hat); // %s is format specifier
1081  if (norm_dot_R_hat < 0) {continue;}
1082 
1083  // Calculate the complex sums that appear in the integral
1084  kR_inv_sum1 = kR_inv * (-j + t_direction*kR_inv + j*kR_inv*kR_inv);
1085  //printf("kR_inv_sum1: %.16g+%16gi\n", kR_inv_sum1.real(), kR_inv_sum1.imag()); // %s is format specifier
1086 
1087  kR_inv_sum2 = kR_inv * (j - t_direction*Three*kR_inv - Three*j * kR_inv*kR_inv);
1088  //printf("kR_inv_sum2: %.16g+%16gi\n", kR_inv_sum2.real(), kR_inv_sum2.imag()); // %s is format specifier
1089 
1090  kR_inv_sum3 = t_direction * (j*kR_inv - t_direction*kR_inv*kR_inv);
1091  //printf("kR_inv_sum3: %.16g+%16gi\n", kR_inv_sum3.real(), kR_inv_sum3.imag()); // %s is format specifier
1092 
1093 
1094  // e-field
1095  ut.dot(js, R_hat, js_dot_R);
1096  ut.s_mult(R_hat, js_dot_R, js_dot_R_R);
1097  ut.ext(ms, R_hat, ms_cross_R);
1098 
1099  // h-field
1100  ut.dot(ms, R_hat, ms_dot_R);
1101  ut.s_mult(R_hat, ms_dot_R, ms_dot_R_R);
1102  ut.ext(js, R_hat, js_cross_R);
1103 
1104  Green = prefactor * exp(t_direction * j * k * R) * cs->area[i];
1105  //printf("%.16g, %.16g\n", Green.real(), Green.imag()); // %s is format specifier
1106 
1107  for( int n=0; n<3; n++)
1108  {
1109  if ((gmode != 1) && ((yv==0) || (yv==ncy-1)))
1110  {
1111  ye_field[n] += half*(js[n] * kR_inv_sum1 + js_dot_R_R[n] * kR_inv_sum2 + ms_cross_R[n] * kR_inv_sum3) * Green;
1112  yh_field[n] += half*(ms[n] * kR_inv_sum1 + ms_dot_R_R[n] * kR_inv_sum2 - js_cross_R[n] * kR_inv_sum3) * Green;
1113  }
1114  else
1115  {
1116  ye_field[n] += (js[n] * kR_inv_sum1 + js_dot_R_R[n] * kR_inv_sum2 + ms_cross_R[n] * kR_inv_sum3) * Green;
1117  yh_field[n] += (ms[n] * kR_inv_sum1 + ms_dot_R_R[n] * kR_inv_sum2 - js_cross_R[n] * kR_inv_sum3) * Green;
1118  }
1119  }
1120  } // End of yv loop
1121  for (int n=0; n<3; n++)
1122  {
1123  if ((xu==0) || (xu==ncx-1)) // Only add half the point value
1124  {
1125  e_field[n] += half*ye_field[n];
1126  h_field[n] += half*yh_field[n];
1127  }
1128  else
1129  {
1130  e_field[n] += ye_field[n];
1131  h_field[n] += yh_field[n];
1132  }
1133  }
1134 
1135  } // End of xu loop
1136 
1137  // Pack e and h together in single container
1138  e_h_field[0] = e_field;
1139  e_h_field[1] = h_field;
1140 
1141  //std::cout << ut.abs(e_field) << std::endl;
1142 
1143  return e_h_field;
1144 }
1145 
1146 
1147 /**
1148  * Calculate field on target.
1149  *
1150  * Calculate integrated E and H fields on a target point.
1151  *
1152  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
1153  * @param currents Pointer to c2Bundle or c2Bundlef object containing currents on source.
1154  * @param point_target Array of 3 double/float containing target point co-ordinate.
1155  *
1156  * @see reflcontainer
1157  * @see reflcontainerf
1158  * @see c2Bundle
1159  * @see c2Bundlef
1160  */
1161 template <class T, class U, class V, class W>
1162 std::array<std::array<std::complex<T>, 3>, 2> Propagation<T, U, V, W>::fieldAtPoint_Old(V *cs,
1163  W *currents, const std::array<T, 3> &point_target)
1164 {
1165  // Scalars (T & complex T)
1166  T r; // Distance between source and target points
1167  T r_inv; // 1 / r
1168  T omega; // Angular frequency of field
1169  T norm_dot_k_hat; // Source normal dotted with wavevector direction
1170  std::complex<T> Green; // Container for Green's function
1171  std::complex<T> r_in_s; // Container for inner products between wavevctor and currents
1172 
1173  // Arrays of Ts
1174  std::array<T, 3> source_point; // Container for xyz co-ordinates
1175  std::array<T, 3> source_norm; // Container for xyz normals.
1176  std::array<T, 3> r_vec; // Distance vector between source and target points
1177  std::array<T, 3> k_hat; // Unit wavevctor
1178  std::array<T, 3> k_arr; // Wavevector
1179 
1180  // Arrays of complex Ts
1181  std::array<std::complex<T>, 3> e_field; // Electric field on target
1182  std::array<std::complex<T>, 3> h_field; // Magnetic field on target
1183  std::array<std::complex<T>, 3> js; // Electric current at source point
1184  std::array<std::complex<T>, 3> ms; // Magnetic current at source point
1185  std::array<std::complex<T>, 3> e_vec_thing; // Electric current contribution to e-field
1186  std::array<std::complex<T>, 3> h_vec_thing; // Magnetic current contribution to h-field
1187  std::array<std::complex<T>, 3> k_out_ms; // Outer product between k and ms
1188  std::array<std::complex<T>, 3> k_out_js; // Outer product between k and js
1189  std::array<std::complex<T>, 3> temp; // Temporary container for intermediate values
1190 
1191  // Return container
1192  std::array<std::array<std::complex<T>, 3>, 2> e_h_field; // Return container containing e and h-fields
1193 
1194  e_field.fill(z0);
1195  h_field.fill(z0);
1196 
1197  omega = C_L * k;
1198 
1199  for(int i=0; i<gs; i++)
1200  {
1201  source_point[0] = cs->x[i];
1202  source_point[1] = cs->y[i];
1203  source_point[2] = cs->z[i];
1204 
1205  source_norm[0] = cs->nx[i];
1206  source_norm[1] = cs->ny[i];
1207  source_norm[2] = cs->nz[i];
1208 
1209  js[0] = {currents->r1x[i], currents->i1x[i]};
1210  js[1] = {currents->r1y[i], currents->i1y[i]};
1211  js[2] = {currents->r1z[i], currents->i1z[i]};
1212 
1213  ms[0] = {currents->r2x[i], currents->i2x[i]};
1214  ms[1] = {currents->r2y[i], currents->i2y[i]};
1215  ms[2] = {currents->r2z[i], currents->i2z[i]};
1216 
1217  ut.diff(point_target, source_point, r_vec);
1218  ut.abs(r_vec, r);
1219 
1220  r_inv = 1 / r;
1221 
1222  ut.s_mult(r_vec, r_inv, k_hat);
1223 
1224  // Check if source point illuminates target point or not.
1225  ut.dot(source_norm, k_hat, norm_dot_k_hat);
1226  if (norm_dot_k_hat < 0) {continue;}
1227 
1228  ut.s_mult(k_hat, k, k_arr);
1229 
1230  // e-field
1231  ut.dot(k_hat, js, r_in_s);
1232  ut.s_mult(k_hat, r_in_s, temp);
1233  ut.diff(js, temp, e_vec_thing);
1234 
1235  ut.ext(k_arr, ms, k_out_ms);
1236 
1237  // h-field
1238  ut.dot(k_hat, ms, r_in_s);
1239  ut.s_mult(k_hat, r_in_s, temp);
1240  ut.diff(ms, temp, h_vec_thing);
1241 
1242  ut.ext(k_arr, js, k_out_js);
1243 
1244  //printf("%.16g\n", r);
1245 
1246  Green = exp(this->t_direction * j * k * r) / (4 * PIf * r) * cs->area[i] * j;
1247 
1248  for( int n=0; n<3; n++)
1249  {
1250  e_field[n] += (-omega * MU_0 * e_vec_thing[n] + k_out_ms[n]) * Green;
1251  h_field[n] += (-omega * EPS * h_vec_thing[n] - k_out_js[n]) * Green;
1252  }
1253  //printf("%.16g, %.16g\n", Green.real(), Green.imag()); // %s is format specifier
1254  }
1255 
1256  // Pack e and h together in single container
1257  e_h_field[0] = e_field;
1258  e_h_field[1] = h_field;
1259 
1260  //std::cout << ut.abs(e_field) << std::endl;
1261 
1262  return e_h_field;
1263 }
1264 
1265 /**
1266  * Calculate scalar field on target.
1267  *
1268  * Calculate integrated scalar field on a target point.
1269  *
1270  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
1271  * @param field Pointer to arrC1 or arrC1f object containing currents on source.
1272  * @param point_target Array of 3 double/float containing target point co-ordinate.
1273  *
1274  * @see reflcontainer
1275  * @see reflcontainerf
1276  * @see arrC1
1277  * @see arrC1f
1278  */
1279 template <class T, class U, class V, class W>
1280 std::complex<T> Propagation<T, U, V, W>::fieldScalarAtPoint(V *cs,
1281  W *field,
1282  int gmode, int ncx, int ncy,
1283  const std::array<T, 3> &point_target)
1284 {
1285  int i;
1286  std::complex<T> out(0., 0.);
1287  std::complex<T> j(0., 1.);
1288  std::complex<T> _field;
1289  std::complex<T> _yfield;
1290  T half = 0.5;
1291 
1292  T r;
1293  std::array<T, 3> r_vec;
1294  std::array<T, 3> source_point;
1295 
1296  for(int xu=0; xu<ncx; xu++)
1297  {
1298  _yfield = z0; // Intermediate electric field due to integral over y/v
1299 
1300  for(int yv=0; yv<ncy; yv++)
1301  {
1302  i = xu*ncy + yv;
1303  source_point[0] = cs->x[i];
1304  source_point[1] = cs->y[i];
1305  source_point[2] = cs->z[i];
1306 
1307  _field = {field->x[i], field->y[i]};
1308 
1309  ut.diff(point_target, source_point, r_vec);
1310  ut.abs(r_vec, r);
1311 
1312  if ((gmode != 1) && ((yv==0) || (yv==ncy-1)))
1313  {
1314  _yfield += - half * k * k * _field * exp(this->t_direction * -j * k * r) / (4 * PIf * r) * cs->area[i];
1315  }
1316  else
1317  {
1318  _yfield += - k * k * _field * exp(this->t_direction * -j * k * r) / (4 * PIf * r) * cs->area[i];
1319  }
1320  } // End of yv loop
1321  if ((xu==0) || (xu==ncx-1)) // Only add half the point value
1322  {
1323  out += half*_yfield;
1324  }
1325  else
1326  {
1327  out += _yfield;
1328  }
1329  }
1330  return out;
1331 }
1332 
1333 /**
1334  * Calculate JM parallel.
1335  *
1336  * Run the propagateBeam_JM function in parallel blocks.
1337  *
1338  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
1339  * @param ct Pointer to reflcontainer or reflcontainerf object containing target grids.
1340  * @param currents Pointer to c2Bundle or c2Bundlef object containing source currents.
1341  * @param res Pointer to c2Bundle or c2Bundlef object to be filled with calculation results.
1342  *
1343  * @see reflcontainer
1344  * @see reflcontainerf
1345  * @see c2Bundle
1346  * @see c2Bundlef
1347  */
1348 template <class T, class U, class V, class W>
1349 void Propagation<T, U, V, W>::parallelProp_JM(V *cs, V *ct,
1350  W *currents,
1351  int gmode, int ncx, int ncy,
1352  U *res)
1353 {
1354  int final_step;
1355 
1356  for(int n=0; n<numThreads; n++)
1357  {
1358  //std::cout << n << std::endl;
1359  if(n == (numThreads-1))
1360  {
1361  final_step = gt;
1362  }
1363 
1364  else
1365  {
1366  final_step = (n+1) * step;
1367  }
1368 
1369  //std::cout << final_step << std::endl;
1370 
1371  threadPool[n] = std::thread(&Propagation::propagateBeam_JM_PEC,
1372  this, n * step, final_step,
1373  cs, ct, currents,
1374  gmode, ncx, ncy,
1375  res);
1376  }
1377  joinThreads();
1378 }
1379 
1380 /**
1381  * Calculate EH parallel.
1382  *
1383  * Run the propagateBeam_EH function in parallel blocks.
1384  *
1385  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
1386  * @param ct Pointer to reflcontainer or reflcontainerf object containing target grids.
1387  * @param currents Pointer to c2Bundle or c2Bundlef object containing source currents.
1388  * @param res Pointer to c2Bundle or c2Bundlef object to be filled with calculation results.
1389  *
1390  * @see reflcontainer
1391  * @see reflcontainerf
1392  * @see c2Bundle
1393  * @see c2Bundlef
1394  */
1395 template <class T, class U, class V, class W>
1396 void Propagation<T, U, V, W>::parallelProp_EH(V *cs, V *ct,
1397  W *currents,
1398  int gmode, int ncx, int ncy,
1399  U *res)
1400 {
1401  int final_step;
1402 
1403  for(int n=0; n<numThreads; n++)
1404  {
1405  //std::cout << n << std::endl;
1406  if(n == (numThreads-1))
1407  {
1408  final_step = gt;
1409  }
1410 
1411  else
1412  {
1413  final_step = (n+1) * step;
1414  }
1415 
1416  //std::cout << final_step << std::endl;
1417 
1418  threadPool[n] = std::thread(&Propagation::propagateBeam_EH,
1419  this, n * step, final_step,
1420  cs, ct, currents,
1421  gmode, ncx, ncy,
1422  res);
1423  }
1424  joinThreads();
1425 }
1426 
1427 /**
1428  * Calculate JM and EH parallel.
1429  *
1430  * Run the propagateBeam_JMEH function in parallel blocks.
1431  *
1432  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
1433  * @param ct Pointer to reflcontainer or reflcontainerf object containing target grids.
1434  * @param currents Pointer to c2Bundle or c2Bundlef object containing source currents.
1435  * @param res Pointer to c4Bundle or c4Bundlef object to be filled with calculation results.
1436  *
1437  * @see reflcontainer
1438  * @see reflcontainerf
1439  * @see c2Bundle
1440  * @see c2Bundlef
1441  * @see c4Bundle
1442  * @see c4Bundlef
1443  */
1444 template <class T, class U, class V, class W>
1445 void Propagation<T, U, V, W>::parallelProp_JMEH(V *cs, V *ct,
1446  W *currents,
1447  int gmode, int ncx, int ncy,
1448  U *res)
1449 {
1450  int final_step;
1451 
1452  for(int n=0; n<numThreads; n++)
1453  {
1454  if(n == (numThreads-1))
1455  {
1456  final_step = gt;
1457  }
1458 
1459  else
1460  {
1461  final_step = (n+1) * step;
1462  }
1463 
1464  threadPool[n] = std::thread(&Propagation::propagateBeam_JMEH_PEC,
1465  this, n * step, final_step,
1466  cs, ct, currents,
1467  gmode, ncx, ncy,
1468  res);
1469  }
1470  joinThreads();
1471 }
1472 
1473 /**
1474  * Calculate reflected EH and P parallel.
1475  *
1476  * Run the propagateBeam_EHP function in parallel blocks.
1477  *
1478  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
1479  * @param ct Pointer to reflcontainer or reflcontainerf object containing target grids.
1480  * @param currents Pointer to c2Bundle or c2Bundlef object containing source currents.
1481  * @param res Pointer to c2rBundle or c2rBundlef object to be filled with calculation results.
1482  *
1483  * @see reflcontainer
1484  * @see reflcontainerf
1485  * @see c2Bundle
1486  * @see c2Bundlef
1487  * @see c2rBundle
1488  * @see c2rBundlef
1489  */
1490 template <class T, class U, class V, class W>
1491 void Propagation<T, U, V, W>::parallelProp_EHP(V *cs, V *ct,
1492  W *currents,
1493  int gmode, int ncx, int ncy,
1494  U *res)
1495 {
1496  int final_step;
1497 
1498  for(int n=0; n<numThreads; n++)
1499  {
1500  //std::cout << n << std::endl;
1501  if(n == (numThreads-1))
1502  {
1503  final_step = gt;
1504  }
1505 
1506  else
1507  {
1508  final_step = (n+1) * step;
1509  }
1510 
1511  threadPool[n] = std::thread(&Propagation::propagateBeam_EHP,
1512  this, n * step, final_step,
1513  cs, ct, currents,
1514  gmode, ncx, ncy,
1515  res);
1516  }
1517  joinThreads();
1518 }
1519 
1520 /**
1521  * Calculate scalar field parallel.
1522  *
1523  * Run the propagateScalarBeam function in parallel blocks.
1524  *
1525  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
1526  * @param ct Pointer to reflcontainer or reflcontainerf object containing target grids.
1527  * @param field Pointer to arrC1 or arrC1f object containing source currents.
1528  * @param res Pointer to arrC1 or arrC1f object to be filled with calculation results.
1529  *
1530  * @see reflcontainer
1531  * @see reflcontainerf
1532  * @see arrC1
1533  * @see arrC1f
1534  */
1535 template <class T, class U, class V, class W>
1536 void Propagation<T, U, V, W>::parallelPropScalar(V *cs, V *ct,
1537  W *field,
1538  int gmode, int ncx, int ncy,
1539  U *res)
1540 {
1541  int final_step;
1542 
1543  for(int n=0; n<numThreads; n++)
1544  {
1545  if(n == (numThreads-1))
1546  {
1547  final_step = gt;
1548  }
1549 
1550  else
1551  {
1552  final_step = (n+1) * step;
1553  }
1554 
1555  threadPool[n] = std::thread(&Propagation::propagateScalarBeam,
1556  this, n * step, final_step,
1557  cs, ct, field,
1558  gmode, ncx, ncy,
1559  res);
1560  }
1561  joinThreads();
1562 }
1563 
1564 /**
1565  * Calculate far-field on target.
1566  *
1567  * Calculate integrated E and H fields on a far-field target point.
1568  *
1569  * @param start Index of first loop iteration in parallel block.
1570  * @param stop Index of last loop iteration in parallel block.
1571  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
1572  * @param ct Pointer to reflcontainer or reflcontainerf object containing target grids.
1573  * @param currents Pointer to c2Bundle or c2Bundlef object containing currents on source.
1574  * @param res Pointer to c2Bundle or c2Bundlef object, to be filled with calculation results.
1575  *
1576  * @see reflcontainer
1577  * @see reflcontainerf
1578  * @see c2Bundle
1579  * @see c2Bundlef
1580  */
1581 template <class T, class U, class V, class W>
1582 void Propagation<T, U, V, W>::propagateToFarField(int start, int stop,
1583  V *cs, V *ct,
1584  W *currents,
1585  int gmode, int ncx, int ncy,
1586  U *res)
1587 {
1588  // Scalars (T & complex T)
1589  T theta;
1590  T phi;
1591  T cosEl;
1592 
1593  std::complex<T> e_th;
1594  std::complex<T> e_ph;
1595 
1596  std::complex<T> e_Az;
1597  std::complex<T> e_El;
1598 
1599  // Arrays of Ts
1600  std::array<T, 2> point_ff; // Angular point on far-field
1601  std::array<T, 3> r_hat; // Unit vector in far-field point direction
1602 
1603  // Arrays of complex Ts
1604  std::array<std::array<std::complex<T>, 3>, 2> eh; // far-field EH-field
1605 
1606  int jc = 0;
1607  for(int i=start; i<stop; i++)
1608  {
1609  phi = ct->x[i];
1610  theta = ct->y[i];
1611  cosEl = std::sqrt(1 - sin(theta) * sin(phi) * sin(theta) * sin(phi));
1612 
1613  r_hat[0] = cos(theta) * sin(phi);
1614  r_hat[1] = sin(theta) * sin(phi);
1615  r_hat[2] = cos(phi);
1616 
1617  // Calculate total incoming E and H field at point on target
1618  eh = farfieldAtPoint(cs, currents,
1619  gmode, ncx, ncy,
1620  r_hat);
1621 
1622  res->r1x[i] = eh[0][0].real();
1623  res->r1y[i] = eh[0][1].real();
1624  res->r1z[i] = eh[0][2].real();
1625 
1626  res->i1x[i] = eh[0][0].imag();
1627  res->i1y[i] = eh[0][1].imag();
1628  res->i1z[i] = eh[0][2].imag();
1629 
1630  // TODO: Calculate H far fields
1631  res->r2x[i] = eh[1][0].real();
1632  res->r2y[i] = eh[1][1].real();
1633  res->r2z[i] = eh[1][2].real();
1634 
1635  res->i2x[i] = eh[1][0].imag();
1636  res->i2y[i] = eh[1][1].imag();
1637  res->i2z[i] = eh[1][2].imag();
1638  }
1639 }
1640 
1641 
1642 /**
1643  * Calculate far-field on target.
1644  *
1645  * Calculate integrated E and H fields on a far-field target point.
1646  *
1647  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
1648  * @param currents Pointer to c2Bundle or c2Bundlef object containing currents on source.
1649  * @param r_hat Array of 3 double/float containing target direction angles.
1650  *
1651  * @see reflcontainer
1652  * @see reflcontainerf
1653  * @see c2Bundle
1654  * @see c2Bundlef
1655  */
1656 template <class T, class U, class V, class W>
1657 std::array<std::array<std::complex<T>, 3>, 2> Propagation<T, U, V, W>::farfieldAtPoint(V *cs,
1658  W *currents,
1659  int gmode, int ncx, int ncy,
1660  const std::array<T, 3> &r_hat)
1661 {
1662  int i;
1663  // Scalars (T & complex T)
1664  T norm_dot_R_hat; // Source normal dotted with wavevector direction
1665  T half = 0.5;
1666 
1667  std::complex<T> Green; // Container for Green's function
1668  std::complex<T> js_dot_R; // Magnitude of electric currents on to R_hat
1669  std::complex<T> ms_dot_R; // Magnitude of magnetic currents on to R_hat
1670  T Rprime_dot_R_hat; // Magnitude of source point projected onto R_hat
1671 
1672  // Arrays of Ts
1673  std::array<T, 3> source_point; // Container for xyz co-ordinates
1674  std::array<T, 3> source_norm; // Container for xyz normals.
1675 
1676  // Arrays of complex Ts
1677  std::array<std::complex<T>, 3> e_field; // Electric field surface integral at point
1678  std::array<std::complex<T>, 3> h_field; // Magnetic field surface integral at point
1679  std::array<std::complex<T>, 3> ye_field; // Electric field surface integral at point
1680  std::array<std::complex<T>, 3> yh_field; // Magnetic field surface integral at point
1681  std::array<std::complex<T>, 3> js; // Electric current at source point
1682  std::array<std::complex<T>, 3> ms; // Magnetic current at source point
1683  std::array<std::complex<T>, 3> js_dot_R_R; // Electric currents projected on to R_hat
1684  std::array<std::complex<T>, 3> R_cross_ms; // Cross product of magnetic currents and R_hat
1685  std::array<std::complex<T>, 3> ms_dot_R_R; // Magnetic currents projected on to R_hat
1686  std::array<std::complex<T>, 3> R_cross_js; // Cross product of electric currents and R_hat
1687 
1688 
1689  // Return container
1690  std::array<std::array<std::complex<T>, 3>, 2> e_h_field; // Return container containing e and h-fields
1691 
1692  e_field.fill(z0);
1693  h_field.fill(z0);
1694 
1695  for(int xu=0; xu<ncx; xu++)
1696  {
1697  for (int n=0; n<3; n++)
1698  {
1699  ye_field[n] = z0; // Intermediate electric field due to integral over y/v
1700  yh_field[n] = z0; // Intermediate magnetic field due to integral over y/v
1701  }
1702 
1703  for(int yv=0; yv<ncy; yv++)
1704  {
1705  i = xu*ncy + yv;
1706  source_point[0] = cs->x[i];
1707  source_point[1] = cs->y[i];
1708  source_point[2] = cs->z[i];
1709 
1710  source_norm[0] = cs->nx[i];
1711  source_norm[1] = cs->ny[i];
1712  source_norm[2] = cs->nz[i];
1713 
1714  js[0] = {currents->r1x[i], currents->i1x[i]};
1715  js[1] = {currents->r1y[i], currents->i1y[i]};
1716  js[2] = {currents->r1z[i], currents->i1z[i]};
1717 
1718  ms[0] = {currents->r2x[i], currents->i2x[i]};
1719  ms[1] = {currents->r2y[i], currents->i2y[i]};
1720  ms[2] = {currents->r2z[i], currents->i2z[i]};
1721 
1722  // Check if source point illuminates target point or not.
1723  ut.dot(source_norm, r_hat, norm_dot_R_hat);
1724  //printf("norm_dot_R_hat: %.16g\n", norm_dot_R_hat); // %s is format specifier
1725  if ((norm_dot_R_hat < 0)) {continue;}
1726 
1727  // e-field
1728  ut.dot(js, r_hat, js_dot_R);
1729  ut.s_mult(r_hat, js_dot_R, js_dot_R_R);
1730  ut.ext(r_hat, ms, R_cross_ms);
1731 
1732  // h-field
1733  ut.dot(ms, r_hat, ms_dot_R);
1734  ut.s_mult(r_hat, ms_dot_R, ms_dot_R_R);
1735  ut.ext(r_hat, js, R_cross_js);
1736 
1737  ut.dot(source_point, r_hat, Rprime_dot_R_hat);
1738 
1739  Green = -prefactor * j * exp(-t_direction * j * k * Rprime_dot_R_hat) * cs->area[i];
1740  //printf("%.16g, %.16g\n", Green.real(), Green.imag()); // %s is format specifier
1741 
1742  for( int n=0; n<3; n++)
1743  {
1744  if ((gmode != 1) && ((yv==0) || (yv==ncy-1)))
1745  {
1746  ye_field[n] += half*(js[n] - js_dot_R_R[n] + t_direction * R_cross_ms[n]) * Green;
1747  yh_field[n] += half*(ms[n] - ms_dot_R_R[n] - t_direction * R_cross_js[n]) * Green;
1748  }
1749  else
1750  {
1751  ye_field[n] += (js[n] - js_dot_R_R[n] + t_direction * R_cross_ms[n]) * Green;
1752  yh_field[n] += (ms[n] - ms_dot_R_R[n] - t_direction * R_cross_js[n]) * Green;
1753  }
1754  }
1755 
1756  } // End of yv loop
1757  for (int n=0; n<3; n++)
1758  {
1759  if ((xu==0) || (xu==ncx-1)) // Only add half the point value
1760  {
1761  e_field[n] += half*ye_field[n];
1762  h_field[n] += half*yh_field[n];
1763  }
1764  else
1765  {
1766  e_field[n] += ye_field[n];
1767  h_field[n] += yh_field[n];
1768  }
1769  }
1770 
1771  } // End of xu loop
1772 
1773 
1774  // Pack e and h together in single container
1775  e_h_field[0] = e_field;
1776  e_h_field[1] = h_field;
1777 
1778  //std::cout << ut.abs(e_field) << std::endl;
1779 
1780  return e_h_field;
1781 }
1782 
1783 /**
1784  * Calculate far-field on target.
1785  *
1786  * Calculate integrated E and H fields on a far-field target point.
1787  *
1788  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
1789  * @param currents Pointer to c2Bundle or c2Bundlef object containing currents on source.
1790  * @param r_hat Array of 3 double/float containing target direction angles.
1791  *
1792  * @see reflcontainer
1793  * @see reflcontainerf
1794  * @see c2Bundle
1795  * @see c2Bundlef
1796  */
1797 template <class T, class U, class V, class W>
1798 std::array<std::array<std::complex<T>, 3>, 2> Propagation<T, U, V, W>::farfieldAtPoint_Old(V *cs,
1799  W *currents,
1800  const std::array<T, 3> &r_hat)
1801 {
1802  // Scalars (T & complex T)
1803  T omega_mu; // Angular frequency of field times mu
1804  T omega_eps; // Angular frequency of field times eps
1805  T r_hat_in_rp; // r_hat dot product r_prime
1806  std::complex<T> r_in_s; // Container for inner products between wavevctor and currents
1807  std::complex<T> expo;
1808  T area;
1809 
1810  // Arrays of Ts
1811  std::array<T, 3> source_point; // Container for xyz co-ordinates
1812 
1813  // Arrays of complex Ts
1814  std::array<std::complex<T>, 3> e; // Electric field on far-field point
1815  std::array<std::complex<T>, 3> h; // Magnetic field on far-field point
1816  std::array<std::complex<T>, 3> _js; // Temporary Electric current at source point
1817  std::array<std::complex<T>, 3> _ms; // Temporary Magnetic current at source point
1818 
1819  std::array<std::complex<T>, 3> js; // Build radiation integral
1820  std::array<std::complex<T>, 3> ms; // Build radiation integral
1821 
1822  std::array<std::complex<T>, 3> _ctemp;
1823  std::array<std::complex<T>, 3> js_tot_factor;
1824  std::array<std::complex<T>, 3> ms_tot_factor;
1825  std::array<std::complex<T>, 3> js_tot_factor_h;
1826  std::array<std::complex<T>, 3> ms_tot_factor_h;
1827 
1828  // Output array
1829  std::array<std::array<std::complex<T>, 3>, 2> out;
1830 
1831  // Matrices
1832  std::array<std::array<T, 3>, 3> rr_dyad; // Dyadic product between r_hat - r_hat
1833  std::array<std::array<T, 3>, 3> eye_min_rr; // I - rr
1834 
1835  e.fill(z0);
1836  h.fill(z0);
1837  js.fill(z0);
1838  ms.fill(z0);
1839 
1840  omega_mu = C_L * k * MU_0;
1841 
1842  ut.dyad(r_hat, r_hat, rr_dyad);
1843  ut.matDiff(this->eye, rr_dyad, eye_min_rr);
1844 
1845  for(int i=0; i<gs; i++)
1846  {
1847  source_point[0] = cs->x[i];
1848  source_point[1] = cs->y[i];
1849  source_point[2] = cs->z[i];
1850 
1851  ut.dot(r_hat, source_point, r_hat_in_rp);
1852 
1853  expo = exp(j * k * r_hat_in_rp);
1854  area = cs->area[i];
1855 
1856  _js[0] = {currents->r1x[i], currents->i1x[i]};
1857  _js[1] = {currents->r1y[i], currents->i1y[i]};
1858  _js[2] = {currents->r1z[i], currents->i1z[i]};
1859 
1860  _ms[0] = {currents->r2x[i], currents->i2x[i]};
1861  _ms[1] = {currents->r2y[i], currents->i2y[i]};
1862  _ms[2] = {currents->r2z[i], currents->i2z[i]};
1863 
1864  for (int n=0; n<3; n++)
1865  {
1866  js[n] += _js[n] * expo * area;
1867  ms[n] += _ms[n] * expo * area;
1868  }
1869  }
1870 
1871  ut.matVec(eye_min_rr, js, _ctemp);
1872  ut.s_mult(_ctemp, omega_mu, js_tot_factor);
1873 
1874  ut.ext(r_hat, ms, _ctemp);
1875  ut.s_mult(_ctemp, k, ms_tot_factor);
1876 
1877  omega_eps = C_L * k * EPS;
1878 
1879  ut.matVec(eye_min_rr, ms, _ctemp);
1880  ut.s_mult(_ctemp, omega_eps, ms_tot_factor_h);
1881 
1882  ut.ext(r_hat, js, _ctemp);
1883  ut.s_mult(_ctemp, k, js_tot_factor_h);
1884 
1885  for (int n=0; n<3; n++)
1886  {
1887  e[n] = -js_tot_factor[n] + ms_tot_factor[n];
1888  h[n] = -ms_tot_factor[n] - js_tot_factor[n];
1889  }
1890 
1891  out[0] = e;
1892  out[1] = h;
1893  return out;
1894 }
1895 
1896 /**
1897  * Calculate far-field parallel.
1898  *
1899  * Run the propagateToFarField function in parallel blocks.
1900  *
1901  * @param cs Pointer to reflcontainer or reflcontainerf object containing source grids.
1902  * @param ct Pointer to reflcontainer or reflcontainerf object containing target grids.
1903  * @param currents Pointer to c2Bundle or c2Bundlef object containing source currents.
1904  * @param res Pointer to c2Bundle or c2Bundlef object to be filled with calculation results.
1905  *
1906  * @see reflcontainer
1907  * @see reflcontainerf
1908  * @see c2Bundle
1909  * @see c2Bundlef
1910  */
1911 template <class T, class U, class V, class W>
1912 void Propagation<T, U, V, W>::parallelFarField(V *cs, V *ct,
1913  W *currents,
1914  int gmode, int ncx, int ncy,
1915  U *res)
1916 {
1917  int final_step;
1918 
1919  for(int n=0; n<numThreads; n++)
1920  {
1921  if(n == (numThreads-1))
1922  {
1923  final_step = gt;
1924  }
1925 
1926  else
1927  {
1928  final_step = (n+1) * step;
1929  }
1930 
1931  threadPool[n] = std::thread(&Propagation::propagateToFarField,
1932  this, n * step, final_step,
1933  cs, ct, currents,
1934  gmode, ncx, ncy,
1935  res);
1936  }
1937  joinThreads();
1938 }
1939 
1940 template <class T, class U, class V, class W>
1941 void Propagation<T, U, V, W>::joinThreads()
1942 {
1943  for (std::thread &t : threadPool)
1944  {
1945  if (t.joinable())
1946  {
1947  t.join();
1948  }
1949  }
1950 }
1951 
1952 template <class T, class U, class V, class W>
1953 void Propagation<T, U, V, W>::_debugArray(T *arr, int idx)
1954 {
1955  T toPrint = arr[idx];
1956  std::cout << "Value of c-style array, element " << idx << ", is : " << toPrint << std::endl;
1957 }
1958 
1959 template <class T, class U, class V, class W>
1960 void Propagation<T, U, V, W>::_debugArray(std::array<T, 3> arr)
1961 {
1962  std::cout << arr[0] << ", " << arr[1] << ", " << arr[2] << std::endl;
1963 }
1964 
1965 template <class T, class U, class V, class W>
1966 void Propagation<T, U, V, W>::_debugArray(std::array<std::complex<T>, 3> arr)
1967 {
1968  std::cout << arr[0].real() << " + " << arr[0].imag() << "j"
1969  << ", " << arr[1].real() << " + " << arr[1].imag() << "j"
1970  << ", " << arr[2].real() << " + " << arr[2].imag() << "j"
1971  << std::endl;
1972 }
1973 
1974 
1975 #endif
Propagation::propagateScalarBeam
void propagateScalarBeam(int start, int stop, V *cs, V *ct, W *field, int gmode, int ncx, int ncy, U *res)
Definition: Propagation.h:934
_debugArray
__host__ __device__ void _debugArray(cuFloatComplex arr[3])
Definition: Debug.h:23
Propagation::parallelProp_EH
void parallelProp_EH(V *cs, V *ct, W *currents, int gmode, int ncx, int ncy, U *res)
Definition: Propagation.h:1396
Propagation::threadPool
std::vector< std::thread > threadPool
Definition: Propagation.h:69
Propagation::propagateBeam_EHP
void propagateBeam_EHP(int start, int stop, V *cs, V *ct, W *currents, int gmode, int ncx, int ncy, U *res)
Definition: Propagation.h:806
Propagation::propagateBeam_JM_PEC
void propagateBeam_JM_PEC(int start, int stop, V *cs, V *ct, W *currents, int gmode, int ncx, int ncy, U *res)
Definition: Propagation.h:270
Propagation::farfieldAtPoint_Old
std::array< std::array< std::complex< T >, 3 >, 2 > farfieldAtPoint_Old(V *cs, W *currents, const std::array< T, 3 > &point_target)
Definition: Propagation.h:1798
Propagation::parallelFarField
void parallelFarField(V *cs, V *ct, W *currents, int gmode, int ncx, int ncy, U *res)
Definition: Propagation.h:1912
Propagation::propagateBeam_JM
void propagateBeam_JM(int start, int stop, V *cs, V *ct, W *currents, int gmode, int ncx, int ncy, U *res)
Definition: Propagation.h:343
Propagation::parallelProp_JMEH
void parallelProp_JMEH(V *cs, V *ct, W *currents, int gmode, int ncx, int ncy, U *res)
Definition: Propagation.h:1445
Propagation::parallelProp_EHP
void parallelProp_EHP(V *cs, V *ct, W *currents, int gmode, int ncx, int ncy, U *res)
Definition: Propagation.h:1491
Structs.h
Structs used within PyPO.
Utils
Definition: Utils.h:22
Propagation::Propagation
Propagation(T omega, int numThreads, int gs, int gt, T epsilon, T t_direction, bool verbose=false)
Definition: Propagation.h:194
Propagation::parallelPropScalar
void parallelPropScalar(V *cs, V *ct, W *field, int gmode, int ncx, int ncy, U *res)
Definition: Propagation.h:1536
Propagation::propagateBeam_JMEH
void propagateBeam_JMEH(int start, int stop, V *cs, V *ct, W *currents, int gmode, int ncx, int ncy, U *res)
Definition: Propagation.h:639
Utils.h
Linear algebra functions for the CPU version of PyPO.
Propagation::fieldScalarAtPoint
std::complex< T > fieldScalarAtPoint(V *cs, W *field, int gmode, int ncx, int ncy, const std::array< T, 3 > &point_target)
Definition: Propagation.h:1280
Propagation::parallelProp_JM
void parallelProp_JM(V *cs, V *ct, W *currents, int gmode, int ncx, int ncy, U *res)
Definition: Propagation.h:1349
Propagation::fieldAtPoint_Old
std::array< std::array< std::complex< T >, 3 >, 2 > fieldAtPoint_Old(V *cs, W *currents, const std::array< T, 3 > &point_target)
Definition: Propagation.h:1162
Propagation::farfieldAtPoint
std::array< std::array< std::complex< T >, 3 >, 2 > farfieldAtPoint(V *cs, W *currents, int gmode, int ncx, int ncy, const std::array< T, 3 > &point_target)
Definition: Propagation.h:1657
farfieldAtPoint
__device__ void farfieldAtPoint(float *d_xs, float *d_ys, float *d_zs, float *d_nxs, float *d_nys, float *d_nzs, cuFloatComplex *d_Jx, cuFloatComplex *d_Jy, cuFloatComplex *d_Jz, cuFloatComplex *d_Mx, cuFloatComplex *d_My, cuFloatComplex *d_Mz, float(&r_hat)[3], float *d_A, int gmode, int ncx, int ncy, cuFloatComplex(&e)[3], cuFloatComplex(&h)[3])
Definition: Kernelsf.cu:1115
Propagation::ut
Utils< T > ut
Definition: Propagation.h:74
fieldAtPoint
__device__ void fieldAtPoint(float *d_xs, float *d_ys, float *d_zs, float *d_nxs, float *d_nys, float *d_nzs, cuFloatComplex *d_Jx, cuFloatComplex *d_Jy, cuFloatComplex *d_Jz, cuFloatComplex *d_Mx, cuFloatComplex *d_My, cuFloatComplex *d_Mz, float(&point)[3], float *d_A, int gmode, int ncx, int ncy, cuFloatComplex(&d_ei)[3], cuFloatComplex(&d_hi)[3])
Definition: Kernelsf.cu:116
Propagation
Definition: Propagation.h:34
Propagation::propagateBeam_JMEH_PEC
void propagateBeam_JMEH_PEC(int start, int stop, V *cs, V *ct, W *currents, int gmode, int ncx, int ncy, U *res)
Definition: Propagation.h:483
Propagation::propagateToFarField
void propagateToFarField(int start, int stop, V *cs, V *ct, W *currents, int gmode, int ncx, int ncy, U *res)
Definition: Propagation.h:1582
InterfaceReflector.h
Header for reflector generation interface.
Propagation::propagateBeam_EH
void propagateBeam_EH(int start, int stop, V *cs, V *ct, W *currents, int gmode, int ncx, int ncy, U *res)
Definition: Propagation.h:575
Propagation::fieldAtPoint
std::array< std::array< std::complex< T >, 3 >, 2 > fieldAtPoint(V *cs, W *currents, int gmode, int ncx, int ncy, const std::array< T, 3 > &point_target)
Definition: Propagation.h:978