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

sasmodels/models/core_shell_ellipsoid.c

@@ -1 @@
-double form_volume(double radius_equat_core,
-                   double polar_core,
-                   double equat_shell,
-                   double polar_shell);
-double Iq(double q,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld);
 
+// Converted from Igor function gfn4, using the same pattern as ellipsoid
+// for evaluating the parts of the integral.
+//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN
+//                       BY (53) & (58-59) IN CHEN AND
+//                       KOTLARCHYK REFERENCE
+//
+//       <OBLATE ELLIPSOID>
+static double
+_cs_ellipsoid_kernel(double qab, double qc,
+    double equat_core, double polar_core,
+    double equat_shell, double polar_shell,
+    double sld_core_shell, double sld_shell_solvent)
+{
+    const double qr_core = sqrt(square(equat_core*qab) + square(polar_core*qc));
+    const double si_core = sas_3j1x_x(qr_core);
+    const double volume_core = M_4PI_3*equat_core*equat_core*polar_core;
+    const double fq_core = si_core*volume_core*sld_core_shell;
 
-double Iqxy(double qx, double qy,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double theta,
-          double phi);
+    const double qr_shell = sqrt(square(equat_shell*qab) + square(polar_shell*qc));
+    const double si_shell = sas_3j1x_x(qr_shell);
+    const double volume_shell = M_4PI_3*equat_shell*equat_shell*polar_shell;
+    const double fq_shell = si_shell*volume_shell*sld_shell_solvent;
 
+    return fq_core + fq_shell;
+}
 
-double form_volume(double radius_equat_core,
-                   double x_core,
-                   double thick_shell,
-                   double x_polar_shell)
+static double
+form_volume(double radius_equat_core,
+    double x_core,
+    double thick_shell,
+    double x_polar_shell)
 {
     const double equat_shell = radius_equat_core + thick_shell;
@@ -37 +39 @@
 
 static double
-core_shell_ellipsoid_xt_kernel(double q,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld)
+Iq(double q,
+    double radius_equat_core,
+    double x_core,
+    double thick_shell,
+    double x_polar_shell,
+    double core_sld,
+    double shell_sld,
+    double solvent_sld)
 {
-    const double lolim = 0.0;
-    const double uplim = 1.0;
-
-
-    const double delpc = core_sld - shell_sld; //core - shell
-    const double delps = shell_sld - solvent_sld; //shell - solvent
-
+    const double sld_core_shell = core_sld - shell_sld;
+    const double sld_shell_solvent = shell_sld - solvent_sld;
 
     const double polar_core = radius_equat_core*x_core;
@@ -58 +55 @@
     const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell;
 
-    double summ = 0.0;   //initialize intergral
+    // translate from [-1, 1] => [0, 1]
+    const double m = 0.5;
+    const double b = 0.5;
+    double total = 0.0;     //initialize intergral
     for(int i=0;i<76;i++) {
-        double zi = 0.5*( Gauss76Z[i]*(uplim-lolim) + uplim + lolim );
-        double yyy = gfn4(zi, radius_equat_core, polar_core, equat_shell,
-                          polar_shell, delpc, delps, q);
-        summ += Gauss76Wt[i] * yyy;
+        const double cos_theta = Gauss76Z[i]*m + b;
+        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
+        double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta,
+            radius_equat_core, polar_core,
+            equat_shell, polar_shell,
+            sld_core_shell, sld_shell_solvent);
+        total += Gauss76Wt[i] * fq * fq;
     }
-    summ *= 0.5*(uplim-lolim);
+    total *= m;
 
     // convert to [cm-1]
-    return 1.0e-4 * summ;
+    return 1.0e-4 * total;
 }
 
 static double
-core_shell_ellipsoid_xt_kernel_2d(double qx, double qy,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double theta,
-          double phi)
+Iqxy(double qx, double qy,
+    double radius_equat_core,
+    double x_core,
+    double thick_shell,
+    double x_polar_shell,
+    double core_sld,
+    double shell_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);
+    const double qab = q*sin_alpha;
+    const double qc = q*cos_alpha;
 
-    const double sldcs = core_sld - shell_sld;
-    const double sldss = shell_sld- solvent_sld;
+    const double sld_core_shell = core_sld - shell_sld;
+    const double sld_shell_solvent = shell_sld - solvent_sld;
 
     const double polar_core = radius_equat_core*x_core;
@@ -93 +98 @@
     const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell;
 
-    // Call the IGOR library function to get the kernel:
-    // MUST use gfn4 not gf2 because of the def of params.
-    double answer = gfn4(cos_alpha,
-                  radius_equat_core,
-                  polar_core,
-                  equat_shell,
-                  polar_shell,
-                  sldcs,
-                  sldss,
-                  q);
+    double fq = _cs_ellipsoid_kernel(qab, qc,
+                  radius_equat_core, polar_core,
+                  equat_shell, polar_shell,
+                  sld_core_shell, sld_shell_solvent);
 
     //convert to [cm-1]
-    answer *= 1.0e-4;
-
-    return answer;
+    return 1.0e-4 * fq * fq;
 }
-
-double Iq(double q,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld)
-{
-    double intensity = core_shell_ellipsoid_xt_kernel(q,
-           radius_equat_core,
-           x_core,
-           thick_shell,
-           x_polar_shell,
-           core_sld,
-           shell_sld,
-           solvent_sld);
-
-    return intensity;
-}
-
-
-double Iqxy(double qx, double qy,
-          double radius_equat_core,
-          double x_core,
-          double thick_shell,
-          double x_polar_shell,
-          double core_sld,
-          double shell_sld,
-          double solvent_sld,
-          double theta,
-          double phi)
-{
-    double intensity = core_shell_ellipsoid_xt_kernel_2d(qx, qy,
-                       radius_equat_core,
-                       x_core,
-                       thick_shell,
-                       x_polar_shell,
-                       core_sld,
-                       shell_sld,
-                       solvent_sld,
-                       theta,
-                       phi);
-
-    return intensity;
-}