Changes in sasmodels/models/core_shell_bicelle_elliptical.c [592343f:dedcf34] in sasmodels

File:

• sasmodels/models/core_shell_bicelle_elliptical.c (modified) (6 diffs)

sasmodels/models/core_shell_bicelle_elliptical.c

@@ -1 @@
-double form_volume(double radius, double x_core, double thick_rim, double thick_face, double length);
-double Iq(double q,
-          double radius,
-          double x_core,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld);
-
-
-double Iqxy(double qx, double qy,
-          double radius,
-          double x_core,
-          double thick_rim,
-          double thick_face,
-          double length,
-          double core_sld,
-          double face_sld,
-          double rim_sld,
-          double solvent_sld,
-          double theta,
-          double phi,
-          double psi);
-
 // NOTE that "length" here is the full height of the core!
-double form_volume(double radius, double x_core, double thick_rim, double thick_face, double length)
+static double
+form_volume(double r_minor,
+        double x_core,
+        double thick_rim,
+        double thick_face,
+        double length)
 {
-    return M_PI*(radius+thick_rim)*(radius*x_core+thick_rim)*(length+2.0*thick_face);
+    return M_PI*(r_minor+thick_rim)*(r_minor*x_core+thick_rim)*(length+2.0*thick_face);
 }
 
-double 
-                Iq(double qq,
-                   double rad,
-                   double x_core,
-                   double radthick,
-                   double facthick,
-                   double length,
-                   double rhoc,
-                   double rhoh,
-                   double rhor,
-                   double rhosolv)
+static double
+Iq(double q,
+        double r_minor,
+        double x_core,
+        double thick_rim,
+        double thick_face,
+        double length,
+        double rhoc,
+        double rhoh,
+        double rhor,
+        double rhosolv)
 {
     double si1,si2,be1,be2;
      // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle
-     // tested against limiting cases of cylinder, elliptical_cylinder and core_shell_bicelle
+     // tested against limiting cases of cylinder, elliptical_cylinder, stacked_discs, and core_shell_bicelle
      //    const double uplim = M_PI_4;
     const double halfheight = 0.5*length;
@@ -55 +33 @@
     //const double vbj=M_PI;
 
-    const double radius_major = rad * x_core;
-    const double rA = 0.5*(square(radius_major) + square(rad));
-    const double rB = 0.5*(square(radius_major) - square(rad));
-    const double dr1 = (rhoc-rhoh)   *M_PI*rad*radius_major*(2.0*halfheight);;
-    const double dr2 = (rhor-rhosolv)*M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick);
-    const double dr3 = (rhoh-rhor)   *M_PI*rad*radius_major*2.0*(halfheight+facthick);
-    //const double vol1 = M_PI*rad*radius_major*(2.0*halfheight);
-    //const double vol2 = M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick);
-    //const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick);
+    const double r_major = r_minor * x_core;
+    const double r2A = 0.5*(square(r_major) + square(r_minor));
+    const double r2B = 0.5*(square(r_major) - square(r_minor));
+    const double dr1 = (rhoc-rhoh)   *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+thick_face);
+    const double dr3 = (rhoh-rhor)   *M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
+    //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+thick_face);
+    //const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
 
     //initialize integral
@@ -74 +52 @@
         const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha);
         double inner_sum=0;
-        double sinarg1 = qq*halfheight*cos_alpha;
-        double sinarg2 = qq*(halfheight+facthick)*cos_alpha;
+        double sinarg1 = q*halfheight*cos_alpha;
+        double sinarg2 = q*(halfheight+thick_face)*cos_alpha;
         si1 = sas_sinx_x(sinarg1);
         si2 = sas_sinx_x(sinarg2);
@@ -82 +60 @@
             //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;
             const double beta = ( Gauss76Z[j] +1.0)*M_PI_2;
-            const double rr = sqrt(rA - rB*cos(beta));
-            double besarg1 = qq*rr*sin_alpha;
-            double besarg2 = qq*(rr+radthick)*sin_alpha;
+            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);
@@ -98 +76 @@
 }
 
-double 
+static double
 Iqxy(double qx, double qy,
-          double rad,
+          double r_minor,
           double x_core,
-          double radthick,
-          double facthick,
+          double thick_rim,
+          double thick_face,
           double length,
           double rhoc,
@@ -114 +92 @@
 {
        // THIS NEEDS TESTING
-    double qq, cos_val, cos_mu, cos_nu;
-    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, qq, cos_val, cos_mu, cos_nu);
+    double q, xhat, yhat, zhat;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
     const double dr1 = rhoc-rhoh;
     const double dr2 = rhor-rhosolv;
     const double dr3 = rhoh-rhor;
-    const double radius_major = rad*x_core;
+    const double r_major = r_minor*x_core;
     const double halfheight = 0.5*length;
-    const double vol1 = M_PI*rad*radius_major*length;
-    const double vol2 = M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick);
-    const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick);
+    const double vol1 = M_PI*r_minor*r_major*length;
+    const double vol2 = M_PI*(r_minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);
+    const double vol3 = M_PI*r_minor*r_major*2.0*(halfheight+thick_face);
 
-    // Compute:  r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2)
-    // Given:    radius_major = r_ratio * radius_minor  
-    // ASSUME the sin_alpha is included in the separate integration over orientation of rod angle
-    const double r = rad*sqrt(square(x_core*cos_nu) + cos_mu*cos_mu);
-    const double be1 = sas_2J1x_x( qq*r );
-    const double be2 = sas_2J1x_x( qq*(r + radthick ) );
-    const double si1 = sas_sinx_x( qq*halfheight*cos_val );
-    const double si2 = sas_sinx_x( qq*(halfheight + facthick)*cos_val );
+    // Compute effective radius in rotated coordinates
+    const double r_hat = sqrt(square(r_major*xhat) + square(r_minor*yhat));
+    const double rshell_hat = sqrt(square((r_major+thick_rim)*xhat)
+                                   + square((r_minor+thick_rim)*yhat));
+    const double be1 = sas_2J1x_x( q*r_hat );
+    const double be2 = sas_2J1x_x( q*rshell_hat );
+    const double si1 = sas_sinx_x( q*halfheight*zhat );
+    const double si2 = sas_sinx_x( q*(halfheight + thick_face)*zhat );
     const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 +  vol3*dr3*si2*be1);
-    //const double vol = form_volume(radius_minor, r_ratio, length);
     return 1.0e-4 * Aq;
 }