Changeset 6831fa0 in sasmodels

Timestamp:

Oct 14, 2016 3:18:34 PM (8 years ago)

Author:

Paul Kienzle <pkienzle@…>

Branches:

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

Children:

b716cc6

Parents:

87bc707

Message:

stacked_disks: enable 2D modeling for stacked disks

Location:

sasmodels

Files:

• compare.py (modified) (3 diffs)
• conversion_table.py (modified) (1 diff)
• models/stacked_disks.c (modified) (8 diffs)

sasmodels/compare.py

@@ -38 +38 @@
 from . import core
 from . import kerneldll
-from . import weights
+from . import exception
 from .data import plot_theory, empty_data1D, empty_data2D
 from .direct_model import DirectModel
@@ -45 +45 @@
 try:
     from typing import Optional, Dict, Any, Callable, Tuple
-except:
+except Exception:
     pass
 else:
@@ -466 +466 @@
             if k.endswith('.type'):
                 par = k[:-5]
+                if v == 'gaussian': continue
                 cls = dispersers[v if v != 'rectangle' else 'rectangula']
                 handle = cls()
                 model[0].disperser_handles[par] = handle
-                model[0].set_dispersion(par, handle)
+                try:
+                    model[0].set_dispersion(par, handle)
+                except Exception:
+                    exception.annotate_exception("while setting %s to %r"
+                                                 %(par, v))
+                    raise
+
 
         #print("sasview pars",oldpars)

sasmodels/conversion_table.py

@@ -729 +729 @@
             "theta": "axis_theta",
             "sld_solvent": "solvent_sld",
-            "n_stacking": "n_stacking"
+            "n_stacking": "n_stacking",
+            "thick_layer": "layer_thick",
+            "thick_core": "core_thick",
         }
     ],

sasmodels/models/stacked_disks.c

@@ -14 +14 @@
           double solvent_sld);
 
+double Iqxy(double qx, double qy,
+          double thick_core,
+          double thick_layer,
+          double radius,
+          double n_stacking,
+          double sigma_dnn,
+          double core_sld,
+          double layer_sld,
+          double solvent_sld,
+          double theta,
+          double phi);
+
 static
 double _kernel(double qq,
@@ -22 +34 @@
                double halfheight,
                double thick_layer,
-               double zi,
+               double sin_alpha,
+               double cos_alpha,
                double sigma_dnn,
                double d,
@@ -34 +47 @@
         // zi is the dummy variable for the integration (x in Feigin's notation)
 
-        const double besarg1 = qq*radius*sin(zi);
-        const double besarg2 = qq*radius*sin(zi);
+        const double besarg1 = qq*radius*sin_alpha;
+        //const double besarg2 = qq*radius*sin_alpha;
 
-        const double sinarg1 = qq*halfheight*cos(zi);
-        const double sinarg2 = qq*(halfheight+thick_layer)*cos(zi);
+        const double sinarg1 = qq*halfheight*cos_alpha;
+        const double sinarg2 = qq*(halfheight+thick_layer)*cos_alpha;
 
         const double be1 = sas_J1c(besarg1);
-        const double be2 = sas_J1c(besarg2);
-        const double si1 = sin(sinarg1)/sinarg1;
-        const double si2 = sin(sinarg2)/sinarg2;
+        //const double be2 = sas_J1c(besarg2);
+        const double be2 = be1;
+        const double si1 = sinc(sinarg1);
+        const double si2 = sinc(sinarg2);
 
         const double dr1 = (core_sld-solvent_sld);
         const double dr2 = (layer_sld-solvent_sld);
         const double area = M_PI*radius*radius;
-        const double totald=2.0*(thick_layer+halfheight);
+        const double totald = 2.0*(thick_layer+halfheight);
 
         const double t1 = area*(2.0*halfheight)*dr1*(si1)*(be1);
@@ -54 +68 @@
 
 
-        double retval =((t1+t2)*(t1+t2))*sin(zi);
+        double retval =((t1+t2)*(t1+t2));
 
         // loop for the structure facture S(q)
         double sqq=0.0;
         for(int kk=1;kk<n_stacking;kk+=1) {
-                double dexpt=qq*cos(zi)*qq*cos(zi)*d*d*sigma_dnn*sigma_dnn*kk/2.0;
-                sqq=sqq+(n_stacking-kk)*cos(qq*cos(zi)*d*kk)*exp(-1.*dexpt);
+                double qd_cos_alpha = qq*d*cos_alpha;
+                double dexpt=square(qd_cos_alpha*sigma_dnn)*kk/2.0;
+                sqq += (n_stacking-kk)*cos(qd_cos_alpha*kk)*exp(-1.*dexpt);
         }
-
         // end of loop for S(q)
         sqq=1.0+2.0*sqq/n_stacking;
 
-        retval *= sqq;
-
-        return(retval);
+        return retval * sqq;
 }
 
@@ -88 +100 @@
         double summ = 0.0;      //initialize integral
 
-        double d=2.0*thick_layer+thick_core;
-        double halfheight = thick_core/2.0;
+        double d = 2.0*thick_layer+thick_core;
+        double halfheight = 0.5*thick_core;
 
-        for(int i=0;i<N_POINTS_76;i++) {
-                double zi = (Gauss76Z[i] + 1.0)*M_PI/4.0;
-                double yyy = Gauss76Wt[i] *
-                    _kernel(q,
+        for(int i=0; i<N_POINTS_76; i++) {
+                double zi = (Gauss76Z[i] + 1.0)*M_PI_4;
+                double sin_alpha, cos_alpha; // slots to hold sincos function output
+                SINCOS(zi, sin_alpha, cos_alpha);
+                double yyy = _kernel(q,
                                    radius,
                                    core_sld,
@@ -101 +114 @@
                                    halfheight,
                                    thick_layer,
-                                   zi,
+                                   sin_alpha,
+                                   cos_alpha,
                                    sigma_dnn,
                                    d,
                                    n_stacking);
-                summ += yyy;
+                summ += Gauss76Wt[i] * yyy * sin_alpha;
         }
 
-        double answer = M_PI/4.0*summ;
+        double answer = M_PI_4*summ;
 
         //Convert to [cm-1]
@@ -116 +130 @@
 }
 
-static double stacked_disks_kernel_2d(double q, double q_x, double q_y,
-                            double thick_core,
-                            double thick_layer,
-                            double radius,
-                            double n_stacking,
-                            double sigma_dnn,
-                            double core_sld,
-                            double layer_sld,
-                            double solvent_sld,
-                            double theta,
-                            double phi)
-{
-
-    double ct, st, cp, sp;
-
-    //convert angle degree to radian
-    theta = theta * M_PI/180.0;
-    phi = phi * M_PI/180.0;
-
-    SINCOS(theta, st, ct);
-    SINCOS(phi, sp, cp);
-
-    // silence compiler warnings about unused variable
-    (void) sp;
-
-    // parallelepiped orientation
-    const double cyl_x = st * cp;
-    const double cyl_y = st * sp;
-
-    // Compute the angle btw vector q and the
-    // axis of the parallelepiped
-    const double cos_val = cyl_x*q_x + cyl_y*q_y;
-
-    // Note: cos(alpha) = 0 and 1 will get an
-    // undefined value from Stackdisc_kern
-    double alpha = acos( cos_val );
-
-    // Call the IGOR library function to get the kernel
-    double d = 2 * thick_layer + thick_core;
-    double halfheight = thick_core/2.0;
-    double answer = _kernel(q,
-                     radius,
-                     core_sld,
-                     layer_sld,
-                     solvent_sld,
-                     halfheight,
-                     thick_layer,
-                     alpha,
-                     sigma_dnn,
-                     d,
-                     n_stacking);
-
-    answer /= sin(alpha);
-    //convert to [cm-1]
-    answer *= 1.0e-4;
-
-    return answer;
-}
-
 double form_volume(double thick_core,
                    double thick_layer,
                    double radius,
                    double n_stacking){
-    double d = 2 * thick_layer + thick_core;
-    return acos(-1.0) * radius * radius * d * n_stacking;
+    double d = 2.0 * thick_layer + thick_core;
+    return M_PI * radius * radius * d * n_stacking;
 }
 
@@ -203 +158 @@
                     solvent_sld);
 }
+
+
+double
+Iqxy(double qx, double qy,
+     double thick_core,
+     double thick_layer,
+     double radius,
+     double n_stacking,
+     double sigma_dnn,
+     double core_sld,
+     double layer_sld,
+     double solvent_sld,
+     double theta,
+     double phi)
+{
+    double q, sin_alpha, cos_alpha;
+    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
+
+    double d = 2.0 * thick_layer + thick_core;
+    double halfheight = 0.5*thick_core;
+    double answer = _kernel(q,
+                     radius,
+                     core_sld,
+                     layer_sld,
+                     solvent_sld,
+                     halfheight,
+                     thick_layer,
+                     sin_alpha,
+                     cos_alpha,
+                     sigma_dnn,
+                     d,
+                     n_stacking);
+
+    //convert to [cm-1]
+    answer *= 1.0e-4;
+
+    return answer;
+}
+