Changeset 904cd9c in sasmodels for sasmodels/models/core_shell_parallelepiped.c

Timestamp:

Oct 30, 2017 3:21:00 PM (6 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:

bcb5594

Parents:

5b5ea20 (diff), ef969d9 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.

Message:

Merge branch 'master' into ticket-776-orientation

File:

• sasmodels/models/core_shell_parallelepiped.c (modified) (9 diffs)

sasmodels/models/core_shell_parallelepiped.c

@@ -26 +26 @@
     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled
     // Did not understand the code completely, it should be rechecked (Miguel Gonzalez)
+    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker)
 
     const double mu = 0.5 * q * length_b;
@@ -40 +41 @@
     const double c_scaled = length_c / length_b;
 
-    // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG)
-    // the a_scaled in the definition of tb was present in CSPPKernel in libCylinder.c,
-    // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions
-    // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!)
-    double ta = (a_scaled + 2.0*thick_rim_a)/length_b;
-    double tb = (a_scaled + 2.0*thick_rim_b)/length_b;
+    double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b;
+    double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b;
+    double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present
 
     double Vin = length_a * length_b * length_c;
@@ -54 +52 @@
     double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc)
     double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc)
+    double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present
 
     // Scale factors (note that drC is not used later)
@@ -59 +58 @@
     const double drhoA = (arim_sld-solvent_sld);
     const double drhoB = (brim_sld-solvent_sld);
-    //const double drC_Vot = (crim_sld-solvent_sld)*Vot;
+    const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot;
+
 
     // Precompute scale factors for combining cross terms from the shape
     const double scale23 = drhoA*V1;
     const double scale14 = drhoB*V2;
-    const double scale12 = drho0*Vin - scale23 - scale14;
+    const double scale24 = drhoC*V3;
+    const double scale11 = drho0*Vin;
+    const double scale12 = drho0*Vin - scale23 - scale14 - scale24;
 
     // outer integral (with gauss points), integration limits = 0, 1
@@ -75 +77 @@
         // inner integral (with gauss points), integration limits = 0, 1
         double inner_total = 0.0;
+        double inner_total_crim = 0.0;
         for(int j=0; j<76; j++) {
             const double uu = 0.5 * ( Gauss76Z[j] + 1.0 );
@@ -80 +83 @@
             SINCOS(M_PI_2*uu, sin_uu, cos_uu);
             const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);
-            const double si2 = sas_sinx_x(mu_proj * cos_uu);
+            const double si2 = sas_sinx_x(mu_proj * cos_uu );
             const double si3 = sas_sinx_x(mu_proj * sin_uu * ta);
             const double si4 = sas_sinx_x(mu_proj * cos_uu * tb);
@@ -86 +89 @@
             // Expression in libCylinder.c (neither drC nor Vot are used)
             const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4;
-
-            // To note also that in csparallelepiped.cpp, there is a function called
-            // cspkernel, which seems to make more sense and has the following comment:
-            //   This expression is different from NIST/IGOR package (I strongly believe the IGOR is wrong!!!). 10/15/2010.
-            //   tmp =( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3)*
-            //   ( dr0*tmp1*tmp2*tmp3*Vin + drA*(tmpt1-tmp1)*tmp2*tmp3*V1+ drB*tmp1*(tmpt2-tmp2)*tmp3*V2 + drC*tmp1*tmp2*(tmpt3-tmp3)*V3);   //  correct FF : square of sum of phase factors
-            // This is the function called by csparallelepiped_analytical_2D_scaled,
-            // while CSParallelepipedModel calls CSParallelepiped in libCylinder.c
+            const double form_crim = scale11*si1*si2;
 
             //  correct FF : sum of square of phase factors
             inner_total += Gauss76Wt[j] * form * form;
+            inner_total_crim += Gauss76Wt[j] * form_crim * form_crim;
         }
         inner_total *= 0.5;
-
+        inner_total_crim *= 0.5;
         // now sum up the outer integral
         const double si = sas_sinx_x(mu * c_scaled * sigma);
-        outer_total += Gauss76Wt[i] * inner_total * si * si;
+        const double si_crim = sas_sinx_x(mu * tc * sigma);
+        outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim);
     }
     outer_total *= 0.5;
@@ -141 +139 @@
 
     // The definitions of ta, tb, tc are not the same as in the 1D case because there is no
-    // the scaling by B. The use of length_a for the 3 of them seems clearly a mistake to me,
-    // but for the moment I let it like this until understanding better the code.
+    // the scaling by B.
     double ta = length_a + 2.0*thick_rim_a;
-    double tb = length_a + 2.0*thick_rim_b;
-    double tc = length_a + 2.0*thick_rim_c;
+    double tb = length_b + 2.0*thick_rim_b;
+    double tc = length_c + 2.0*thick_rim_c;
     //handle arg=0 separately, as sin(t)/t -> 1 as t->0
     double siA = sas_sinx_x(0.5*length_a*qa);
@@ -160 +157 @@
                + drA*(siAt-siA)*siB*siC*V1
                + drB*siA*(siBt-siB)*siC*V2
-               + drC*siA*siB*(siCt*siCt-siC)*V3);
+               + drC*siA*siB*(siCt-siC)*V3);
 
     return 1.0e-4 * f * f;