halfspace.h
2.92 KB
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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
#ifndef RTS_HALFSPACE_H
#define RTS_HALFSPACE_H
namespace rts{
template<typename T> class halfspace;
template<typename T> class efield;
}
template<typename T>
void operator<<(rts::efield<T> &ef, rts::halfspace<T> hs);
namespace rts{
template<class T>
class halfspace
{
private:
rts::plane<T> S; //surface plane splitting the half space
rts::complex<T> n0; //refractive index at the front of the plane
rts::complex<T> n1; //refractive index at the back of the plane
//lists of waves in front (pw0) and behind (pw1) the plane
std::vector< rts::planewave<T> > w0;
std::vector< rts::planewave<T> > w1;
//rts::planewave<T> pi; //incident plane wave
//rts::planewave<T> pr; //plane wave reflected from the surface
//rts::planewave<T> pt; //plane wave transmitted through the surface
void init(){
n0 = 1.0;
n1 = 1.5;
}
//compute which side of the interface is hit by the incoming plane wave (0 = front, 1 = back)
bool facing(planewave<T> p){
if(p.kvec().dot(S.norm()) > 0)
return 1;
else
return 0;
}
T calc_theta_i(vec<T> v){
vec<T> v_hat = v.norm();
//compute the cosine of the angle between k_hat and the plane normal
T cos_theta_i = v_hat.dot(S.norm());
return acos(abs(cos_theta_i));
}
T calc_theta_t(T ni_nt, T theta_i){
T sin_theta_t = ni_nt * sin(theta_i);
return asin(sin_theta_t);
}
public:
//constructors
halfspace(){
init();
}
halfspace(T na, T nb){
init();
n0 = na;
n1 = nb;
}
halfspace(T na, T nb, plane<T> p){
n0 = na;
n1 = nb;
S = p;
}
//compute the transmitted and reflective waves given the incident (vacuum) plane wave p
void incident(rts::planewave<T> p){
planewave<T> r, t;
p.scatter(S, n1.real()/n0.real(), r, t);
//std::cout<<"i+r: "<<p.r()[0] + r.r()[0]<<std::endl;
//std::cout<<"t: "<<t.r()[0]<<std::endl;
if(facing(p)){
w1.push_back(p);
if(r.E().len() != 0)
w1.push_back(r);
if(t.E().len() != 0)
w0.push_back(t);
}
else{
w0.push_back(p);
if(r.E().len() != 0)
w0.push_back(r);
if(t.E().len() != 0)
w1.push_back(t);
}
}
void incident(rts::beam<T> b, unsigned int N = 10000){
//generate a plane wave decomposition for the beam
std::vector< planewave<T> > pw_list = b.mc(N);
//calculate the reflected and refracted waves for each incident wave
for(unsigned int w = 0; w < pw_list.size(); w++){
incident(pw_list[w]);
}
}
std::string str(){
std::stringstream ss;
ss<<"Half Space-------------"<<std::endl;
ss<<"n0: "<<n0<<std::endl;
ss<<"n1: "<<n1<<std::endl;
ss<<S.str();
if(w0.size() != 0 || w1.size() != 0){
ss<<std::endl<<"Plane Waves:"<<std::endl;
for(unsigned int w = 0; w < w0.size(); w++){
ss<<"w0 = "<<w<<": "<<w0[w]<<std::endl;
}
for(unsigned int w = 0; w < w1.size(); w++){
ss<<"w1 = "<<w<<": "<<w1[w]<<std::endl;
}
}
return ss.str();
}
////////FRIENDSHIP
friend void operator<< <> (rts::efield<T> &ef, rts::halfspace<T> hs);
};
}
#endif