b28f312d
 David Mayerich
added code for ev...
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
|
#ifndef TIRA_BEAM
#define TIRA_BEAM
#include "../math/vec3.h"
#include "../math/function.h"
#include "../optics/planewave.h"
#include <vector>
#include <numbers>
#include <random>
namespace tira{
namespace optics{
template<typename T>
class beam {
vec3<T> m_direction;
vec3<T> m_focus;
cvec3<T> m_E;
T m_k;
public:
/// <summary>
/// Beam constructor
/// </summary>
/// <param name="d">Propagation direction of the beam (automatically normalized).</param>
/// <param name="f">Focal point for the beam.</param>
/// <param name="pol">Beam polarization (automatically normalized and orthogonalized). Only works for linear polarization, though.</param>
/// <param name="E">Complex beam amplitude (applied to the polarization to create a linearly polarized beam).</param>
/// <param name="lambda">Wavelength</param>
beam(
vec3<T> d = vec3<T>(0, 0, 1),
vec3<T> f = vec3<T>(0, 0, 0),
vec3<T> pol = vec3<T>(1, 0, 0),
std::complex<T> E = std::complex<T>(1.0, 0.0),
T lambda = 1.0
) {
m_direction = d.direction(); //normalize and set the plane wave direction
pol = (pol - d.dot(pol)).direction(); //orthogonalize and normalize the polarization vector
m_focus = f; //set the focal point
m_E[0] = pol[0] * E; //convert the polarization and complex amplitude to a linearly-polarized E vector
m_E[1] = pol[1] * E;
m_E[2] = pol[2] * E;
m_k = 2.0 * std::numbers::pi * lambda; //calculate and store the wavenumber
}
/// <summary>
/// Generate a focal spot for the beam using Monte-Carlo sampling.
/// </summary>
/// <param name="NA">Numerical aperture of the focusing optics.</param>
/// <param name="samples">Number of Monte-Carlo samples.</param>
/// <param name="NA_in">NA of the internal obscuration.</param>
/// <returns></returns>
std::vector< planewave<T> > mcfocus(T NA, size_t samples, T NA_in = 0) {
//create a rotation matrix to orient the beam along the specified direction vector
tira::matrix_sq<double, 3> R;// = tira::matrix_sq<double, 3>::identity(); //initialize the rotation matrix to the identity matrix
tira::vec3<double> Z(0, 0, 1); //create a vector along the Z direction
tira::quaternion<double> q;
q.CreateRotation(Z, m_direction);
R = q.toMatrix3();
std::vector< planewave<T> > result(samples); //generate an array of N plane waves
double alpha = std::asin(NA); //calculate the angle of the band-pass
double alpha_in = std::asin(NA_in); //calculate the angle of the central obscuration
double cos_alpha = std::cos(alpha); //calculate the cosine of the band-pass and central obscuration
double cos_alpha_in = std::cos(alpha_in);
//random sampling is done using Archimede's method - uniform samples are taken in cylindrical space and projected onto a sphere
std::random_device rd; //create a random number generator
std::default_random_engine eng(rd());
std::uniform_real_distribution<> theta_dist(0.0, 2.0 * std::numbers::pi); //create a distribution for theta (in cylindrical coordinates)
std::uniform_real_distribution<> z_dist(cos_alpha, cos_alpha_in); //create a uniform distribution for sampling the z-axis
//generate a guide plane wave
planewave<T> pw(m_direction[0] * m_k, m_direction[1] * m_k, m_direction[2] * m_k, m_E[0], m_E[1], m_E[2]);
double theta, phi; //variables to store the coordinates for theta/phi on a sphere
tira::vec3<T> dir;
T w = 1.0 / samples; //integration weight
for (size_t n = 0; n < samples; n++) {
theta = theta_dist(eng); //randomly generate a theta coordinate between 0 and 2pi
phi = std::acos(z_dist(eng)); //randomly generate a z coordinate based on the NA and project to phi
dir = R * tira::vec3<double>(1.0, theta, phi).sph2cart(); //convert the value to Cartesian coordinates and store the sample vector
result[n] = pw.refract(dir); //refract the plane wave to align with the sample direction
result[n] = result[n].translate(m_focus[0], m_focus[1], m_focus[2]); //refocus the plane wave to the focal position
result[n] = result[n] * w; //apply the integration weight (1/N)
}
return result; //return the sample array
}
};
} //end namespace optics
} //end namespace tira
#endif
|
