Changeset 8a20be5 in sasmodels for Kernel-Cylinder.cpp

Timestamp:

Jul 10, 2014 3:05:08 PM (10 years ago)

Author:

HMP1 <helen.park@…>

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

8faffcd

Parents:

5378e40

Message:

Added a fit2 (fits two different models at different angles)
(preliminary) Added CoreshellCyl? and CapCyl? Kernels
(preliminary) Updated kernels to include functions

File:

• Kernel-Cylinder.cpp (modified) (3 diffs)

Kernel-Cylinder.cpp

@@ -1 @@
-__kernel void CylinderKernel(__global const float *qx, global const float *qy, __global float *_ptvalue, const float sub,
-const float rr, const float h, const float scale, const float radius_weight, const float length_weight,
-const float theta_weight, const float phi_weight, const float cyl_theta,
-const float cyl_phi, const int count, const int size)
+__kernel void CylinderKernel(__global const real *qx, global const real *qy, __global real *_ptvalue, const real sub,
+const real rr, const real h, const real scale, const real radius_weight, const real length_weight,
+const real theta_weight, const real phi_weight, const real cyl_theta,
+const real cyl_phi, const int count, const int size)
 {
         // qq is the q-value for the calculation (1/A)
@@ -11 +11 @@
     if(i < count)
     {
-        float qq = sqrt(qx[i]*qx[i]+qy[i]*qy[i]);
+        real qq = sqrt(qx[i]*qx[i]+qy[i]*qy[i]);
 
-        float pi = 4.0*atan(1.0);
-        float theta = cyl_theta*pi/180.0;
-        float phi = cyl_phi*pi/180.0;
+        real pi = 4.0*atan(1.0);
+        real theta = cyl_theta*pi/180.0;
+        real phi = cyl_phi*pi/180.0;
 
-        float cyl_x = cos(theta)*cos(phi);
-        float cyl_y = sin(theta);
-        float cos_val = cyl_x*(qx[i]/qq) + cyl_y*(qy[i]/qq);
+        real cyl_x = cos(theta)*cos(phi);
+        real cyl_y = sin(theta);
+        real cos_val = cyl_x*(qx[i]/qq) + cyl_y*(qy[i]/qq);
 
-        float alpha = acos(cos_val);
+        real alpha = acos(cos_val);
         if(alpha == 0.0){
             alpha = 1.0e-26;
         }
-        float besarg = qq*rr*sin(alpha);
-        float siarg = qq*h/2*cos(alpha);
+        real besarg = qq*rr*sin(alpha);
+        real siarg = qq*h/2*cos(alpha);
+        real be=0.0; real si=0.0;
 
-            float xx=0.0; float y=0.0; float bj=0.0; float ans1=0.0; float ans2=0.0; float z=0.0; float answer=0.0;
-            float contrast=0.0; float form=0.0; float be=0.0; float si=0.0;
+        real bj = NR_BessJ1(besarg);
 
-        float ax = fabs(besarg);
+        real d1 = qq*rr*sin(alpha);
 
-        if(ax < 8.0) {
-            y=besarg*besarg;
-            ans1=besarg*(72362614232.0+y*(-7895059235.0+y*(242396853.1+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
-            ans2=144725228442.0+y*(2300535178.0+y*(18583304.74+y*(99447.43394+y*(376.9991397+y*1.0))));
-            bj=ans1/ans2;
+        if (besarg == 0.0){
+            be = sin(alpha);
         }
         else{
-            z=8.0/ax;
-            y=z*z;
-            xx=ax - 2.356194491;
-            ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4+y*(0.2457520174e-5+y*(-0.240337019e-6))));
-            ans2=0.04687499995+y*(-0.2002690873e-3+y*(0.8449199096e-5+y*(-0.88228987e-6+y*0.105787412e-6)));
-            bj=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
-
-            if(besarg < 0.0){bj*=-1;}
+            be = bj*bj*4.0*sin(alpha)/(d1*d1);
+        }
+        if(siarg == 0.0){
+            si = 1.0;
+        }
+        else{
+            si = sin(siarg)*sin(siarg)/(siarg*siarg);
         }
 
-        float d1 = qq*rr*sin(alpha);
-
-        if (besarg == 0.0) {be = sin(alpha);}
-        else {be = bj*bj*4.0*sin(alpha)/(d1*d1);}
-        if(siarg == 0.0) {si = 1.0;}
-        else{si = sin(siarg)*sin(siarg)/(siarg*siarg);}
-
-        form = be*si/sin(alpha);
-        answer = sub*sub*form*acos(-1.0)*rr*rr*h*1.0e8*scale;
+        real form = be*si/sin(alpha);
+        real answer = sub*sub*form*acos(-1.0)*rr*rr*h*1.0e8*scale;
 
         _ptvalue[i] = radius_weight*length_weight*theta_weight*phi_weight*answer*pow(rr,2)*h;
@@ -64 +53 @@
             _ptvalue[i] *= fabs(cos(cyl_theta*pi/180.0));
         }
-}
+    }
 }