Changeset 2a0b2b1 in sasmodels for sasmodels/models/capped_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/capped_cylinder.c (modified) (5 diffs)

sasmodels/models/capped_cylinder.c

@@ -14 +14 @@
 //   radius_cap is the radius of the lens
 //   length is the cylinder length, or the separation between the lens halves
-//   alpha is the angle of the cylinder wrt q.
+//   theta is the angle of the cylinder wrt q.
 static double
-_cap_kernel(double q, double h, double radius_cap,
-                      double half_length, double sin_alpha, double cos_alpha)
+_cap_kernel(double qab, double qc, double h, double radius_cap,
+            double half_length)
 {
     // translate a point in [-1,1] to a point in [lower,upper]
@@ -26 +26 @@
 
     // cos term in integral is:
-    //    cos (q (R t - h + L/2) cos(alpha))
+    //    cos (q (R t - h + L/2) cos(theta))
     // so turn it into:
     //    cos (m t + b)
     // where:
-    //    m = q R cos(alpha)
-    //    b = q(L/2-h) cos(alpha)
-    const double m = q*radius_cap*cos_alpha; // cos argument slope
-    const double b = q*(half_length-h)*cos_alpha; // cos argument intercept
-    const double qrst = q*radius_cap*sin_alpha; // Q*R*sin(theta)
+    //    m = q R cos(theta)
+    //    b = q(L/2-h) cos(theta)
+    const double m = radius_cap*qc; // cos argument slope
+    const double b = (half_length-h)*qc; // cos argument intercept
+    const double qab_r = radius_cap*qab; // Q*R*sin(theta)
     double total = 0.0;
     for (int i=0; i<76 ;i++) {
         const double t = Gauss76Z[i]*zm + zb;
         const double radical = 1.0 - t*t;
-        const double bj = sas_2J1x_x(qrst*sqrt(radical));
+        const double bj = sas_2J1x_x(qab_r*sqrt(radical));
         const double Fq = cos(m*t + b) * radical * bj;
         total += Gauss76Wt[i] * Fq;
@@ -50 +50 @@
 
 static double
-_fq(double q, double h, double radius_cap, double radius, double half_length,
-    double sin_alpha, double cos_alpha)
+_fq(double qab, double qc, double h, double radius_cap, double radius, double half_length)
 {
-    const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha);
-    const double bj = sas_2J1x_x(q*radius*sin_alpha);
-    const double si = sas_sinx_x(q*half_length*cos_alpha);
+    const double cap_Fq = _cap_kernel(qab, qc, h, radius_cap, half_length);
+    const double bj = sas_2J1x_x(radius*qab);
+    const double si = sas_sinx_x(half_length*qc);
     const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si;
     const double Aq = cap_Fq + cyl_Fq;
@@ -101 +100 @@
     double total = 0.0;
     for (int i=0; i<76 ;i++) {
-        const double alpha= Gauss76Z[i]*zm + zb;
-        double sin_alpha, cos_alpha; // slots to hold sincos function output
-        SINCOS(alpha, sin_alpha, cos_alpha);
-
-        const double Aq = _fq(q, h, radius_cap, radius, half_length, sin_alpha, cos_alpha);
-        // sin_alpha for spherical coord integration
-        total += Gauss76Wt[i] * Aq * Aq * sin_alpha;
+        const double theta = Gauss76Z[i]*zm + zb;
+        double sin_theta, cos_theta; // slots to hold sincos function output
+        SINCOS(theta, sin_theta, cos_theta);
+        const double qab = q*sin_theta;
+        const double qc = q*cos_theta;
+        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length);
+        // scale by sin_theta for spherical coord integration
+        total += Gauss76Wt[i] * Aq * Aq * sin_theta;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -125 +125 @@
     double q, sin_alpha, cos_alpha;
     ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
     const double h = sqrt(radius_cap*radius_cap - radius*radius);
-    const double Aq = _fq(q, h, radius_cap, radius, 0.5*length, sin_alpha, cos_alpha);
+    const double Aq = _fq(qab, qc, h, radius_cap, radius, 0.5*length);
 
     // Multiply by contrast^2 and convert to cm-1