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
|
#include <vector>
|
8a86bd56
David Mayerich
changed rts names...
|
9
|
namespace stim{
|
a9275be5
David Mayerich
added vector fiel...
|
10
11
|
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
|
{
public:
enum beam_type {Uniform, Bartlett, Hamming, Hanning};
private:
|
559d0fcb
David Mayerich
wave interactions...
|
19
|
P _na[2]; //numerical aperature of the focusing optics
|
a9275be5
David Mayerich
added vector fiel...
|
20
|
vec<P> f; //focal point
|
a9275be5
David Mayerich
added vector fiel...
|
21
|
function<P, P> apod; //apodization function
|
559d0fcb
David Mayerich
wave interactions...
|
22
|
unsigned int apod_res; //resolution of apodization filter functions
|
a9275be5
David Mayerich
added vector fiel...
|
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
|
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
|
beam(
|
559d0fcb
David Mayerich
wave interactions...
|
72
|
vec<P> k = rts::vec<P>(0, 0, rtsTAU),
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
73
74
|
vec<P> _E0 = rts::vec<P>(1, 0, 0),
beam_type _apod = Uniform)
|
559d0fcb
David Mayerich
wave interactions...
|
75
|
: planewave<P>(k, _E0)
|
a9275be5
David Mayerich
added vector fiel...
|
76
|
{
|
559d0fcb
David Mayerich
wave interactions...
|
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
|
}
|
559d0fcb
David Mayerich
wave interactions...
|
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
|
beam<P> refract(rts::vec<P> kn) const{
beam<P> new_beam;
new_beam._na[0] = _na[0];
new_beam._na[1] = _na[1];
rts::planewave<P> pw = planewave<P>::bend(kn);
//std::cout<<pw.str()<<std::endl;
new_beam.k = pw.kvec();
new_beam.E0 = pw.E();
return new_beam;
}
|
a9275be5
David Mayerich
added vector fiel...
|
100
|
///Numerical Aperature functions
|
559d0fcb
David Mayerich
wave interactions...
|
101
|
void NA(P na)
|
a9275be5
David Mayerich
added vector fiel...
|
102
|
{
|
559d0fcb
David Mayerich
wave interactions...
|
103
104
|
_na[0] = (P)0;
_na[1] = na;
|
a9275be5
David Mayerich
added vector fiel...
|
105
|
}
|
559d0fcb
David Mayerich
wave interactions...
|
106
|
void NA(P na0, P na1)
|
a9275be5
David Mayerich
added vector fiel...
|
107
|
{
|
559d0fcb
David Mayerich
wave interactions...
|
108
109
|
_na[0] = na0;
_na[1] = na1;
|
a9275be5
David Mayerich
added vector fiel...
|
110
111
|
}
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
112
113
114
115
116
117
118
|
/*string str() :
{
stringstream ss;
ss<<"Beam Center: "<<k<<std::endl;
return ss.str();
}*/
|
a9275be5
David Mayerich
added vector fiel...
|
119
120
|
//Monte-Carlo decomposition into plane waves
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
121
|
std::vector< planewave<P> > mc(unsigned int N = 100000, unsigned int seed = 0) const
|
a9275be5
David Mayerich
added vector fiel...
|
122
123
124
125
|
{
/*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...
|
126
|
seed = seed for the random number generator
|
a9275be5
David Mayerich
added vector fiel...
|
127
|
*/
|
a9275be5
David Mayerich
added vector fiel...
|
128
129
|
srand(seed); //seed the random number generator
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
130
131
|
vec<P> k_hat = beam::k.norm();
|
a9275be5
David Mayerich
added vector fiel...
|
132
|
///compute the rotation operator to transform (0, 0, 1) to k
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
133
|
P cos_angle = k_hat.dot(rts::vec<P>(0, 0, 1));
|
a9275be5
David Mayerich
added vector fiel...
|
134
|
rts::matrix<P, 3> rotation;
|
559d0fcb
David Mayerich
wave interactions...
|
135
136
137
138
139
140
|
//if the cosine of the angle is -1, the rotation is just a flip across the z axis
if(cos_angle == -1){
rotation(2, 2) = -1;
}
else if(cos_angle != 1.0)
|
a9275be5
David Mayerich
added vector fiel...
|
141
|
{
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
142
|
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...
|
143
144
145
146
147
148
149
150
|
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];
|
559d0fcb
David Mayerich
wave interactions...
|
151
152
|
PHI[0] = (P)asin(_na[0]);
PHI[1] = (P)asin(_na[1]);
|
a9275be5
David Mayerich
added vector fiel...
|
153
154
155
156
157
158
159
160
161
|
//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...
|
162
163
164
165
166
167
168
169
170
171
172
173
174
175
|
//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
|
559d0fcb
David Mayerich
wave interactions...
|
176
177
|
//std::cout<<"k prime: "<<rotation<<std::endl;
samples.push_back(planewave<P>::refract(k_prime) * apod(phi/PHI[1]));
|
a9275be5
David Mayerich
added vector fiel...
|
178
179
180
181
182
|
}
return samples;
}
|
559d0fcb
David Mayerich
wave interactions...
|
183
184
185
186
|
std::string str()
{
std::stringstream ss;
ss<<"Beam:"<<std::endl;
|
ecfd14df
David Mayerich
implemented D fie...
|
187
188
|
//ss<<" Central Plane Wave: "<<beam::E0<<" e^i ( "<<beam::k<<" . r )"<<std::endl;
ss<<" Central Plane Wave: "<<beam::k<<std::endl;
|
559d0fcb
David Mayerich
wave interactions...
|
189
190
191
192
193
194
195
196
197
198
|
if(_na[0] == 0)
ss<<" NA: "<<_na[1];
else
ss<<" NA: "<<_na[0]<<" -- "<<_na[1];
return ss.str();
}
|
a9275be5
David Mayerich
added vector fiel...
|
199
200
201
202
|
};
}
|
8a86bd56
David Mayerich
changed rts names...
|
203
|
#endif
|