Changes in sasmodels/models/core_shell_bicelle_elliptical.c [dedcf34:592343f] 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!
-static double
-form_volume(double r_minor,
-        double x_core,
-        double thick_rim,
-        double thick_face,
-        double length)
+double form_volume(double radius, double x_core, double thick_rim, double thick_face, double length)
 {
-    return M_PI*(r_minor+thick_rim)*(r_minor*x_core+thick_rim)*(length+2.0*thick_face);
+    return M_PI*(radius+thick_rim)*(radius*x_core+thick_rim)*(length+2.0*thick_face);
 }
 
-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 
+                Iq(double qq,
+                   double rad,
+                   double x_core,
+                   double radthick,
+                   double facthick,
+                   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, stacked_discs, and core_shell_bicelle
+     // tested against limiting cases of cylinder, elliptical_cylinder and core_shell_bicelle
      //    const double uplim = M_PI_4;
     const double halfheight = 0.5*length;
@@ -33 +55 @@
     //const double vbj=M_PI;
 
-    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);
+    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);
 
     //initialize integral
@@ -52 +74 @@
         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;
+        double sinarg1 = qq*halfheight*cos_alpha;
+        double sinarg2 = qq*(halfheight+facthick)*cos_alpha;
         si1 = sas_sinx_x(sinarg1);
         si2 = sas_sinx_x(sinarg2);
@@ -60 +82 @@
             //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(r2A - r2B*cos(beta));
-            double besarg1 = q*rr*sin_alpha;
-            double besarg2 = q*(rr+thick_rim)*sin_alpha;
+            const double rr = sqrt(rA - rB*cos(beta));
+            double besarg1 = qq*rr*sin_alpha;
+            double besarg2 = qq*(rr+radthick)*sin_alpha;
             be1 = sas_2J1x_x(besarg1);
             be2 = sas_2J1x_x(besarg2);
@@ -76 +98 @@
 }
 
-static double
+double 
 Iqxy(double qx, double qy,
-          double r_minor,
+          double rad,
           double x_core,
-          double thick_rim,
-          double thick_face,
+          double radthick,
+          double facthick,
           double length,
           double rhoc,
@@ -92 +114 @@
 {
        // THIS NEEDS TESTING
-    double q, xhat, yhat, zhat;
-    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);
+    double qq, cos_val, cos_mu, cos_nu;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, qq, cos_val, cos_mu, cos_nu);
     const double dr1 = rhoc-rhoh;
     const double dr2 = rhor-rhosolv;
     const double dr3 = rhoh-rhor;
-    const double r_major = r_minor*x_core;
+    const double radius_major = rad*x_core;
     const double halfheight = 0.5*length;
-    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);
+    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);
 
-    // 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 );
+    // 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 );
     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;
 }