Changeset 2a0b2b1 in sasmodels for sasmodels/models/cylinder.c

Timestamp:

Apr 14, 2017 8:30:29 AM (7 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

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

Children:

fdd56a1

Parents:

9901384

Message:

restructure all 2D models to work with (qa,qb,qc) = rotate(qx,qy) rather than working with angles directly in preparation for revised jitter algorithm

File:

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

sasmodels/models/cylinder.c

@@ -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,
-    double radius, double length, double theta, double phi);
-
 #define INVALID(v) (v.radius<0 || v.length<0)
 
-double form_volume(double radius, double length)
+static double
+form_volume(double radius, double length)
 {
     return M_PI*radius*radius*length;
 }
 
-double fq(double q, double sn, double cn, double radius, double length)
+static double
+fq(double qab, double qc, 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_2J1x_x(qr*sn) * sas_sinx_x(qh*cn);
+    return sas_2J1x_x(qab*radius) * sas_sinx_x(qc*0.5*length);
 }
 
-double orient_avg_1D(double q, double radius, double length)
+static 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; 
+    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] * square( fq(q, sn, cn, radius, length) ) * sn;
+        const double theta = Gauss76Z[i]*zm + zb;
+        double sin_theta, cos_theta; // slots to hold sincos function output
+        // theta (theta,phi) the projection of the cylinder on the detector plane
+        SINCOS(theta , sin_theta, cos_theta);
+        const double form = fq(q*sin_theta, q*cos_theta, radius, length);
+        total += Gauss76Wt[i] * form * form * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -39 +33 @@
 }
 
-double Iq(double q,
+static double
+Iq(double q,
     double sld,
     double solvent_sld,
@@ -49 +44 @@
 }
 
-
-double Iqxy(double qx, double qy,
+static double
+Iqxy(double qx, double qy,
     double sld,
     double solvent_sld,
@@ -60 +55 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
-    //printf("sn: %g cn: %g\n", sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
     const double s = (sld-solvent_sld) * form_volume(radius, length);
-    const double form = fq(q, sin_alpha, cos_alpha, radius, length);
+    const double form = fq(qab, qc, radius, length);
     return 1.0e-4 * square(s * form);
 }