Changeset 0d6e865 in sasmodels for sasmodels/models/cylinder.c

Timestamp:

Oct 11, 2016 11:42:00 AM (8 years ago)

Author:

dirk

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

1fdb555

Parents:

19e7ca7

Message:

Rewriting selected models in spherical coordinates. This fixes ticket #492 for these models.

File:

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

sasmodels/models/cylinder.c

@@ -1 +1 @@
 double form_volume(double radius, double length);
+double fq(double q, double sn, double cn,double radius, double length);
+double orient_avg_1D(double q, double radius, double length);
 double Iq(double q, double sld, double solvent_sld, double radius, double length);
 double Iqxy(double qx, double qy, double sld, double solvent_sld,
@@ -11 +13 @@
 }
 
+double fq(double q, double sn, double cn, double radius, double length)
+{
+    // precompute qr and qh to save time in the loop
+    const double qr = q*radius;
+    const double qh = q*0.5*length; 
+    return  sas_J1c(qr*sn) * sinc(qh*cn) ;
+}
+
+double orient_avg_1D(double q, double radius, double 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;
+    for (int i=0; i<76 ;i++) {
+        const double alpha = Gauss76Z[i]*zm + zb;
+        double sn, cn; // slots to hold sincos function output
+        // alpha(theta,phi) the projection of the cylinder on the detector plane
+        SINCOS(alpha, sn, cn);
+        total += Gauss76Wt[i] * fq(q, sn, cn, radius, length)* fq(q, sn, cn, radius, length) * sn;
+    }
+    // translate dx in [-1,1] to dx in [lower,upper]
+    return total*zm;
+}
+
 double Iq(double q,
     double sld,
@@ -17 +45 @@
     double length)
 {
-    // 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;
-    for (int i=0; i<76 ;i++) {
-        const double alpha = Gauss76Z[i]*zm + zb;
-        double sn, cn;
-        SINCOS(alpha, sn, cn);
-        const double fq = sinc(qh*cn) * sas_J1c(qr*sn);
-        total += Gauss76Wt[i] * fq*fq * sn;
-    }
-    // translate dx in [-1,1] to dx in [lower,upper]
-    const double form = total*zm;
     const double s = (sld - solvent_sld) * form_volume(radius, length);
-    return 1.0e-4 * s * s * form;
+    return 1.0e-4 * s * s * orient_avg_1D(q, radius, length);
 }
 
@@ -51 +61 @@
 
     // Compute angle alpha between q and the cylinder axis
-    SINCOS(theta*M_PI_180, sn, cn);
+    SINCOS(phi*M_PI_180, sn, cn);
     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 cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q);
+
     const double alpha = acos(cos_val);
 
     SINCOS(alpha, sn, cn);
-    const double fq = sinc(q*0.5*length*cn) * sas_J1c(q*radius*sn);
     const double s = (sld-solvent_sld) * form_volume(radius, length);
-    return 1.0e-4 * square(s * fq);
+    return 1.0e-4 * square(s * fq(q, sn, cn, radius, length));
 }