a9275be5
David Mayerich
added vector fiel...
|
1
2
3
4
5
|
#ifndef RTS_BEAM
#define RTS_BEAM
#include "../math/vector.h"
#include "../math/function.h"
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
6
|
#include "../optics/planewave.h"
|
a9275be5
David Mayerich
added vector fiel...
|
7
8
9
10
11
|
#include <vector>
namespace rts{
template<typename P>
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
12
|
class beam : public planewave<P>
|
a9275be5
David Mayerich
added vector fiel...
|
13
14
15
16
17
18
19
20
|
{
public:
enum beam_type {Uniform, Bartlett, Hamming, Hanning};
private:
P na[2]; //numerical aperature of the focusing optics
vec<P> f; //focal point
|
a9275be5
David Mayerich
added vector fiel...
|
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
|
function<P, P> apod; //apodization function
unsigned int apod_res; //resolution of complex apodization filters
void apod_uniform()
{
apod = (P)1;
}
void apod_bartlett()
{
apod = (P)1;
apod.insert((P)1, (P)0);
}
void apod_hanning()
{
apod = (P)0;
P x, y;
for(unsigned int n=0; n<apod_res; n++)
{
x = (P)n/(P)apod_res;
y = pow( cos( ((P)3.14159 * x) / 2 ), 2);
apod.insert(x, y);
}
}
void apod_hamming()
{
apod = (P)0;
P x, y;
for(unsigned int n=0; n<apod_res; n++)
{
x = (P)n/(P)apod_res;
y = (P)27/(P)50 + ( (P)23/(P)50 ) * cos((P)3.14159 * x);
apod.insert(x, y);
}
}
void set_apod(beam_type type)
{
if(type == Uniform)
apod_uniform();
if(type == Bartlett)
apod_bartlett();
if(type == Hanning)
apod_hanning();
if(type == Hamming)
apod_hamming();
}
public:
///constructor: build a default beam (NA=1.0)
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
71
72
73
74
75
|
beam(
vec<P> _k = rts::vec<P>(0, 0, TAU),
vec<P> _E0 = rts::vec<P>(1, 0, 0),
beam_type _apod = Uniform)
: planewave<P>(_k, _E0)
|
a9275be5
David Mayerich
added vector fiel...
|
76
77
78
|
{
na[0] = (P)0.0;
na[1] = (P)1.0;
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
79
|
f = vec<P>( (P)0, (P)0, (P)0 );
|
a9275be5
David Mayerich
added vector fiel...
|
80
81
|
apod_res = 256; //set the default resolution for apodization filters
set_apod(_apod); //set the apodization function type
|
a9275be5
David Mayerich
added vector fiel...
|
82
83
84
85
86
87
88
89
90
91
92
93
94
95
|
}
///Numerical Aperature functions
void NA(P _na)
{
na[0] = (P)0;
na[1] = _na;
}
void NA(P _na0, P _na1)
{
na[0] = _na0;
na[1] = _na1;
}
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
96
97
98
99
100
101
102
|
/*string str() :
{
stringstream ss;
ss<<"Beam Center: "<<k<<std::endl;
return ss.str();
}*/
|
a9275be5
David Mayerich
added vector fiel...
|
103
104
|
//Monte-Carlo decomposition into plane waves
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
105
|
std::vector< planewave<P> > mc(unsigned int N = 100000, unsigned int seed = 0) const
|
a9275be5
David Mayerich
added vector fiel...
|
106
107
108
109
|
{
/*Create Monte-Carlo samples of a cassegrain objective by performing uniform sampling
of a sphere and projecting these samples onto an inscribed sphere.
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
110
|
seed = seed for the random number generator
|
a9275be5
David Mayerich
added vector fiel...
|
111
|
*/
|
a9275be5
David Mayerich
added vector fiel...
|
112
113
|
srand(seed); //seed the random number generator
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
114
115
|
vec<P> k_hat = beam::k.norm();
|
a9275be5
David Mayerich
added vector fiel...
|
116
|
///compute the rotation operator to transform (0, 0, 1) to k
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
117
|
P cos_angle = k_hat.dot(rts::vec<P>(0, 0, 1));
|
a9275be5
David Mayerich
added vector fiel...
|
118
119
120
|
rts::matrix<P, 3> rotation;
if(cos_angle != 1.0)
{
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
121
|
rts::vec<P> r_axis = rts::vec<P>(0, 0, 1).cross(k_hat).norm(); //compute the axis of rotation
|
a9275be5
David Mayerich
added vector fiel...
|
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
|
P angle = acos(cos_angle); //compute the angle of rotation
rts::quaternion<P> quat; //create a quaternion describing the rotation
quat.CreateRotation(angle, r_axis);
rotation = quat.toMatrix3(); //compute the rotation matrix
}
//find the phi values associated with the cassegrain ring
P PHI[2];
PHI[0] = (P)asin(na[0]);
PHI[1] = (P)asin(na[1]);
//calculate the z-axis cylinder coordinates associated with these angles
P Z[2];
Z[0] = cos(PHI[0]);
Z[1] = cos(PHI[1]);
P range = Z[0] - Z[1];
std::vector< planewave<P> > samples; //create a vector of plane waves
|
a9275be5
David Mayerich
added vector fiel...
|
141
142
143
144
145
146
147
148
149
150
151
152
153
154
|
//draw a distribution of random phi, z values
P z, phi, theta;
for(int i=0; i<N; i++) //for each sample
{
z = ((P)rand() / (P)RAND_MAX) * range + Z[1]; //find a random position on the surface of a cylinder
theta = ((P)rand() / (P)RAND_MAX) * 2 * (P)3.14159;
phi = acos(z); //project onto the sphere, computing phi in spherical coordinates
//compute and store cartesian coordinates
rts::vec<P> spherical(1, theta, phi); //convert from spherical to cartesian coordinates
rts::vec<P> cart = spherical.sph2cart();
vec<P> k_prime = rotation * cart; //create a sample vector
//store a wave refracted along the given direction
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
155
|
samples.push_back(beam::refract(k_prime) * apod(phi/PHI[1]));
|
a9275be5
David Mayerich
added vector fiel...
|
156
157
158
159
160
161
162
163
164
165
|
}
return samples;
}
};
}
#endif
|