Changeset 50e1e40 in sasmodels for sasmodels/models/cylinder.c

Timestamp:

Mar 1, 2016 7:49:00 PM (8 years ago)

Author:

Paul Kienzle <pkienzle@…>

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:

ad90df9

Parents:

a4a7308

Message:

use lib functions to simplify barbell, capped_cylinder, cylinder, ellipsoid, triaxial_ellipsoid; fix barbell calculation

File:

• sasmodels/models/cylinder.c (modified) (3 diffs)

sasmodels/models/cylinder.c

@@ -3 +3 @@
 double Iqxy(double qx, double qy, double sld, double solvent_sld,
     double radius, double length, double theta, double phi);
-
-// twovd = 2 * volume * delta_rho
-// besarg = q * R * sin(alpha)
-// siarg = q * L/2 * cos(alpha)
-double _cyl(double besarg, double siarg);
-double _cyl(double besarg, double siarg)
-{
-    const double bj = (besarg == 0.0 ? 0.5 : J1(besarg)/besarg);
-    const double si = (siarg == 0.0 ? 1.0 : sin(siarg)/siarg);
-    return si*bj;
-}
 
 double form_volume(double radius, double length)
@@ -26 +15 @@
     double length)
 {
+    // TODO: return NaN if radius<0 or length<0?
+    // precompute qr and qh to save time in the loop
     const double qr = q*radius;
     const double qh = q*0.5*length;
+
+    // translate a point in [-1,1] to a point in [0, pi/2]
+    const double zm = M_PI_4;
+    const double zb = M_PI_4;
+
     double total = 0.0;
-    // double lower=0, upper=M_PI_2;
     for (int i=0; i<76 ;i++) {
-        // translate a point in [-1,1] to a point in [lower,upper]
-        //const double alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
-        const double alpha = M_PI_4*(Gauss76Z[i] + 1.0);
+        const double alpha = Gauss76Z[i]*zm + zb;
         double sn, cn;
         SINCOS(alpha, sn, cn);
-        // For a bit of efficiency, we are moving the 2 V delta rho constant
-        // factor, 2Vdrho, out of the loop, so this is fq/2Vdrho rather than fq.
-        const double fq = _cyl(qr*sn, qh*cn);
-        total += Gauss76Wt[i] * fq * fq * sn;
+        const double fq = sinc(qh*cn) * J1c(qr*sn);
+        total += Gauss76Wt[i] * fq*fq * sn;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
-    //const double form = (upper-lower)/2.0*total;
-    const double twoVdrho = 2.0*(sld-solvent_sld)*form_volume(radius, length);
-    return 1.0e-4 * twoVdrho * twoVdrho * total * M_PI_4;
+    const double form = total*zm;
+    const double s = (sld - solvent_sld) * form_volume(radius, length);
+    return 1.0e-4 * s * s * form;
 }
 
@@ -56 +47 @@
     double phi)
 {
-    // TODO: check that radius<0 and length<0 give zero scattering.
-    // This should be the case since the polydispersity weight vector should
-    // be zero length, and this function never called.
+    // TODO: return NaN if radius<0 or length<0?
     double sn, cn; // slots to hold sincos function output
 
     // Compute angle alpha between q and the cylinder axis
     SINCOS(theta*M_PI_180, sn, cn);
-    const double q = sqrt(qx*qx+qy*qy);
+    const double q = sqrt(qx*qx + qy*qy);
     const double cos_val = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q);
     const double alpha = acos(cos_val);
+
     SINCOS(alpha, sn, cn);
-    //sn = sqrt(1.0 - cos_val*cos_val);
-    //sn = 1.0 - 0.5*cos_val*cos_val;  // if cos_val is very small
-    //cn = cos_val;
-
-    const double twovd = 2.0*(sld-solvent_sld)*form_volume(radius, length);
-    const double fq = twovd * _cyl(q*radius*sn, q*0.5*length*cn);
-    return 1.0e-4 * fq * fq;
+    const double fq = sinc(q*0.5*length*cn) * J1c(q*radius*sn);
+    const double s = (sld-solvent_sld) * form_volume(radius, length);
+    return 1.0e-4 * square(s * fq);
 }