| 1 | __kernel void CapCylinderKernel(__global const real *qx, __global const real *qy, __global real *_ptvalue, __global real *vol_i, |
|---|
| 2 | const real rad_cyl, const real rad_cap, const real length, const real thet, const real ph, const real sub, |
|---|
| 3 | const real scale, const real phi_weight, const real theta_float, const real rad_cap_weight, const real rad_cyl_weight, |
|---|
| 4 | const real length_weight, const int total, const int size, __const real Gauss76Wt, __const real Gauss76Z) |
|---|
| 5 | //ph is phi, sub is sldc-slds, thet is theta |
|---|
| 6 | { |
|---|
| 7 | int i = get_global_id(0); |
|---|
| 8 | if(i < total) |
|---|
| 9 | { |
|---|
| 10 | real q = sqrt(qx[i]*qx[i] + qy[i]*qy[i]); |
|---|
| 11 | real pi = 4.0*atan(1.0); |
|---|
| 12 | real theta = thet*pi/180.0; |
|---|
| 13 | real phi = ph*pi/180.0; |
|---|
| 14 | real cyl_x = cos(theta)*cos(phi); |
|---|
| 15 | real cyl_y = sin(theta); |
|---|
| 16 | real cos_val = cyl_x*qx[i]/q + cyl_y*qy[i]/q; |
|---|
| 17 | real alpha = acos(cos_val); |
|---|
| 18 | real yyy=0; real ans1=0; real ans2=0; real y=0; real xx=0; real ans=0; real zij=0; real be=0; real summj=0; |
|---|
| 19 | |
|---|
| 20 | real hDist = -1.0*sqrt(fabs(rad_cap*rad_cap-rad_cyl*rad_cyl)); |
|---|
| 21 | vol_i[i] = pi*rad_cyl*rad_cyl*length+2.0*pi/3.0*((rad_cap-hDist)*(rad_cap-hDist)*(2*rad_cap+hDist)); |
|---|
| 22 | real vaj = -1.0*hDist/rad_cap; |
|---|
| 23 | |
|---|
| 24 | for(int j=0;j<76;j++) //the 76 corresponds to the Gauss constants |
|---|
| 25 | { |
|---|
| 26 | zij = (Gauss76Z[j]*(1.0-vaj)+vaj+1.0)/2.0; |
|---|
| 27 | yyy = Gauss76Wt[j]*ConvLens_kernel(length,rad_cyl,rad_cap,q,zij,alpha); |
|---|
| 28 | summj += yyy; |
|---|
| 29 | } |
|---|
| 30 | real inner = (1.0-vaj)/2.0*summj*4.0*pi*rad_cap*rad_cap*rad_cap; |
|---|
| 31 | real arg1 = q*length/2.0*cos(alpha); |
|---|
| 32 | real arg2 = q*rad_cyl*sin(alpha); |
|---|
| 33 | yyy = inner; |
|---|
| 34 | |
|---|
| 35 | if(arg2 == 0) {be = 0.5;} |
|---|
| 36 | else { |
|---|
| 37 | be = NR_BessJ1(arg2)/arg2; |
|---|
| 38 | } |
|---|
| 39 | |
|---|
| 40 | if(arg1 == 0.0) { |
|---|
| 41 | yyy += pi*rad_cyl*rad_cyl*length*2.0*be; |
|---|
| 42 | } |
|---|
| 43 | else { |
|---|
| 44 | yyy += pi*rad_cyl*rad_cyl*length*sin(arg1)/arg1*2.0*be; |
|---|
| 45 | } |
|---|
| 46 | real answer=yyy*yyy*1.0e8*sub*sub*scale/pi*rad_cyl*rad_cyl*length+2.0*pi*(2.0*rad_cap*rad_cap*rad_cap/3.0+rad_cap*rad_cap*hDist-hDist*hDist*hDist/3.0); |
|---|
| 47 | answer/=sin(alpha); |
|---|
| 48 | |
|---|
| 49 | _ptvalue[i] = rad_cyl_weight*length_weight*rad_cap_weight*theta_weight*phi_weight*vol_i[i]*answer |
|---|
| 50 | if (_ptvalue[i] == INFINITY || _ptvalue == NAN){ |
|---|
| 51 | _ptvalue[i] = 0.0; |
|---|
| 52 | } |
|---|
| 53 | // if (size>1) { |
|---|
| 54 | // _ptvalue[i] *= fabs(cos(thet*pi/180.0)); |
|---|
| 55 | // } |
|---|
| 56 | } |
|---|
| 57 | } |
|---|
