esphere.cuh
5.33 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
150
151
152
153
154
155
156
157
158
159
160
161
162
#ifndef RTS_ESPHERE
#define RTS_ESPHERE
#include "../math/complex.h"
#include "../math/bessel.h"
#include "../visualization/colormap.h"
#include "../optics/planewave.h"
#include "../cuda/devices.h"
#include "../optics/efield.cuh"
namespace stim{
/* This class implements a discrete representation of an electromagnetic field
in 2D scattered by a sphere. This class implements Mie scattering.
*/
template<typename P>
class esphere : public efield<P>
{
private:
stim::complex<P> n; //sphere refractive index
P a; //sphere radius
//parameters dependent on wavelength
unsigned int Nl; //number of orders for the calculation
rts::complex<P>* A; //internal scattering coefficients
rts::complex<P>* B; //external scattering coefficients
void calcNl(P kmag)
{
//return ceil( ((P)6.282 * a) / lambda + 4 * pow( ((P)6.282 * a) / lambda, ((P)1/(P)3)) + 2);
Nl = ceil( kmag*a + 4 * pow(kmag * a, (P)1/(P)3) + 2);
}
void calcAB(P k, unsigned int Nl, rts::complex<P>* A, rts::complex<P>* B)
{
/* These calculations require double precision, so they are computed
using doubles and converted to P at the end.
Input:
k = magnitude of the k vector (tau/lambda)
Nl = highest order coefficient ([0 Nl] are computed)
*/
//clear the previous coefficients
rts::complex<P>* cpuA = (rts::complex<P>*)malloc(sizeof(rts::complex<P>) * (Nl+1));
rts::complex<P>* cpuB = (rts::complex<P>*)malloc(sizeof(rts::complex<P>) * (Nl+1));
//convert to an std complex value
complex<double> nc = (rts::complex<double>)n;
//compute the magnitude of the k vector
complex<double> kna = nc * k * (double)a;
//compute the arguments k*a and k*n*a
complex<double> ka(k * a, 0.0);
//allocate space for the Bessel functions of the first and second kind (and derivatives)
unsigned int bytes = sizeof(complex<double>) * (Nl + 1);
complex<double>* cjv_ka = (complex<double>*)malloc(bytes);
complex<double>* cyv_ka = (complex<double>*)malloc(bytes);
complex<double>* cjvp_ka = (complex<double>*)malloc(bytes);
complex<double>* cyvp_ka = (complex<double>*)malloc(bytes);
complex<double>* cjv_kna = (complex<double>*)malloc(bytes);
complex<double>* cyv_kna = (complex<double>*)malloc(bytes);
complex<double>* cjvp_kna = (complex<double>*)malloc(bytes);
complex<double>* cyvp_kna = (complex<double>*)malloc(bytes);
//allocate space for the spherical Hankel functions and derivative
complex<double>* chv_ka = (complex<double>*)malloc(bytes);
complex<double>* chvp_ka = (complex<double>*)malloc(bytes);
//compute the bessel functions using the CPU-based algorithm
double vm;
cbessjyva_sph(Nl, ka, vm, cjv_ka, cyv_ka, cjvp_ka, cyvp_ka);
cbessjyva_sph(Nl, kna, vm, cjv_kna, cyv_kna, cjvp_kna, cyvp_kna);
//compute A for each order
complex<double> i(0, 1);
complex<double> a, b, c, d;
complex<double> An, Bn;
for(int l=0; l<=Nl; l++)
{
//compute the Hankel functions from j and y
chv_ka[l] = cjv_ka[l] + i * cyv_ka[l];
chvp_ka[l] = cjvp_ka[l] + i * cyvp_ka[l];
//Compute A (internal scattering coefficient)
//compute the numerator and denominator for A
a = cjv_ka[l] * chvp_ka[l] - cjvp_ka[l] * chv_ka[l];
b = cjv_kna[l] * chvp_ka[l] - chv_ka[l] * cjvp_kna[l] * nc;
//calculate A and add it to the list
rts::complex<double> An = (2.0 * l + 1.0) * pow(i, l) * (a / b);
cpuA[l] = (rts::complex<P>)An;
//Compute B (external scattering coefficient)
c = cjv_ka[l] * cjvp_kna[l] * nc - cjv_kna[l] * cjvp_ka[l];
d = cjv_kna[l] * chvp_ka[l] - chv_ka[l] * cjvp_kna[l] * nc;
//calculate B and add it to the list
rts::complex<double> Bn = (2.0 * l + 1.0) * pow(i, l) * (c / d);
cpuB[l] = (rts::complex<P>)Bn;
std::cout<<"A: "<<cpuA[l]<<" B: "<<cpuB[l]<<std::endl;
}
if(A != NULL) cudaFree(A); //free any previous coefficients
if(B != NULL) cudaFree(B);
cudaMalloc(&A, sizeof(rts::complex<P>) * (Nl+1)); //allocate memory for new coefficients
cudaMalloc(&B, sizeof(rts::complex<P>) * (Nl+1));
//copy the calculations from the CPU to the GPU
cudaMemcpy(A, cpuA, sizeof(rts::complex<P>) * (Nl+1), cudaMemcpyDeviceToHost);
cudaMemcpy(B, cpuB, sizeof(rts::complex<P>) * (Nl+1), cudaMemcpyDeviceToHost);
}
public:
esphere(unsigned int res0, unsigned int res1, P _a, rts::complex<P> _n, bool _scalar = false) :
efield<P>(res0, res1, _scalar)
{
std::cout<<"Sphere scattered field created."<<std::endl;
n = _n; //save the refractive index
a = _a; //save the radius
}
//assignment operator: build an electric field from a plane wave
efield<P> & operator= (const planewave<P> & rhs)
{
calcNl(rhs.kmag()); //compute the number of scattering coefficients
std::cout<<"Nl: "<<Nl<<std::endl;
calcAB(rhs.kmag(), Nl, A, B); //compute the scattering coefficients
//determine important parameters for the scattering domain
unsigned int sR = ceil(sqrt( (P)(pow((P)esphere::R[0],2) + pow((P)esphere::R[1],2))) );
//unsigned int thetaR = 256;
/////////////////////continue scattering code here/////////////////////////
esphere::clear(); //clear any previous field data
from_planewave(rhs); //create a field from the planewave
return *this; //return the current object
}
string str()
{
stringstream ss;
ss<<"Mie Scattered Field"<<std::endl;
ss<<(*this).efield<P>::str()<<std::endl;
ss<<"a = "<<a<<std::endl;
ss<<"n = "<<n;
return ss.str();
}
};
} //end namespace rts
#endif