Changeset f4cf580 in sasmodels for sasmodels

Timestamp:

Sep 2, 2014 1:15:40 PM (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:

87985ca

Parents:

5d4777d

Message:

resolve remaining differences between sasview and sasmodels

Location:

sasmodels

Files:

• convert.py (modified) (2 diffs)
• gen.py (modified) (4 diffs)
• gpu.py (modified) (2 diffs)
• models/capped_cylinder.c (modified) (4 diffs)
• models/capped_cylinder.py (modified) (2 diffs)
• models/sphere.py (modified) (1 diff)

sasmodels/convert.py

@@ -37 +37 @@
         name='sphere',
         sldSph='sld', sldSolv='solvent_sld',
-        #radius='radius',  # DO NOT LIST IDENTICAL PARAMETERS!
+        radius='radius',  # listing identical parameters is optional
         ),
     }
@@ -59 +59 @@
     newpars = pars.copy()
     for old,new in mapping.items():
+        if old == new: continue
         # Bumps style parameter names
         for variant in ("", "_pd", "_pd_n", "_pd_nsigma", "_pd_type"):

sasmodels/gen.py

@@ -352 +352 @@
   const real I = %(fn)s(%(qcall)s, %(pcall)s);
   if (I>=REAL(0.0)) { // scattering cannot be negative
-    ret += weight*I;
+    ret += weight*I%(sasview_spherical)s;
     norm += weight;
     %(volume_norm)s
@@ -364 +364 @@
 SPHERICAL_CORRECTION="""\
 // Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta
-real spherical_correction = (Ntheta>1 ? fabs(cos(M_PI_180*phi)) : REAL(1.0));\
+real spherical_correction = (Ntheta>1 ? fabs(sin(M_PI_180*theta)) : REAL(1.0));\
+"""
+# Use this to reproduce sasview behaviour
+SASVIEW_SPHERICAL_CORRECTION="""\
+// Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta
+real spherical_correction = (Ntheta>1 ? fabs(cos(M_PI_180*theta))*M_PI_2 : REAL(1.0));\
 """
 
@@ -472 +477 @@
     fq_pars = [p[0] for p in info['parameters'][len(COMMON_PARAMETERS):]
                if p[0] in set(fixed_pars+pd_pars)]
-    if False and "phi" in pd_pars:
+    if False and "theta" in pd_pars:
         spherical_correction = [indent(SPHERICAL_CORRECTION, depth)]
         weights = [p+"_w" for p in pd_pars]+['spherical_correction']
+        sasview_spherical = ""
+    elif "theta" in pd_pars:
+        spherical_correction = [indent(SASVIEW_SPHERICAL_CORRECTION,depth)]
+        weights = [p+"_w" for p in pd_pars]
+        sasview_spherical = "*spherical_correction"
     else:
         spherical_correction = []
         weights = [p+"_w" for p in pd_pars]
+        sasview_spherical = ""
     subst = {
         'weight_product': "*".join(weights),
@@ -484 +495 @@
         'qcall': q_pars['qcall'],
         'pcall': ", ".join(fq_pars), # skip scale and background
+        'sasview_spherical': sasview_spherical,
         }
     loop_body = [indent(LOOP_BODY%subst, depth)]

sasmodels/gpu.py

@@ -22 +22 @@
 devices, where it can be combined with other structure factors and form
 factors and have instrumental resolution effects applied.
-
-
 """
 import numpy as np
@@ -59 +57 @@
     source, info = gen.make(kernel_module)
     ## for debugging, save source to a .cl file, edit it, and reload as model
-    open(info['name']+'.cl','w').write(source)
+    #open(info['name']+'.cl','w').write(source)
     #source = open(info['name']+'.cl','r').read()
     return GpuModel(source, info, dtype)

sasmodels/models/capped_cylinder.c

@@ -22 +22 @@
     const real upper = REAL(1.0);
     const real lower = h/cap_radius; // integral lower bound
-    const real m = q*cos_alpha*cap_radius; // cos argument slope
-    const real b = q*cos_alpha*(REAL(0.5)*length-h); // cos argument intercept
-    const real qrst = q*sin_alpha*cap_radius; // Q*R*sin(theta)
+    // cos term in integral is:
+    //    cos (q (R t - h + L/2) cos(alpha))
+    // so turn it into:
+    //    cos (m t + b)
+    // where:
+    //    m = q R cos(alpha)
+    //    b = q(L/2-h) cos(alpha)
+    const real m = q*cap_radius*cos_alpha; // cos argument slope
+    const real b = q*(REAL(0.5)*length-h)*cos_alpha; // cos argument intercept
+    const real qrst = q*cap_radius*sin_alpha; // Q*R*sin(theta)
     real total = REAL(0.0);
     for (int i=0; i<76 ;i++) {
@@ -31 +38 @@
         const real t = REAL(0.5)*(Gauss76Z[i]*(upper-lower)+upper+lower);
         const real radical = REAL(1.0) - t*t;
-        const real caparg = qrst*sqrt(radical); // cap bessel function arg
-        const real be = (caparg == REAL(0.0) ? REAL(0.5) : J1(caparg)/caparg);
+        const real arg = qrst*sqrt(radical); // cap bessel function arg
+        const real be = (arg == REAL(0.0) ? REAL(0.5) : J1(arg)/arg);
         const real Fq = cos(m*t + b) * radical * be;
         total += Gauss76Wt[i] * Fq;
@@ -84 +91 @@
         // faster since we don't multiply and divide sin(alpha).
         const real besarg = q*radius*sn;
-        const real siarg = REAL(0.5)*q*length*cn;
+        const real siarg = q*REAL(0.5)*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);
@@ -131 +138 @@
     const real cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn);
 
-    // 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 multiply and divide sin(alpha).
     const real besarg = q*radius*sn;
-    const real siarg = REAL(0.5)*q*length*cn;
+    const real siarg = q*REAL(0.5)*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);

sasmodels/models/capped_cylinder.py

@@ -1 +1 @@
 r"""
 Calculates the scattering from a cylinder with spherical section end-caps.
-This model simply becomes the `convex_lens`_ when the length of the cylinder
+This model simply becomes the a convex lens when the length of the cylinder
 $L=0$, that is, a sphereocylinder with end caps that have a radius larger
 than that of the cylinder and the center of the end cap radius lies within
@@ -128 +128 @@
     [ "radius", "Ang",  20, [0, inf], "volume",
       "Cylinder radius" ],
+    # TODO: use an expression for cap radius with fixed bounds.
+    # The current form requires cap radius R bigger than cylinder radius r.
+    # Could instead use R/r in [1,inf], r/R in [0,1], or the angle between
+    # cylinder and cap in [0,90].  The problem is similar for the barbell
+    # model.  Propose r/R in [0,1] in both cases, with the model specifying
+    # cylinder radius in the capped cylinder model and sphere radius in the
+    # barbell model.  This leads to the natural value of zero for no cap
+    # in the capped cylinder, and zero for no bar in the barbell model.  In
+    # both models, one would be a pill.
     [ "cap_radius", "Ang", 20, [0, inf], "volume",
       "Cap radius" ],

sasmodels/models/sphere.py

@@ -30 +30 @@
 
 Iq = """
+    const real qr = q*radius;
     real sn, cn;
-    const real qr = q*radius;
     SINCOS(qr, sn, cn);
-    const real bes = (qr==REAL(0.0) ? REAL(1.0) : (REAL(3.0)*(sn-qr*cn))/(qr*qr*qr));
-    const real f = bes * (sld - solvent_sld) * form_volume(radius);
-    return REAL(1.0e-4)*f*f;
+    const real bes = qr==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(sn-qr*cn)/(qr*qr*qr);
+    const real fq = bes * (sld - solvent_sld) * form_volume(radius);
+    return REAL(1.0e-4)*fq*fq;
     """