Changeset 5d4777d in sasmodels for sasmodels/models/cylinder.c

Timestamp:

Sep 2, 2014 1:24:38 AM (10 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:

f4cf580

Parents:

ff7119b

Message:

reorganize, check and update models

File:

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

sasmodels/models/cylinder.c

@@ -1 +1 @@
 real form_volume(real radius, real length);
 real Iq(real q, real sld, real solvent_sld, real radius, real length);
-real Iqxy(real qx, real qy, real sld, real solvent_sld, real radius, real length, real theta, real phi);
+real Iqxy(real qx, real qy, real sld, real solvent_sld,
+    real radius, real length, real theta, real phi);
+
+// twovd = 2 * volume * delta_rho
+// besarg = q * R * sin(alpha)
+// siarg = q * L/2 * cos(alpha)
+real _cyl(real twovd, real besarg, real siarg);
+real _cyl(real twovd, real besarg, real siarg)
+{
+    const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg);
+    const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg);
+    return twovd*si*bj;
+}
 
 real form_volume(real radius, real length)
@@ -14 +26 @@
     real length)
 {
-    const real halflength = REAL(0.5)*length;
-    real summ = REAL(0.0);
+    const real qr = q*radius;
+    const real qh = q*REAL(0.5)*length;
+    const real twovd = REAL(2.0)*(sld-solvent_sld)*form_volume(radius, length);
+    real total = REAL(0.0);
     // real 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 real zi = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
-        const real zi = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2);
-        summ += Gauss76Wt[i] * CylKernel(q, radius, halflength, zi);
+        //const real alpha = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0;
+        const real alpha = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2);
+        real sn, cn;
+        SINCOS(alpha, sn, cn);
+        const real fq = _cyl(twovd, qr*sn, qh*cn);
+        total += Gauss76Wt[i] * fq * fq * sn;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
-    //const real form = (upper-lower)/2.0*summ;
-    const real form = summ * M_PI_4;
-
-    // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1
-    // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
-    // The additional volume factor is for polydisperse volume normalization.
-    const real s = (sld - solvent_sld) * form_volume(radius, length);
-    return REAL(1.0e-4) * form * s * s;
+    //const real form = (upper-lower)/2.0*total;
+    return REAL(1.0e-4) * total * M_PI_4;
 }
 
@@ -43 +54 @@
     real 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.
     real sn, cn; // slots to hold sincos function output
 
@@ -54 +68 @@
     const real alpha = acos(cos_val);
 
-    // The following is CylKernel() / sin(alpha), but we are doing it in place
-    // to avoid sin(alpha)/sin(alpha) for alpha = 0.  It is also a teensy bit
-    // faster since we don't mulitply and divide sin(alpha).
+    const real twovd = REAL(2.0)*(sld-solvent_sld)*form_volume(radius, length);
     SINCOS(alpha, sn, cn);
-    const real besarg = q*radius*sn;
-    const real siarg = REAL(0.5)*q*length*cn;
-    // lim_{x->0} J1(x)/x = 1/2,   lim_{x->0} sin(x)/x = 1
-    const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg);
-    const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg);
-    const real form = REAL(4.0)*bj*bj*si*si;
-
-    // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1
-    // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
-    // The additional volume factor is for polydisperse volume normalization.
-    const real s = (sld - solvent_sld) * form_volume(radius, length);
-    return REAL(1.0e-4) * form * s * s; // * correction;
+    const real fq = _cyl(twovd, q*radius*sn, q*REAL(0.5)*length*cn);
+    return REAL(1.0e-4) * fq * fq; // * correction;
 }