sphere.cpp
3.35 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
#include "sphere.h"
#include "rts/math/complex.h"
#include <complex>
#include <stdlib.h>
using namespace rts;
using namespace std;
int cbessjyva(double v,complex<double> z,double &vm,complex<double>*cjv,
complex<double>*cyv,complex<double>*cjvp,complex<double>*cyvp);
int cbessjyva_sph(int v,complex<double> z,double &vm,complex<double>*cjv,
complex<double>*cyv,complex<double>*cjvp,complex<double>*cyvp);
void sphere::calcCoeff(ptype lambda, rtsComplex<ptype> ri)
{
/* These calculations are done at high-precision on the CPU
since they are only required once for each wavelength.
Input:
lambda = wavelength of the incident field
n = complex refractive index of the sphere
*/
//clear the previous coefficients
A.clear();
B.clear();
//convert to an std complex value
complex<double> nc(ri.real(), ri.imag());
n = ri;
//compute the magnitude of the k vector
double k = 2 * PI / lambda;
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)
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);
//cout<<"Begin Sphere---------"<<endl;
//cout<<"Nl = "<<Nl<<endl;
//cout<<"ka = "<<ka<<endl;
//cout<<"kna = "<<kna<<endl;
//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
An = (2.0 * l + 1.0) * pow(i, l) * (a / b);
A.push_back(bsComplex(An.real(), An.imag()));
//cout<<"A: "<<An<<endl;
//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
Bn = (2.0 * l + 1.0) * pow(i, l) * (c / d);
B.push_back(bsComplex(Bn.real(), Bn.imag()));
//cout<<"B: "<<Bn<<endl;
}
}