821513c8
 David Mayerich
ERROR plane wave ...
|
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
|
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
|
