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