__device__ complexf besselh11Func complexf complexf cii complexf cone

 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
__device__ complexf besselh11Func(complexf z, complexf *cii, complexf *cone, complexf *czero)
{
complexf z1, z2, cr, cp, cs, cp0, cq0, cp1, cq1, ct1, ct2, cu;
complexf cj1, cy1;
float a0, w1;
int k, kz;
/*complexf sharedCii(0.0f, 1.0f);
complexf sharedCone(1.0f, 0.0f);
complexf sharedCzero(0.0f, 0.0f);*/
/*__shared__ complexf sharedCii;
__shared__ complexf sharedCone;
__shared__ complexf sharedCzero;
sharedCii = (*cii);
sharedCone = (*cone);
sharedCzero = (*czero);
__syncthreads();*/
a0 = abs(z);
z2 = z*z;
z1 = z;
if (a0 == 0.0f) {
cj1 = (*czero);
cy1 = complexf(-1e38f, 0.0f);
return 0.0f;
}
if (real(z) < 0.0f) z1 = -z;
if (a0 <= 12.0f) {
cj1 = (*cone);
cr = (*cone);
for (k = 1; k <= 40; k++) {
cr *= -0.25f* z2 / (k*(k + 1.0f));
cj1 += cr;
if (abs(cr) < abs(cj1)*eps) break;
}
cj1 *= 0.5f*z1;
w1 = 0.0f;
cr = (*cone);
cs = (*cone);
for (k = 1; k <= 40; k++) {
w1 += 1.0 / k;
cr *= -0.25f*z2 / (k*(k + 1.0f));
cp = cr * (2.0f*w1 + 1.0f / (k + 1.0f));
cs += cp;
if (abs(cp) < abs(cs)*eps) break;
}
auto tm21 = z1; // log(0.5f * z1);
auto tmp1 = tm21; // + el;
auto tmp2 = cj1 - 1.0f / z1 - 0.25f * z1 * cs;
cy1 = M_2_PI * tmp1 * tmp2;
}
else {
if (a0 >= 50.0f) kz = 8; // can be changed to 10
else if (a0 >= 35.0f) kz = 10; // " " " 12
else kz = 12; // " " " 14
ct1 = z1 - M_PI_4;
cp0 = (*cone);
for (k = 0; k<kz; k++) {
cp0 += a[k] * pow(z1, -2.0f*k - 2.0f);
}
cq0 = -0.125f / z1;
for (k = 0; k<kz; k++) {
cq0 += b[k] * pow(z1, -2.0f*k - 3.0f);
}
cu = sqrt(M_2_PI / z1);
ct2 = z1 - 0.75f*M_PI;
cp1 = (*cone);
for (k = 0; k<kz; k++) {
cp1 += a1[k] * pow(z1, -2.0f*k - 2.0f);
}
cq1 = 0.375f / z1;
for (k = 0; k<kz; k++) {
cq1 += b1[k] * pow(z1, -2.0f*k - 3.0f);
}
cj1 = cu*(cp1*cos(ct2) - cq1*sin(ct2));
cy1 = cu*(cp1*sin(ct2) + cq1*cos(ct2));
}
if (real(z) < 0.0f) {
if (imag(z) < 0.0f) {
cy1 = -(cy1 - 2.0f*(*cii)*cj1);
}
else if (imag(z) > 0.0f) {
cy1 = -(cy1 + 2.0f*(*cii)*cj1);
}
cj1 = -cj1;
}
return (cj1 + (*cii) * cy1);
}