PyPO User Manual
BeamInit.h
Go to the documentation of this file.
1 #include <iostream>
2 #include <array>
3 #include <cmath>
4 #include <complex>
5 
6 #include "Utils.h"
7 #include "Structs.h"
8 #include "InterfaceReflector.h"
9 #define M_PI 3.14159265358979323846
10 
11 #ifndef __BeamInit_h
12 #define __BeamInit_h
13 
14 /*! \file BeamInit.h
15  \brief Initialize beam objects.
16 
17  Initializes ray-trace frames, PO beams, and custom beams.
18 */
19 
20 /**
21  * Initialize ray-trace frame from RTDict or RTDictf.
22  *
23  * Takes an RTDict or RTDictf and generates a frame object, which can be used
24  * to initialize a ray-trace.
25  *
26  * @param rdict RTDict or RTDictf object from which to generate a frame.
27  * @param fr Pointer to cframe or cframef object.
28  *
29  * @see RTDict
30  * @see RTDictf
31  * @see cframe
32  * @see cframef
33  */
34 template<typename T, typename U, typename V>
35 void initFrame(T rdict, U *fr);
36 
37 /**
38  * Initialize Gaussian ray-trace frame from RTDict or RTDictf.
39  *
40  * Takes a GRTDict or GRTDictf and generates a frame object, which can be used
41  * to initialize a Gaussian ray-trace.
42  *
43  * @param grdict GRTDict or GRTDictf object from which to generate a frame.
44  * @param fr Pointer to cframe or cframef object.
45  *
46  * @see GRTDict
47  * @see GRTDictf
48  * @see cframe
49  * @see cframef
50  */
51 template<typename T, typename U, typename V>
52 void initRTGauss(T grdict, U *fr);
53 
54 /**
55  * Probability density function for drawing ray positions and tilts from Gaussian.
56  *
57  * The density function is used for generating Gaussian ray-trace beams.
58  * Using rejection sampling, the method draws samples from the Gaussian pdf.
59  *
60  * @param vars Vector of length 4, containing the xy positions and tilts of the ray to be checked.
61  * @param scales Vector of length 4 containing the scale factors along xy and tilts
62  */
63 template<typename T>
64 T pdfGauss(std::vector<T> vars, std::vector<T> scales);
65 
66 /**
67  * Initialize Gaussian beam from GPODict or GPODictf.
68  *
69  * The generated Gaussian has E-field and magnetic currents only. This _is_ what you want
70  * as a source in PO calculations.
71  *
72  * Takes a GPODict or GPODictf and generates two c2Bundle or c2Bundlef objects, which contain the field and
73  * associated currents and are allocated to passed pointer arguments.
74  *
75  * @param gdict GPODict or GPODictf object from which to generate a Gaussian beam.
76  * @param refldict reflparams or reflparamsf object corresponding to surface on
77  * which to generate the Gaussian beam.
78  * @param res_field Pointer to c2Bundle or c2Bundlef object.
79  * @param res_current Pointer to c2Bundle or c2Bundlef object.
80  *
81  * @see GPODict
82  * @see GPODictf
83  * @see reflparams
84  * @see reflparamsf
85  * @see c2Bundle
86  * @see c2Bundlef
87  */
88 template<typename T, typename U, typename V, typename W, typename G>
89 void initGauss(T gdict, U refldict, V *res_field, V *res_current);
90 
91 
92 /**
93  * Initialize scalar Gaussian beam from GPODict or GPODictf.
94  *
95  * Takes a ScalarGPODict or ScalarGPODictf and generates an arrC1 or arrC1f object.
96  *
97  * @param gdict ScalarGPODict or ScalarGPODictf object from which to generate a Gaussian beam.
98  * @param refldict reflparams or reflparamsf object corresponding to surface on
99  * which to generate the Gaussian beam.
100  * @param res_field Pointer to arrC1 or arrC1f object.
101  *
102  * @see ScalarGPODict
103  * @see ScalarGPODictf
104  * @see reflparams
105  * @see reflparamsf
106  * @see arrC1
107  * @see arrC1f
108  */
109 template<typename T, typename U, typename V, typename W, typename G>
110 void initScalarGauss(T sgdict, U refldict, V *res_field);
111 
112 /**
113  * Calculate currents from electromagnetic field.
114  *
115  * Calculate the J and M vectorial currents given a vectorial E and H field.
116  * Can calculate full currents, PMC and PEC surfaces.
117  *
118  * @param res_field Pointer to c2Bundle or c2Bundlef object.
119  * @param res_current Pointer to c2Bundle or c2Bundlef object.
120  * @param refldict reflparams or reflparamsf object corresponding to surface on
121  * which to calculate currents.
122  * @param mode How to calculate currents. 0 is full currents, 1 is PMC and 2 is PEC.
123  *
124  * @see c2Bundle
125  * @see c2Bundlef
126  * @see reflparams
127  * @see reflparamsf
128  */
129 template<typename T, typename U, typename V, typename W>
130 void calcJM(T *res_field, T *res_current, V refldict, int mode);
131 
132 template<typename T, typename U, typename V>
133 void initFrame(T rdict, U *fr)
134 {
135  std::array<V, 3> nomChief = {0, 0, 1};
136  std::array<V, 3> zero = {0, 0, 0};
137 
138  Utils<V> ut;
139 
140  int nTot = 1 + rdict.nRing * 4 * rdict.nRays;
141 
142  fr->size = nTot;
143 
144  V alpha = 0; // Set first circle ray in the right of beam
145  V d_alpha = 0;
146 
147  if (rdict.nRays > 0) {d_alpha = 2 * M_PI / (4 * rdict.nRays);}
148 
149  int n = 1;
150 
151  std::array<V, 3> rotation;
152  std::array<V, 3> direction;
153 
154  fr->x[0] = 0.;
155  fr->y[0] = 0.;
156  fr->z[0] = 0.;
157 
158  fr->dx[0] = 0.;
159  fr->dy[0] = 0.;
160  fr->dz[0] = 1.;
161 
162  std::array<V, 3> pos;
163 
164  for (int i=1; i<nTot; i++)
165  {
166  fr->x[i] = rdict.x0 * cos(alpha) / rdict.nRing * n;
167  fr->y[i] = rdict.y0 * sin(alpha) / rdict.nRing * n;
168  fr->z[i] = 0.;
169 
170  rotation[0] = rdict.angy0 * sin(alpha) / rdict.nRing * n;
171  rotation[1] = rdict.angx0 * cos(alpha) / rdict.nRing * n;
172  rotation[2] = 2 * alpha;
173 
174  ut.matRot(rotation, nomChief, zero, direction);
175 
176  fr->dx[i] = direction[0];
177  fr->dy[i] = direction[1];
178  fr->dz[i] = direction[2];
179  alpha += d_alpha;
180 
181  if (i == int(nTot / rdict.nRing) * n)
182  {
183  n += 1;
184  alpha = 0;
185  }
186  }
187 }
188 
189 template<typename T, typename U, typename V>
190 void initRTGauss(T grdict, U *fr)
191 {
192  V thres = 3.; // Choose 3-sigma point
193  std::array<V, 3> nomChief = {0, 0, 1};
194  std::array<V, 3> zero = {0, 0, 0};
195 
196  Utils<V> ut;
197  Random<V> rando;
198  if (grdict.seed != -1) {Random<V> rando(grdict.seed);}
199 
200  fr->size = grdict.nRays;
201 
202  std::array<V, 3> rotation;
203  std::array<V, 3> direction;
204  std::array<V, 3> pos;
205 
206  // Initialize scale vector
207  std::vector<V> scales{grdict.x0, grdict.y0, grdict.angx0, grdict.angy0};
208 
209  fr->x[0] = 0.;
210  fr->y[0] = 0.;
211  fr->z[0] = 0.;
212 
213  fr->dx[0] = 0.;
214  fr->dy[0] = 0.;
215  fr->dz[0] = 1.;
216 
217  // Start rejection sampling. Use n_suc as succes counter
218  int n_suc = 1;
219  V yi;
220  std::vector<V> xi(4, 0);
221  V lower = 0.0;
222 
223  while (n_suc < grdict.nRays)
224  {
225  // First, generate y-value between 0 and 1.
226  yi = rando.generateUniform(lower);
227 
228  // Now, generate vector of xi
229  xi = rando.generateUniform(grdict.nRays);
230 
231  for (int k = 0; k<4; k++) {xi[k] = xi[k] * thres * scales[k];}
232  if (pdfGauss<V>(xi, scales) > yi)
233  {
234  // Rotate chief ray by tilt angles found
235  rotation = {xi[3], xi[2], 0};
236  ut.matRot(rotation, nomChief, zero, direction);
237  pos = {xi[0], xi[1], 0};
238 
239  fr->x[n_suc] = pos[0];
240  fr->y[n_suc] = pos[1];
241  fr->z[n_suc] = pos[2];
242 
243  fr->dx[n_suc] = direction[0];
244  fr->dy[n_suc] = direction[1];
245  fr->dz[n_suc] = direction[2];
246  n_suc++;
247  }
248  }
249 }
250 
251 template<typename T>
252 T pdfGauss(std::vector<T> vars, std::vector<T> scales)
253 {
254  T norm = 1 / (M_PI*M_PI * scales[0] * scales[1] * scales[2] * scales[3]);
255 
256  return norm * exp(-vars[0]*vars[0] / (scales[0]*scales[0])) * exp(-vars[1]*vars[1] / (scales[1]*scales[1])) *
257  exp(-vars[2]*vars[2] / (scales[2]*scales[2])) * exp(-vars[3]*vars[3] / (scales[3]*scales[3]));
258 }
259 
260 
261 template<typename T, typename U, typename V, typename W, typename G>
262 void initGauss(T gdict, U refldict, V *res_field, V *res_current)
263 {
264  int nTot = refldict.n_cells[0] * refldict.n_cells[1];
265 
266  W reflc;
267 
268  reflc.size = nTot;
269 
270  reflc.x = new G[nTot];
271  reflc.y = new G[nTot];
272  reflc.z = new G[nTot];
273 
274  reflc.nx = new G[nTot];
275  reflc.ny = new G[nTot];
276  reflc.nz = new G[nTot];
277 
278  reflc.area = new G[nTot];
279 
280  Utils<G> ut;
281 
282  bool transform = true;
283  generateGrid(refldict, &reflc, transform);
284 
285 
286  G zRx = M_PI * gdict.w0x*gdict.w0x * gdict.n / gdict.lam;
287  G zRy = M_PI * gdict.w0y*gdict.w0y * gdict.n / gdict.lam;
288  G k = 2 * M_PI / gdict.lam;
289 
290  G wzx;
291  G wzy;
292  G Rzx_inv;
293  G Rzy_inv;
294  G phizx;
295  G phizy;
296 
297  std::complex<G> j(0, 1);
298 
299  std::complex<G> field_atPoint;
300  std::array<std::complex<G>, 3> efield;
301  std::array<G, 3> n_source;
302  std::array<std::complex<G>, 3> m;
303 
304  for (int i=0; i<nTot; i++)
305  {
306  wzx = gdict.w0x * sqrt(1 + (reflc.z[i] / zRx)*(reflc.z[i] / zRx));
307  wzy = gdict.w0y * sqrt(1 + ((reflc.z[i] - gdict.dxyz) / zRy)*((reflc.z[i] - gdict.dxyz) / zRy));
308  Rzx_inv = reflc.z[i] / (reflc.z[i]*reflc.z[i] + zRx*zRx);
309  Rzy_inv = (reflc.z[i] - gdict.dxyz) / ((reflc.z[i] - gdict.dxyz)*(reflc.z[i] - gdict.dxyz) + zRy*zRy);
310  phizx = atan(reflc.z[i] / zRx);
311  phizy = atan((reflc.z[i] - gdict.dxyz) / zRy);
312 
313  //field_atPoint = gdict.E0 * sqrt(2 / (M_PI * wzx * wzy)) * exp(-(reflc.x[i]/wzx)*(reflc.x[i]/wzx) - (reflc.y[i]/wzy)*(reflc.y[i]/wzy) -
314  // j*M_PI/gdict.lam * (reflc.x[i]*reflc.x[i]*Rzx_inv + reflc.y[i]*reflc.y[i]*Rzy_inv) - j*k*reflc.z[i] + j*(phizx - phizy)*0.5);
315 
316  field_atPoint = gdict.E0 * exp(-(reflc.x[i]/wzx)*(reflc.x[i]/wzx) - (reflc.y[i]/wzy)*(reflc.y[i]/wzy) -
317  j*M_PI/gdict.lam * (reflc.x[i]*reflc.x[i]*Rzx_inv + reflc.y[i]*reflc.y[i]*Rzy_inv) - j*k*reflc.z[i] + j*(phizx - phizy)*0.5);
318 
319  efield[0] = field_atPoint * gdict.pol[0];
320  efield[1] = field_atPoint * gdict.pol[1];
321  efield[2] = field_atPoint * gdict.pol[2];
322 
323  n_source[0] = reflc.nx[i];
324  n_source[1] = reflc.ny[i];
325  n_source[2] = reflc.nz[i];
326 
327  ut.ext(n_source, efield, m);
328 
329  res_field->r1x[i] = efield[0].real();
330  res_field->i1x[i] = efield[0].imag();
331 
332  res_field->r1y[i] = efield[1].real();
333  res_field->i1y[i] = efield[1].imag();
334 
335  res_field->r1z[i] = efield[2].real();
336  res_field->i1z[i] = efield[2].imag();
337 
338  // Set H to zero
339  res_field->r2x[i] = 0;
340  res_field->i2x[i] = 0;
341 
342  res_field->r2y[i] = 0;
343  res_field->i2y[i] = 0;
344 
345  res_field->r2z[i] = 0;
346  res_field->i2z[i] = 0;
347 
348  // Fill currents
349  res_current->r1x[i] = 0;
350  res_current->i1x[i] = 0;
351 
352  res_current->r1y[i] = 0;
353  res_current->i1y[i] = 0;
354 
355  res_current->r1z[i] = 0;
356  res_current->i1z[i] = 0;
357 
358  // Set H to zero
359  res_current->r2x[i] = -2*m[0].real();
360  res_current->i2x[i] = -2*m[0].imag();
361 
362  res_current->r2y[i] = -2*m[1].real();
363  res_current->i2y[i] = -2*m[1].imag();
364 
365  res_current->r2z[i] = -2*m[2].real();
366  res_current->i2z[i] = -2*m[2].imag();
367  }
368 
369  delete reflc.x;
370  delete reflc.y;
371  delete reflc.z;
372 
373  delete reflc.nx;
374  delete reflc.ny;
375  delete reflc.nz;
376 
377  delete reflc.area;
378 }
379 
380 template<typename T, typename U, typename V, typename W, typename G>
381 void initScalarGauss(T sgdict, U refldict, V *res_field)
382 {
383  int nTot = refldict.n_cells[0] * refldict.n_cells[1];
384  W reflc;
385 
386  reflc.size = nTot;
387  reflc.x = new G[nTot];
388  reflc.y = new G[nTot];
389  reflc.z = new G[nTot];
390 
391  reflc.nx = new G[nTot];
392  reflc.ny = new G[nTot];
393  reflc.nz = new G[nTot];
394 
395  reflc.area = new G[nTot];
396  Utils<G> ut;
397 
398  bool transform = true;
399  generateGrid(refldict, &reflc, transform);
400 
401  G zRx = M_PI * sgdict.w0x*sgdict.w0x * sgdict.n / sgdict.lam;
402  G zRy = M_PI * sgdict.w0y*sgdict.w0y * sgdict.n / sgdict.lam;
403  G k = 2 * M_PI / sgdict.lam;
404 
405  G wzx;
406  G wzy;
407  G Rzx_inv;
408  G Rzy_inv;
409  G phizx;
410  G phizy;
411  std::complex<G> j(0, 1);
412 
413  std::complex<G> efield;
414 
415  for (int i=0; i<nTot; i++)
416  {
417  wzx = sgdict.w0x * sqrt(1 + (reflc.z[i] / zRx)*(reflc.z[i] / zRx));
418  wzy = sgdict.w0y * sqrt(1 + ((reflc.z[i] - sgdict.dxyz) / zRy)*((reflc.z[i] - sgdict.dxyz) / zRy));
419  Rzx_inv = reflc.z[i] / (reflc.z[i]*reflc.z[i] + zRx*zRx);
420  Rzy_inv = (reflc.z[i] - sgdict.dxyz) / ((reflc.z[i] - sgdict.dxyz)*(reflc.z[i] - sgdict.dxyz) + zRy*zRy);
421  phizx = atan(reflc.z[i] / zRx);
422  phizy = atan((reflc.z[i] - sgdict.dxyz) / zRy);
423 
424  efield = sgdict.E0 * sqrt(2 / (M_PI * wzx * wzy)) * exp(-(reflc.x[i]/wzx)*(reflc.x[i]/wzx) - (reflc.y[i]/wzy)*(reflc.y[i]/wzy) -
425  j*M_PI/sgdict.lam * (reflc.x[i]*reflc.x[i]*Rzx_inv + reflc.y[i]*reflc.y[i]*Rzy_inv) - j*k*reflc.z[i] + j*(phizx - phizy)*0.5);
426 
427  res_field->x[i] = efield.real();
428  res_field->y[i] = efield.imag();
429  }
430 
431  delete reflc.x;
432  delete reflc.y;
433  delete reflc.z;
434 
435  delete reflc.nx;
436  delete reflc.ny;
437  delete reflc.nz;
438 
439  delete reflc.area;
440 }
441 
442 template<typename T, typename U, typename V, typename W>
443 void calcJM(T *res_field, T *res_current, V refldict, int mode)
444 {
445  W reflc;
446  int nTot = refldict.n_cells[0] * refldict.n_cells[1];
447 
448  reflc.size = nTot;
449 
450  reflc.x = new U[nTot];
451  reflc.y = new U[nTot];
452  reflc.z = new U[nTot];
453 
454  reflc.nx = new U[nTot];
455  reflc.ny = new U[nTot];
456  reflc.nz = new U[nTot];
457 
458  reflc.area = new U[nTot];
459 
460  bool transform = true;
461  generateGrid(refldict, &reflc, transform);
462 
463  Utils<U> ut;
464 
465  std::array<std::complex<U>, 3> field;
466  std::array<U, 3> n_source;
467 
468  std::array<std::complex<U>, 3> js;
469  std::array<std::complex<U>, 3> ms;
470 
471  // full currents
472  if (mode == 0)
473  {
474  for (int i=0; i<nTot; i++)
475  {
476  field[0] = {res_field->r1x[i], res_field->i1x[i]};
477  field[1] = {res_field->r1y[i], res_field->i1y[i]};
478  field[2] = {res_field->r1z[i], res_field->i1z[i]};
479 
480  n_source[0] = reflc.nx[i];
481  n_source[1] = reflc.ny[i];
482  n_source[2] = reflc.nz[i];
483 
484  ut.ext(n_source, field, ms);
485 
486  res_current->r2x[i] = -ms[0].real();
487  res_current->i2x[i] = -ms[0].imag();
488 
489  res_current->r2y[i] = -ms[1].real();
490  res_current->i2y[i] = -ms[1].imag();
491 
492  res_current->r2z[i] = -ms[2].real();
493  res_current->i2z[i] = -ms[2].imag();
494 
495  field[0] = {res_field->r2x[i], res_field->i2x[i]};
496  field[1] = {res_field->r2y[i], res_field->i2y[i]};
497  field[2] = {res_field->r2z[i], res_field->i2z[i]};
498 
499  ut.ext(n_source, field, js);
500 
501  res_current->r1x[i] = js[0].real();
502  res_current->i1x[i] = js[0].imag();
503 
504  res_current->r1y[i] = js[1].real();
505  res_current->i1y[i] = js[1].imag();
506 
507  res_current->r1z[i] = js[2].real();
508  res_current->i1z[i] = js[2].imag();
509  }
510  }
511 
512  // PMC
513  else if (mode == 1)
514  {
515  for (int i=0; i<nTot; i++)
516  {
517  field[0] = {res_field->r1x[i], res_field->i1x[i]};
518  field[1] = {res_field->r1y[i], res_field->i1y[i]};
519  field[2] = {res_field->r1z[i], res_field->i1z[i]};
520 
521  n_source[0] = reflc.nx[i];
522  n_source[1] = reflc.ny[i];
523  n_source[2] = reflc.nz[i];
524 
525  ut.ext(n_source, field, ms);
526 
527  res_current->r2x[i] = -2*ms[0].real();
528  res_current->i2x[i] = -2*ms[0].imag();
529 
530  res_current->r2y[i] = -2*ms[1].real();
531  res_current->i2y[i] = -2*ms[1].imag();
532 
533  res_current->r2z[i] = -2*ms[2].real();
534  res_current->i2z[i] = -2*ms[2].imag();
535 
536  res_current->r1x[i] = 0;
537  res_current->i1x[i] = 0;
538 
539  res_current->r1y[i] = 0;
540  res_current->i1y[i] = 0;
541 
542  res_current->r1z[i] = 0;
543  res_current->i1z[i] = 0;
544  }
545  }
546 
547  // PEC
548  else if (mode == 2)
549  {
550  for (int i=0; i<nTot; i++)
551  {
552  field[0] = {res_field->r2x[i], res_field->i2x[i]};
553  field[1] = {res_field->r2y[i], res_field->i2y[i]};
554  field[2] = {res_field->r2z[i], res_field->i2z[i]};
555 
556  n_source[0] = reflc.nx[i];
557  n_source[1] = reflc.ny[i];
558  n_source[2] = reflc.nz[i];
559 
560  ut.ext(n_source, field, js);
561 
562  res_current->r1x[i] = 2*js[0].real();
563  res_current->i1x[i] = 2*js[0].imag();
564 
565  res_current->r1y[i] = 2*js[1].real();
566  res_current->i1y[i] = 2*js[1].imag();
567 
568  res_current->r1z[i] = 2*js[2].real();
569  res_current->i1z[i] = 2*js[2].imag();
570 
571  res_current->r2x[i] = 0;
572  res_current->i2x[i] = 0;
573 
574  res_current->r2y[i] = 0;
575  res_current->i2y[i] = 0;
576 
577  res_current->r2z[i] = 0;
578  res_current->i2z[i] = 0;
579  }
580  }
581 
582  delete reflc.x;
583  delete reflc.y;
584  delete reflc.z;
585 
586  delete reflc.nx;
587  delete reflc.ny;
588  delete reflc.nz;
589 
590  delete reflc.area;
591 }
592 #endif
initScalarGauss
void initScalarGauss(T sgdict, U refldict, V *res_field)
Definition: BeamInit.h:381
pdfGauss
T pdfGauss(std::vector< T > vars, std::vector< T > scales)
Definition: BeamInit.h:252
initFrame
void initFrame(T rdict, U *fr)
Definition: BeamInit.h:133
calcJM
void calcJM(T *res_field, T *res_current, V refldict, int mode)
Definition: BeamInit.h:443
initGauss
void initGauss(T gdict, U refldict, V *res_field, V *res_current)
Definition: BeamInit.h:262
initRTGauss
void initRTGauss(T grdict, U *fr)
Definition: BeamInit.h:190
generateGrid
void generateGrid(reflparams refl, reflcontainer *container, bool transform, bool spheric)
Definition: InterfaceReflector.cpp:880
InterfaceReflector.h
Header for reflector generation interface.
Structs.h
Structs used within PyPO.
Utils.h
Linear algebra functions for the CPU version of PyPO.
Random
Definition: Random.h:26
Random::generateUniform
T generateUniform(T lower=-1.0)
Definition: Random.h:77
Utils::ext
void ext(const std::array< T, 3 > &v1, const std::array< T, 3 > &v2, std::array< T, 3 > &out)
Definition: Utils.h:170
Utils::matRot
void matRot(const std::array< T, 3 > &rot, const std::array< T, 3 > &v1, const std::array< T, 3 > &cRot, std::array< T, 3 > &out)
Definition: Utils.h:613