Changeset 71b751d in sasmodels for sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

Timestamp:

Aug 14, 2018 12:09:34 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:

86aa992

Parents:

2f8cbb9

Message:

update remaining form factors to use Fq interface

File:

• sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c (modified) (5 diffs)

sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

@@ -11 +11 @@
 }
 
-static double
-Iq(double q,
+static void
+Fq(double q,
+        double *F1,
+        double *F2,
         double r_minor,
         double x_core,
@@ -24 +26 @@
         double sigma)
 {
-    double si1,si2,be1,be2;
      // core_shell_bicelle_elliptical_belt, RKH 5th Oct 2017, core_shell_bicelle_elliptical
      // tested briefly against limiting cases of cylinder, hollow cylinder & elliptical cylinder models
@@ -38 +39 @@
     const double r2A = 0.5*(square(r_major) + square(r_minor));
     const double r2B = 0.5*(square(r_major) - square(r_minor));
+    const double vol1 = M_PI*r_minor*r_major*(2.0*halfheight);
+    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*halfheight;
+    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
     // dr1,2,3 are now for Vcore, Vcore+rim, Vcore+face,
-    const double dr1 = (-rhor - rhoh + rhoc + rhosolv) *M_PI*r_minor*r_major*
-          2.0*halfheight;
-    const double dr2 = (rhor-rhosolv) *M_PI*(r_minor+thick_rim)*(
-         r_major+thick_rim)* 2.0*halfheight;
-    const double dr3 = (rhoh-rhosolv) *M_PI*r_minor*r_major*
-         2.0*(halfheight+thick_face);
+    const double dr1 = vol1*(-rhor - rhoh + rhoc + rhosolv);
+    const double dr2 = vol2*(rhor-rhosolv);
+    const double dr3 = vol3*(rhoh-rhosolv);
+
     //initialize integral
-    double outer_sum = 0.0;
+    double outer_total_F1 = 0.0;
+    double outer_total_F2 = 0.0;
     for(int i=0;i<GAUSS_N;i++) {
         //setup inner integral over the ellipsoidal cross-section
         // since we generate these lots of times, why not store them somewhere?
-        //const double cos_alpha = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
-        const double cos_alpha = ( GAUSS_Z[i] + 1.0 )/2.0;
-        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
-        double inner_sum=0;
-        double sinarg1 = q*halfheight*cos_alpha;
-        double sinarg2 = q*(halfheight+thick_face)*cos_alpha;
-        si1 = sas_sinx_x(sinarg1);
-        si2 = sas_sinx_x(sinarg2);
+        //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0;
+        const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0;
+        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
+        const double qab = q*sin_theta;
+        const double qc = q*cos_theta;
+        const double si1 = sas_sinx_x(halfheight*qc);
+        const double si2 = sas_sinx_x((halfheight+thick_face)*qc);
+        double inner_total_F1 = 0;
+        double inner_total_F2 = 0;
         for(int j=0;j<GAUSS_N;j++) {
             //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe)
@@ -63 +67 @@
             const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2;
             const double rr = sqrt(r2A - r2B*cos(beta));
-            double besarg1 = q*rr*sin_alpha;
-            double besarg2 = q*(rr+thick_rim)*sin_alpha;
-            be1 = sas_2J1x_x(besarg1);
-            be2 = sas_2J1x_x(besarg2);
-            inner_sum += GAUSS_W[j] *square(dr1*si1*be1 +
-                                              dr2*si1*be2 +
-                                              dr3*si2*be1);
+            const double be1 = sas_2J1x_x(rr*qab);
+            const double be2 = sas_2J1x_x((rr+thick_rim)*qab);
+            const double f = dr1*si1*be1 + dr2*si1*be2 + dr3*si2*be1;
+
+            inner_total_F1 += GAUSS_W[j] * f;
+            inner_total_F2 += GAUSS_W[j] * f * f;
         }
         //now calculate outer integral
-        outer_sum += GAUSS_W[i] * inner_sum;
+        outer_total_F1 += GAUSS_W[i] * inner_total_F1;
+        outer_total_F2 += GAUSS_W[i] * inner_total_F2;
     }
+    // now complete change of integration variables (1-0)/(1-(-1))= 0.5
+    outer_total_F1 *= 0.25;
+    outer_total_F2 *= 0.25;
 
-    return outer_sum*2.5e-05*exp(-0.5*square(q*sigma));
+    //convert from [1e-12 A-1] to [cm-1]
+    *F1 = 1e-2*outer_total_F1*exp(-0.25*square(q*sigma));
+    *F2 = 1e-4*outer_total_F2*exp(-0.5*square(q*sigma));
 }
 
@@ -111 +120 @@
     const double si1 = sas_sinx_x( halfheight*qc );
     const double si2 = sas_sinx_x( (halfheight + thick_face)*qc );
-    const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si1*be2 +  vol3*dr3*si2*be1);
-    return 1.0e-4 * Aq*exp(-0.5*(square(qa) + square(qb) + square(qc) )*square(sigma));
+    const double fq = vol1*dr1*si1*be1 + vol2*dr2*si1*be2 +  vol3*dr3*si2*be1;
+    const double atten = exp(-0.5*(qa*qa + qb*qb + qc*qc)*(sigma*sigma));
+    return 1.0e-4 * fq*fq * atten;
 }