Changeset 82c11d3 in sasview for sansmodels/src/c_models

Timestamp:

Jan 5, 2012 6:24:51 PM (13 years ago)

Author:

Mathieu Doucet <doucetm@…>

Branches:

master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload

Children:

20d91bd

Parents:

98fdccd

Message:

refactored bunch of models

Location:

sansmodels/src/c_models

Files:

• binaryHS.cpp (modified) (5 diffs)
• binaryHS_PSF11.cpp (modified) (2 diffs)
• c_models.cpp (modified) (1 diff)
• coreshellcylinder.cpp (modified) (8 diffs)
• ellipsoid.cpp (modified) (8 diffs)
• ellipticalcylinder.cpp (modified) (5 diffs)
• flexcyl_ellipX.cpp (modified) (5 diffs)
• flexiblecylinder.cpp (modified) (1 diff)
• fractal.cpp (modified) (5 diffs)
• hollowcylinder.cpp (modified) (5 diffs)
• lamellar.cpp (modified) (2 diffs)
• lamellarFF_HG.cpp (modified) (5 diffs)
• lamellarPC.cpp (modified) (5 diffs)
• lamellarPS.cpp (modified) (5 diffs)
• lamellarPS_HG.cpp (modified) (2 diffs)
• lorentzian.cpp (added)
• models.hh (deleted)
• multishell.cpp (modified) (5 diffs)
• parallelepiped.cpp (modified) (1 diff)
• polygausscoil.cpp (modified) (1 diff)
• sc.cpp (modified) (1 diff)
• sphere.cpp (modified) (1 diff)
• spheroid.cpp (modified) (6 diffs)
• stackeddisks.cpp (modified) (5 diffs)
• triaxialellipsoid.cpp (modified) (5 diffs)
• vesicle.cpp (modified) (5 diffs)

sansmodels/src/c_models/binaryHS.cpp

@@ -21 +21 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "binaryHS.h"
 
 extern "C" {
-        #include "libSphere.h"
-        #include "binaryHS.h"
+#include "libSphere.h"
 }
 
 BinaryHSModel :: BinaryHSModel() {
 
-        l_radius     = Parameter(160.0, true);
-        l_radius.set_min(0.0);
-        s_radius    = Parameter(25.0, true);
-        s_radius.set_min(0.0);
-        vol_frac_ls  = Parameter(0.2);
-        vol_frac_ss  = Parameter(0.1);
-        ls_sld      = Parameter(3.5e-6);
-        ss_sld     = Parameter(5e-7);
-        solvent_sld   = Parameter(6.36e-6);
-        background = Parameter(0.0);
+  l_radius     = Parameter(160.0, true);
+  l_radius.set_min(0.0);
+  s_radius    = Parameter(25.0, true);
+  s_radius.set_min(0.0);
+  vol_frac_ls  = Parameter(0.2);
+  vol_frac_ss  = Parameter(0.1);
+  ls_sld      = Parameter(3.5e-6);
+  ss_sld     = Parameter(5e-7);
+  solvent_sld   = Parameter(6.36e-6);
+  background = Parameter(0.0);
 }
 
@@ -52 +51 @@
  */
 double BinaryHSModel :: operator()(double q) {
-        double dp[8];
+  double dp[8];
 
-        // Fill parameter array for IGOR library
-        // Add the background after averaging
-        dp[0] = l_radius();
-        dp[1] = s_radius();
-        dp[2] = vol_frac_ls();
-        dp[3] = vol_frac_ss();
-        dp[4] = ls_sld();
-        dp[5] = ss_sld();
-        dp[6] = solvent_sld();
-        dp[7] = 0.0;
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = l_radius();
+  dp[1] = s_radius();
+  dp[2] = vol_frac_ls();
+  dp[3] = vol_frac_ss();
+  dp[4] = ls_sld();
+  dp[5] = ss_sld();
+  dp[6] = solvent_sld();
+  dp[7] = 0.0;
 
 
-        // Get the dispersion points for the large radius
-        vector<WeightPoint> weights_l_radius;
-        l_radius.get_weights(weights_l_radius);
+  // Get the dispersion points for the large radius
+  vector<WeightPoint> weights_l_radius;
+  l_radius.get_weights(weights_l_radius);
 
-        // Get the dispersion points for the small radius
-        vector<WeightPoint> weights_s_radius;
-        s_radius.get_weights(weights_s_radius);
+  // Get the dispersion points for the small radius
+  vector<WeightPoint> weights_s_radius;
+  s_radius.get_weights(weights_s_radius);
 
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
 
-        // Loop over larger radius weight points
-        for(int i=0; i< (int)weights_l_radius.size(); i++) {
-                dp[0] = weights_l_radius[i].value;
+  // Loop over larger radius weight points
+  for(int i=0; i< (int)weights_l_radius.size(); i++) {
+    dp[0] = weights_l_radius[i].value;
 
-                // Loop over small radius weight points
-                for(int j=0; j< (int)weights_s_radius.size(); j++) {
-                        dp[1] = weights_s_radius[j].value;
+    // Loop over small radius weight points
+    for(int j=0; j< (int)weights_s_radius.size(); j++) {
+      dp[1] = weights_s_radius[j].value;
 
 
-                        sum += weights_l_radius[i].weight *weights_s_radius[j].weight * BinaryHS(dp, q);
-                        norm += weights_l_radius[i].weight *weights_s_radius[j].weight;
-                }
-        }
-        return sum/norm + background();
+      sum += weights_l_radius[i].weight *weights_s_radius[j].weight * BinaryHS(dp, q);
+      norm += weights_l_radius[i].weight *weights_s_radius[j].weight;
+    }
+  }
+  return sum/norm + background();
 }
 
@@ -101 +100 @@
  */
 double BinaryHSModel :: operator()(double qx, double qy) {
-        double q = sqrt(qx*qx + qy*qy);
-        return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -113 +112 @@
  */
 double BinaryHSModel :: evaluate_rphi(double q, double phi) {
-        return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 
@@ -121 +120 @@
  */
 double BinaryHSModel :: calculate_ER() {
-        //NOT implemented yet!!!
-        return 0.0;
+  //NOT implemented yet!!!
+  return 0.0;
 }
 

sansmodels/src/c_models/binaryHS_PSF11.cpp

@@ -21 +21 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
@@ -29 +28 @@
         #include "libSphere.h"
 }
+
+class BinaryHSPSF11Model{
+public:
+  // Model parameters
+  Parameter l_radius;
+  Parameter s_radius;
+  Parameter vol_frac_ls;
+  Parameter vol_frac_ss;
+  Parameter ls_sld;
+  Parameter ss_sld;
+  Parameter solvent_sld;
+  Parameter background;
+
+  //Constructor
+  BinaryHSPSF11Model();
+
+  //Operators to get I(Q)
+  double operator()(double q);
+  double operator()(double qx , double qy);
+  double calculate_ER();
+  double evaluate_rphi(double q, double phi);
+};
 
 BinaryHSPSF11Model :: BinaryHSPSF11Model() {

sansmodels/src/c_models/c_models.cpp

@@ -68 +68 @@
         void addDisperser(PyObject *module);
         void addCGaussian(PyObject *module);
-        void addCLorentzian(PyObject *module);
         void addCLogNormal(PyObject *module);
         void addCSchulz(PyObject *module);
 }
+void addCLorentzian(PyObject *module);
 
 /**

sansmodels/src/c_models/coreshellcylinder.cpp

@@ -27 +27 @@
 
 extern "C" {
-        #include "libCylinder.h"
-        #include "libStructureFactor.h"
+#include "libCylinder.h"
+#include "libStructureFactor.h"
 }
 
@@ -63 +63 @@
   double phi = pars->axis_phi * pi/180.0;
 
-    // Cylinder orientation
-    cyl_x = sin(theta) * cos(phi);
-    cyl_y = sin(theta) * sin(phi);
-    cyl_z = cos(theta);
-
-    // q vector
-    q_z = 0;
-
-    // Compute the angle btw vector q and the
-    // axis of the cylinder
-    cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
-
-    // The following test should always pass
-    if (fabs(cos_val)>1.0) {
-      printf("core_shell_cylinder_analytical_2D: Unexpected error: cos(alpha)=%g\n", cos_val);
-      return 0;
-    }
+  // Cylinder orientation
+  cyl_x = sin(theta) * cos(phi);
+  cyl_y = sin(theta) * sin(phi);
+  cyl_z = cos(theta);
+
+  // q vector
+  q_z = 0;
+
+  // Compute the angle btw vector q and the
+  // axis of the cylinder
+  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
+
+  // The following test should always pass
+  if (fabs(cos_val)>1.0) {
+    printf("core_shell_cylinder_analytical_2D: Unexpected error: cos(alpha)=%g\n", cos_val);
+    return 0;
+  }
 
   alpha = acos( cos_val );
@@ -85 +85 @@
   // Call the IGOR library function to get the kernel
   answer = CoreShellCylKernel(q, pars->radius, pars->thickness,
-                pars->core_sld,pars->shell_sld,
-                pars->solvent_sld, pars->length/2.0, alpha) / fabs(sin(alpha));
+      pars->core_sld,pars->shell_sld,
+      pars->solvent_sld, pars->length/2.0, alpha) / fabs(sin(alpha));
 
   //normalize by cylinder volume
   vol=pi*(pars->radius+pars->thickness)
-      *(pars->radius+pars->thickness)
-      *(pars->length+2.0*pars->thickness);
+          *(pars->radius+pars->thickness)
+          *(pars->length+2.0*pars->thickness);
   answer /= vol;
 
@@ -115 +115 @@
   double q;
   q = sqrt(qx*qx+qy*qy);
-    return core_shell_cylinder_analytical_2D_scaled(pars, q, qx/q, qy/q);
+  return core_shell_cylinder_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
 
 CoreShellCylinderModel :: CoreShellCylinderModel() {
-        scale      = Parameter(1.0);
-        radius     = Parameter(20.0, true);
-        radius.set_min(0.0);
-        thickness  = Parameter(10.0, true);
-        thickness.set_min(0.0);
-        length     = Parameter(400.0, true);
-        length.set_min(0.0);
-        core_sld   = Parameter(1.e-6);
-        shell_sld  = Parameter(4.e-6);
-        solvent_sld= Parameter(1.e-6);
-        background = Parameter(0.0);
-        axis_theta = Parameter(90.0, true);
-        axis_phi   = Parameter(0.0, true);
+  scale      = Parameter(1.0);
+  radius     = Parameter(20.0, true);
+  radius.set_min(0.0);
+  thickness  = Parameter(10.0, true);
+  thickness.set_min(0.0);
+  length     = Parameter(400.0, true);
+  length.set_min(0.0);
+  core_sld   = Parameter(1.e-6);
+  shell_sld  = Parameter(4.e-6);
+  solvent_sld= Parameter(1.e-6);
+  background = Parameter(0.0);
+  axis_theta = Parameter(90.0, true);
+  axis_phi   = Parameter(0.0, true);
 }
 
@@ -142 +142 @@
  */
 double CoreShellCylinderModel :: operator()(double q) {
-        double dp[8];
-
-        dp[0] = scale();
-        dp[1] = radius();
-        dp[2] = thickness();
-        dp[3] = length();
-        dp[4] = core_sld();
-        dp[5] = shell_sld();
-        dp[6] = solvent_sld();
-        dp[7] = 0.0;
-
-        // Get the dispersion points for the radius
-        vector<WeightPoint> weights_rad;
-        radius.get_weights(weights_rad);
-
-        // Get the dispersion points for the thickness
-        vector<WeightPoint> weights_thick;
-        thickness.get_weights(weights_thick);
-
-        // Get the dispersion points for the length
-        vector<WeightPoint> weights_len;
-        length.get_weights(weights_len);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double vol = 0.0;
-
-        // Loop over radius weight points
-        for(size_t i=0; i<weights_rad.size(); i++) {
-                dp[1] = weights_rad[i].value;
-
-                // Loop over length weight points
-                for(size_t j=0; j<weights_len.size(); j++) {
-                        dp[3] = weights_len[j].value;
-
-                        // Loop over thickness weight points
-                        for(size_t k=0; k<weights_thick.size(); k++) {
-                                dp[2] = weights_thick[k].value;
-                                //Un-normalize by volume
-                                sum += weights_rad[i].weight
-                                        * weights_len[j].weight
-                                        * weights_thick[k].weight
-                                        * CoreShellCylinder(dp, q)
-                                        * pow(weights_rad[i].value+weights_thick[k].value,2)
-                                        *(weights_len[j].value+2.0*weights_thick[k].value);
-                                //Find average volume
-                                vol += weights_rad[i].weight
-                                        * weights_len[j].weight
-                                        * weights_thick[k].weight
-                                        * pow(weights_rad[i].value+weights_thick[k].value,2)
-                                        *(weights_len[j].value+2.0*weights_thick[k].value);
-                                norm += weights_rad[i].weight
-                                * weights_len[j].weight
-                                * weights_thick[k].weight;
-                        }
-                }
-        }
-
-        if (vol != 0.0 && norm != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm);}
-
-        return sum/norm + background();
+  double dp[8];
+
+  dp[0] = scale();
+  dp[1] = radius();
+  dp[2] = thickness();
+  dp[3] = length();
+  dp[4] = core_sld();
+  dp[5] = shell_sld();
+  dp[6] = solvent_sld();
+  dp[7] = 0.0;
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad;
+  radius.get_weights(weights_rad);
+
+  // Get the dispersion points for the thickness
+  vector<WeightPoint> weights_thick;
+  thickness.get_weights(weights_thick);
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_len;
+  length.get_weights(weights_len);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp[1] = weights_rad[i].value;
+
+    // Loop over length weight points
+    for(size_t j=0; j<weights_len.size(); j++) {
+      dp[3] = weights_len[j].value;
+
+      // Loop over thickness weight points
+      for(size_t k=0; k<weights_thick.size(); k++) {
+        dp[2] = weights_thick[k].value;
+        //Un-normalize by volume
+        sum += weights_rad[i].weight
+            * weights_len[j].weight
+            * weights_thick[k].weight
+            * CoreShellCylinder(dp, q)
+        * pow(weights_rad[i].value+weights_thick[k].value,2)
+        *(weights_len[j].value+2.0*weights_thick[k].value);
+        //Find average volume
+        vol += weights_rad[i].weight
+            * weights_len[j].weight
+            * weights_thick[k].weight
+            * pow(weights_rad[i].value+weights_thick[k].value,2)
+        *(weights_len[j].value+2.0*weights_thick[k].value);
+        norm += weights_rad[i].weight
+            * weights_len[j].weight
+            * weights_thick[k].weight;
+      }
+    }
+  }
+
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -215 +215 @@
  */
 double CoreShellCylinderModel :: operator()(double qx, double qy) {
-        CoreShellCylinderParameters dp;
-        // Fill parameter array
-        dp.scale      = scale();
-        dp.radius     = radius();
-        dp.thickness  = thickness();
-        dp.length     = length();
-        dp.core_sld   = core_sld();
-        dp.shell_sld  = shell_sld();
-        dp.solvent_sld= solvent_sld();
-        dp.background = 0.0;
-        dp.axis_theta = axis_theta();
-        dp.axis_phi   = axis_phi();
-
-        // Get the dispersion points for the radius
-        vector<WeightPoint> weights_rad;
-        radius.get_weights(weights_rad);
-
-        // Get the dispersion points for the thickness
-        vector<WeightPoint> weights_thick;
-        thickness.get_weights(weights_thick);
-
-        // Get the dispersion points for the length
-        vector<WeightPoint> weights_len;
-        length.get_weights(weights_len);
-
-        // Get angular averaging for theta
-        vector<WeightPoint> weights_theta;
-        axis_theta.get_weights(weights_theta);
-
-        // Get angular averaging for phi
-        vector<WeightPoint> weights_phi;
-        axis_phi.get_weights(weights_phi);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double norm_vol = 0.0;
-        double vol = 0.0;
-        double pi = 4.0*atan(1.0);
-        // Loop over radius weight points
-        for(size_t i=0; i<weights_rad.size(); i++) {
-                dp.radius = weights_rad[i].value;
-
-
-                // Loop over length weight points
-                for(size_t j=0; j<weights_len.size(); j++) {
-                        dp.length = weights_len[j].value;
-
-                        // Loop over thickness weight points
-                        for(size_t m=0; m<weights_thick.size(); m++) {
-                                dp.thickness = weights_thick[m].value;
-
-                        // Average over theta distribution
-                        for(size_t k=0; k<weights_theta.size(); k++) {
-                                dp.axis_theta = weights_theta[k].value;
-
-                                // Average over phi distribution
-                                for(size_t l=0; l<weights_phi.size(); l++) {
-                                        dp.axis_phi = weights_phi[l].value;
-                                        //Un-normalize by volume
-                                        double _ptvalue = weights_rad[i].weight
-                                                * weights_len[j].weight
-                                                * weights_thick[m].weight
-                                                * weights_theta[k].weight
-                                                * weights_phi[l].weight
-                                                * core_shell_cylinder_analytical_2DXY(&dp, qx, qy)
-                                                * pow(weights_rad[i].value+weights_thick[m].value,2)
-                                                *(weights_len[j].value+2.0*weights_thick[m].value);
-
-                                        if (weights_theta.size()>1) {
-                                                _ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
-                                        }
-                                        sum += _ptvalue;
-
-                                        //Find average volume
-                                        vol += weights_rad[i].weight
-                                                        * weights_len[j].weight
-                                                        * weights_thick[m].weight
-                                                        * pow(weights_rad[i].value+weights_thick[m].value,2)
-                                                        *(weights_len[j].value+2.0*weights_thick[m].value);
-                                        //Find norm for volume
-                                        norm_vol += weights_rad[i].weight
-                                                                * weights_len[j].weight
-                                                                * weights_thick[m].weight;
-
-                                        norm += weights_rad[i].weight
-                                                * weights_len[j].weight
-                                                * weights_thick[m].weight
-                                                * weights_theta[k].weight
-                                                * weights_phi[l].weight;
-
-                                }
-                        }
-                        }
-                }
-        }
-        // Averaging in theta needs an extra normalization
-        // factor to account for the sin(theta) term in the
-        // integration (see documentation).
-        if (weights_theta.size()>1) norm = norm / asin(1.0);
-
-        if (vol != 0.0 && norm_vol != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm_vol);}
-
-        return sum/norm + background();
+  CoreShellCylinderParameters dp;
+  // Fill parameter array
+  dp.scale      = scale();
+  dp.radius     = radius();
+  dp.thickness  = thickness();
+  dp.length     = length();
+  dp.core_sld   = core_sld();
+  dp.shell_sld  = shell_sld();
+  dp.solvent_sld= solvent_sld();
+  dp.background = 0.0;
+  dp.axis_theta = axis_theta();
+  dp.axis_phi   = axis_phi();
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad;
+  radius.get_weights(weights_rad);
+
+  // Get the dispersion points for the thickness
+  vector<WeightPoint> weights_thick;
+  thickness.get_weights(weights_thick);
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_len;
+  length.get_weights(weights_len);
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_theta;
+  axis_theta.get_weights(weights_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_phi;
+  axis_phi.get_weights(weights_phi);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double norm_vol = 0.0;
+  double vol = 0.0;
+  double pi = 4.0*atan(1.0);
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp.radius = weights_rad[i].value;
+
+
+    // Loop over length weight points
+    for(size_t j=0; j<weights_len.size(); j++) {
+      dp.length = weights_len[j].value;
+
+      // Loop over thickness weight points
+      for(size_t m=0; m<weights_thick.size(); m++) {
+        dp.thickness = weights_thick[m].value;
+
+        // Average over theta distribution
+        for(size_t k=0; k<weights_theta.size(); k++) {
+          dp.axis_theta = weights_theta[k].value;
+
+          // Average over phi distribution
+          for(size_t l=0; l<weights_phi.size(); l++) {
+            dp.axis_phi = weights_phi[l].value;
+            //Un-normalize by volume
+            double _ptvalue = weights_rad[i].weight
+                * weights_len[j].weight
+                * weights_thick[m].weight
+                * weights_theta[k].weight
+                * weights_phi[l].weight
+                * core_shell_cylinder_analytical_2DXY(&dp, qx, qy)
+            * pow(weights_rad[i].value+weights_thick[m].value,2)
+            *(weights_len[j].value+2.0*weights_thick[m].value);
+
+            if (weights_theta.size()>1) {
+              _ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
+            }
+            sum += _ptvalue;
+
+            //Find average volume
+            vol += weights_rad[i].weight
+                * weights_len[j].weight
+                * weights_thick[m].weight
+                * pow(weights_rad[i].value+weights_thick[m].value,2)
+            *(weights_len[j].value+2.0*weights_thick[m].value);
+            //Find norm for volume
+            norm_vol += weights_rad[i].weight
+                * weights_len[j].weight
+                * weights_thick[m].weight;
+
+            norm += weights_rad[i].weight
+                * weights_len[j].weight
+                * weights_thick[m].weight
+                * weights_theta[k].weight
+                * weights_phi[l].weight;
+
+          }
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_theta.size()>1) norm = norm / asin(1.0);
+
+  if (vol != 0.0 && norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm_vol);}
+
+  return sum/norm + background();
 }
 
@@ -331 +331 @@
  */
 double CoreShellCylinderModel :: evaluate_rphi(double q, double phi) {
-        double qx = q*cos(phi);
-        double qy = q*sin(phi);
-        return (*this).operator()(qx, qy);
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
 }
 /**
@@ -340 +340 @@
  */
 double CoreShellCylinderModel :: calculate_ER() {
-        CoreShellCylinderParameters dp;
-
-        dp.radius     = radius();
-        dp.thickness  = thickness();
-        dp.length     = length();
-        double rad_out = 0.0;
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-
-        // Get the dispersion points for the length
-        vector<WeightPoint> weights_length;
-        length.get_weights(weights_length);
-
-        // Get the dispersion points for the thickness
-        vector<WeightPoint> weights_thickness;
-        thickness.get_weights(weights_thickness);
-
-        // Get the dispersion points for the radius
-        vector<WeightPoint> weights_radius ;
-        radius.get_weights(weights_radius);
-
-        // Loop over major shell weight points
-        for(int i=0; i< (int)weights_length.size(); i++) {
-                dp.length = weights_length[i].value;
-                for(int j=0; j< (int)weights_thickness.size(); j++) {
-                        dp.thickness = weights_thickness[j].value;
-                        for(int k=0; k< (int)weights_radius.size(); k++) {
-                                dp.radius = weights_radius[k].value;
-                                //Note: output of "DiamCyl( )" is DIAMETER.
-                                sum +=weights_length[i].weight * weights_thickness[j].weight
-                                        * weights_radius[k].weight*DiamCyl(dp.length+2.0*dp.thickness,dp.radius+dp.thickness)/2.0;
-                                norm += weights_length[i].weight* weights_thickness[j].weight* weights_radius[k].weight;
-                        }
-                }
-        }
-        if (norm != 0){
-                //return the averaged value
-                rad_out =  sum/norm;}
-        else{
-                //return normal value
-                //Note: output of "DiamCyl()" is DIAMETER.
-                rad_out = DiamCyl(dp.length+2.0*dp.thickness,dp.radius+dp.thickness)/2.0;}
-
-        return rad_out;
-}
+  CoreShellCylinderParameters dp;
+
+  dp.radius     = radius();
+  dp.thickness  = thickness();
+  dp.length     = length();
+  double rad_out = 0.0;
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_length;
+  length.get_weights(weights_length);
+
+  // Get the dispersion points for the thickness
+  vector<WeightPoint> weights_thickness;
+  thickness.get_weights(weights_thickness);
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_radius ;
+  radius.get_weights(weights_radius);
+
+  // Loop over major shell weight points
+  for(int i=0; i< (int)weights_length.size(); i++) {
+    dp.length = weights_length[i].value;
+    for(int j=0; j< (int)weights_thickness.size(); j++) {
+      dp.thickness = weights_thickness[j].value;
+      for(int k=0; k< (int)weights_radius.size(); k++) {
+        dp.radius = weights_radius[k].value;
+        //Note: output of "DiamCyl( )" is DIAMETER.
+        sum +=weights_length[i].weight * weights_thickness[j].weight
+            * weights_radius[k].weight*DiamCyl(dp.length+2.0*dp.thickness,dp.radius+dp.thickness)/2.0;
+        norm += weights_length[i].weight* weights_thickness[j].weight* weights_radius[k].weight;
+      }
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    //Note: output of "DiamCyl()" is DIAMETER.
+    rad_out = DiamCyl(dp.length+2.0*dp.thickness,dp.radius+dp.thickness)/2.0;}
+
+  return rad_out;
+}

sansmodels/src/c_models/ellipsoid.cpp

@@ -27 +27 @@
 
 extern "C" {
-        #include "libCylinder.h"
-        #include "libStructureFactor.h"
+#include "libCylinder.h"
+#include "libStructureFactor.h"
 }
 
@@ -60 +60 @@
   double phi = pars->axis_phi * pi/180.0;
 
-    // Ellipsoid orientation
-    cyl_x = sin(theta) * cos(phi);
-    cyl_y = sin(theta) * sin(phi);
-    cyl_z = cos(theta);
-
-    // q vector
-    q_z = 0.0;
-
-    // Compute the angle btw vector q and the
-    // axis of the cylinder
-    cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
-
-    // The following test should always pass
-    if (fabs(cos_val)>1.0) {
-      printf("ellipsoid_ana_2D: Unexpected error: cos(alpha)>1\n");
-      return 0;
-    }
-
-    // Angle between rotation axis and q vector
+  // Ellipsoid orientation
+  cyl_x = sin(theta) * cos(phi);
+  cyl_y = sin(theta) * sin(phi);
+  cyl_z = cos(theta);
+
+  // q vector
+  q_z = 0.0;
+
+  // Compute the angle btw vector q and the
+  // axis of the cylinder
+  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
+
+  // The following test should always pass
+  if (fabs(cos_val)>1.0) {
+    printf("ellipsoid_ana_2D: Unexpected error: cos(alpha)>1\n");
+    return 0;
+  }
+
+  // Angle between rotation axis and q vector
   alpha = acos( cos_val );
 
@@ -88 +88 @@
 
   //normalize by cylinder volume
-    vol = 4.0/3.0 * acos(-1.0) * pars->radius_b * pars->radius_b * pars->radius_a;
+  vol = 4.0/3.0 * acos(-1.0) * pars->radius_b * pars->radius_b * pars->radius_a;
   answer *= vol;
 
@@ -112 +112 @@
   double q;
   q = sqrt(qx*qx+qy*qy);
-    return ellipsoid_analytical_2D_scaled(pars, q, qx/q, qy/q);
+  return ellipsoid_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
 EllipsoidModel :: EllipsoidModel() {
-        scale      = Parameter(1.0);
-        radius_a   = Parameter(20.0, true);
-        radius_a.set_min(0.0);
-        radius_b   = Parameter(400.0, true);
-        radius_b.set_min(0.0);
-        sldEll   = Parameter(4.e-6);
-        sldSolv   = Parameter(1.e-6);
-        background = Parameter(0.0);
-        axis_theta  = Parameter(57.325, true);
-        axis_phi    = Parameter(0.0, true);
+  scale      = Parameter(1.0);
+  radius_a   = Parameter(20.0, true);
+  radius_a.set_min(0.0);
+  radius_b   = Parameter(400.0, true);
+  radius_b.set_min(0.0);
+  sldEll   = Parameter(4.e-6);
+  sldSolv   = Parameter(1.e-6);
+  background = Parameter(0.0);
+  axis_theta  = Parameter(57.325, true);
+  axis_phi    = Parameter(0.0, true);
 }
 
@@ -135 +135 @@
  */
 double EllipsoidModel :: operator()(double q) {
-        double dp[6];
-
-        // Fill parameter array for IGOR library
-        // Add the background after averaging
-        dp[0] = scale();
-        dp[1] = radius_a();
-        dp[2] = radius_b();
-        dp[3] = sldEll();
-        dp[4] = sldSolv();
-        dp[5] = 0.0;
-
-        // Get the dispersion points for the radius_a
-        vector<WeightPoint> weights_rad_a;
-        radius_a.get_weights(weights_rad_a);
-
-        // Get the dispersion points for the radius_b
-        vector<WeightPoint> weights_rad_b;
-        radius_b.get_weights(weights_rad_b);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double vol = 0.0;
-
-        // Loop over radius_a weight points
-        for(size_t i=0; i<weights_rad_a.size(); i++) {
-                dp[1] = weights_rad_a[i].value;
-
-                // Loop over radius_b weight points
-                for(size_t j=0; j<weights_rad_b.size(); j++) {
-                        dp[2] = weights_rad_b[j].value;
-                        //Un-normalize  by volume
-                        sum += weights_rad_a[i].weight
-                                * weights_rad_b[j].weight * EllipsoidForm(dp, q)
-                                * pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
-
-                        //Find average volume
-                        vol += weights_rad_a[i].weight
-                                * weights_rad_b[j].weight
-                                * pow(weights_rad_b[j].value,2)
-                                * weights_rad_a[i].value;
-                        norm += weights_rad_a[i].weight
-                                * weights_rad_b[j].weight;
-                }
-        }
-
-        if (vol != 0.0 && norm != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm);}
-
-        return sum/norm + background();
+  double dp[6];
+
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = radius_a();
+  dp[2] = radius_b();
+  dp[3] = sldEll();
+  dp[4] = sldSolv();
+  dp[5] = 0.0;
+
+  // Get the dispersion points for the radius_a
+  vector<WeightPoint> weights_rad_a;
+  radius_a.get_weights(weights_rad_a);
+
+  // Get the dispersion points for the radius_b
+  vector<WeightPoint> weights_rad_b;
+  radius_b.get_weights(weights_rad_b);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over radius_a weight points
+  for(size_t i=0; i<weights_rad_a.size(); i++) {
+    dp[1] = weights_rad_a[i].value;
+
+    // Loop over radius_b weight points
+    for(size_t j=0; j<weights_rad_b.size(); j++) {
+      dp[2] = weights_rad_b[j].value;
+      //Un-normalize  by volume
+      sum += weights_rad_a[i].weight
+          * weights_rad_b[j].weight * EllipsoidForm(dp, q)
+      * pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
+
+      //Find average volume
+      vol += weights_rad_a[i].weight
+          * weights_rad_b[j].weight
+          * pow(weights_rad_b[j].value,2)
+      * weights_rad_a[i].value;
+      norm += weights_rad_a[i].weight
+          * weights_rad_b[j].weight;
+    }
+  }
+
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -195 +195 @@
  */
 double EllipsoidModel :: operator()(double qx, double qy) {
-        EllipsoidParameters dp;
-        // Fill parameter array
-        dp.scale      = scale();
-        dp.radius_a   = radius_a();
-        dp.radius_b   = radius_b();
-        dp.sldEll   = sldEll();
-        dp.sldSolv   = sldSolv();
-        dp.background = 0.0;
-        dp.axis_theta = axis_theta();
-        dp.axis_phi   = axis_phi();
-
-        // Get the dispersion points for the radius_a
-        vector<WeightPoint> weights_rad_a;
-        radius_a.get_weights(weights_rad_a);
-
-        // Get the dispersion points for the radius_b
-        vector<WeightPoint> weights_rad_b;
-        radius_b.get_weights(weights_rad_b);
-
-        // Get angular averaging for theta
-        vector<WeightPoint> weights_theta;
-        axis_theta.get_weights(weights_theta);
-
-        // Get angular averaging for phi
-        vector<WeightPoint> weights_phi;
-        axis_phi.get_weights(weights_phi);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double norm_vol = 0.0;
-        double vol = 0.0;
-        double pi = 4.0*atan(1.0);
-        // Loop over radius weight points
-        for(size_t i=0; i<weights_rad_a.size(); i++) {
-                dp.radius_a = weights_rad_a[i].value;
-
-
-                // Loop over length weight points
-                for(size_t j=0; j<weights_rad_b.size(); j++) {
-                        dp.radius_b = weights_rad_b[j].value;
-
-                        // Average over theta distribution
-                        for(size_t k=0; k<weights_theta.size(); k++) {
-                                dp.axis_theta = weights_theta[k].value;
-
-                                // Average over phi distribution
-                                for(size_t l=0; l<weights_phi.size(); l++) {
-                                        dp.axis_phi = weights_phi[l].value;
-                                        //Un-normalize by volume
-                                        double _ptvalue = weights_rad_a[i].weight
-                                                * weights_rad_b[j].weight
-                                                * weights_theta[k].weight
-                                                * weights_phi[l].weight
-                                                * ellipsoid_analytical_2DXY(&dp, qx, qy)
-                                                * pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
-                                        if (weights_theta.size()>1) {
-                                                _ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
-                                        }
-                                        sum += _ptvalue;
-                                        //Find average volume
-                                        vol += weights_rad_a[i].weight
-                                                * weights_rad_b[j].weight
-                                                * pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
-                                        //Find norm for volume
-                                        norm_vol += weights_rad_a[i].weight
-                                                * weights_rad_b[j].weight;
-
-                                        norm += weights_rad_a[i].weight
-                                                * weights_rad_b[j].weight
-                                                * weights_theta[k].weight
-                                                * weights_phi[l].weight;
-
-                                }
-                        }
-                }
-        }
-        // Averaging in theta needs an extra normalization
-        // factor to account for the sin(theta) term in the
-        // integration (see documentation).
-        if (weights_theta.size()>1) norm = norm / asin(1.0);
-
-        if (vol != 0.0 && norm_vol != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm_vol);}
-
-        return sum/norm + background();
+  EllipsoidParameters dp;
+  // Fill parameter array
+  dp.scale      = scale();
+  dp.radius_a   = radius_a();
+  dp.radius_b   = radius_b();
+  dp.sldEll   = sldEll();
+  dp.sldSolv   = sldSolv();
+  dp.background = 0.0;
+  dp.axis_theta = axis_theta();
+  dp.axis_phi   = axis_phi();
+
+  // Get the dispersion points for the radius_a
+  vector<WeightPoint> weights_rad_a;
+  radius_a.get_weights(weights_rad_a);
+
+  // Get the dispersion points for the radius_b
+  vector<WeightPoint> weights_rad_b;
+  radius_b.get_weights(weights_rad_b);
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_theta;
+  axis_theta.get_weights(weights_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_phi;
+  axis_phi.get_weights(weights_phi);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double norm_vol = 0.0;
+  double vol = 0.0;
+  double pi = 4.0*atan(1.0);
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad_a.size(); i++) {
+    dp.radius_a = weights_rad_a[i].value;
+
+
+    // Loop over length weight points
+    for(size_t j=0; j<weights_rad_b.size(); j++) {
+      dp.radius_b = weights_rad_b[j].value;
+
+      // Average over theta distribution
+      for(size_t k=0; k<weights_theta.size(); k++) {
+        dp.axis_theta = weights_theta[k].value;
+
+        // Average over phi distribution
+        for(size_t l=0; l<weights_phi.size(); l++) {
+          dp.axis_phi = weights_phi[l].value;
+          //Un-normalize by volume
+          double _ptvalue = weights_rad_a[i].weight
+              * weights_rad_b[j].weight
+              * weights_theta[k].weight
+              * weights_phi[l].weight
+              * ellipsoid_analytical_2DXY(&dp, qx, qy)
+          * pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
+          if (weights_theta.size()>1) {
+            _ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
+          }
+          sum += _ptvalue;
+          //Find average volume
+          vol += weights_rad_a[i].weight
+              * weights_rad_b[j].weight
+              * pow(weights_rad_b[j].value,2) * weights_rad_a[i].value;
+          //Find norm for volume
+          norm_vol += weights_rad_a[i].weight
+              * weights_rad_b[j].weight;
+
+          norm += weights_rad_a[i].weight
+              * weights_rad_b[j].weight
+              * weights_theta[k].weight
+              * weights_phi[l].weight;
+
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_theta.size()>1) norm = norm / asin(1.0);
+
+  if (vol != 0.0 && norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm_vol);}
+
+  return sum/norm + background();
 }
 
@@ -292 +292 @@
  */
 double EllipsoidModel :: evaluate_rphi(double q, double phi) {
-        double qx = q*cos(phi);
-        double qy = q*sin(phi);
-        return (*this).operator()(qx, qy);
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
 }
 
@@ -302 +302 @@
  */
 double EllipsoidModel :: calculate_ER() {
-        EllipsoidParameters dp;
-
-        dp.radius_a = radius_a();
-        dp.radius_b = radius_b();
-
-        double rad_out = 0.0;
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-
-        // Get the dispersion points for the major shell
-        vector<WeightPoint> weights_radius_a;
-        radius_a.get_weights(weights_radius_a);
-
-        // Get the dispersion points for the minor shell
-        vector<WeightPoint> weights_radius_b;
-        radius_b.get_weights(weights_radius_b);
-
-        // Loop over major shell weight points
-        for(int i=0; i< (int)weights_radius_b.size(); i++) {
-                dp.radius_b = weights_radius_b[i].value;
-                for(int k=0; k< (int)weights_radius_a.size(); k++) {
-                        dp.radius_a = weights_radius_a[k].value;
-                        sum +=weights_radius_b[i].weight
-                                * weights_radius_a[k].weight*DiamEllip(dp.radius_a,dp.radius_b)/2.0;
-                        norm += weights_radius_b[i].weight* weights_radius_a[k].weight;
-                }
-        }
-        if (norm != 0){
-                //return the averaged value
-                rad_out =  sum/norm;}
-        else{
-                //return normal value
-                rad_out = DiamEllip(dp.radius_a,dp.radius_b)/2.0;}
-
-        return rad_out;
-}
+  EllipsoidParameters dp;
+
+  dp.radius_a = radius_a();
+  dp.radius_b = radius_b();
+
+  double rad_out = 0.0;
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the major shell
+  vector<WeightPoint> weights_radius_a;
+  radius_a.get_weights(weights_radius_a);
+
+  // Get the dispersion points for the minor shell
+  vector<WeightPoint> weights_radius_b;
+  radius_b.get_weights(weights_radius_b);
+
+  // Loop over major shell weight points
+  for(int i=0; i< (int)weights_radius_b.size(); i++) {
+    dp.radius_b = weights_radius_b[i].value;
+    for(int k=0; k< (int)weights_radius_a.size(); k++) {
+      dp.radius_a = weights_radius_a[k].value;
+      sum +=weights_radius_b[i].weight
+          * weights_radius_a[k].weight*DiamEllip(dp.radius_a,dp.radius_b)/2.0;
+      norm += weights_radius_b[i].weight* weights_radius_a[k].weight;
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    rad_out = DiamEllip(dp.radius_a,dp.radius_b)/2.0;}
+
+  return rad_out;
+}

sansmodels/src/c_models/ellipticalcylinder.cpp

@@ -18 +18 @@
  *   sansmodels/src/libigor
  *
- *      TODO: refactor so that we pull in the old sansmodels.c_extensions
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
+#include <stdlib.h>
 using namespace std;
+#include "elliptical_cylinder.h"
 
 extern "C" {
-        #include "libCylinder.h"
-        #include "libStructureFactor.h"
-        #include "elliptical_cylinder.h"
+#include "libCylinder.h"
+#include "libStructureFactor.h"
+}
+
+typedef struct {
+  double scale;
+  double r_minor;
+  double r_ratio;
+  double length;
+  double sldCyl;
+  double sldSolv;
+  double background;
+  double cyl_theta;
+  double cyl_phi;
+  double cyl_psi;
+} EllipticalCylinderParameters;
+
+
+static double elliptical_cylinder_kernel(EllipticalCylinderParameters *pars, double q, double alpha, double nu) {
+  double qr;
+  double qL;
+  double Be,Si;
+  double r_major;
+  double kernel;
+
+  r_major = pars->r_ratio * pars->r_minor;
+
+  qr = q*sin(alpha)*sqrt( r_major*r_major*sin(nu)*sin(nu) + pars->r_minor*pars->r_minor*cos(nu)*cos(nu) );
+  qL = q*pars->length*cos(alpha)/2.0;
+
+  if (qr==0){
+    Be = 0.5;
+  }else{
+    Be = NR_BessJ1(qr)/qr;
+  }
+  if (qL==0){
+    Si = 1.0;
+  }else{
+    Si = sin(qL)/qL;
+  }
+
+
+  kernel = 2.0*Be * Si;
+  return kernel*kernel;
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double elliptical_cylinder_analytical_2D_scaled(EllipticalCylinderParameters *pars, double q, double q_x, double q_y) {
+  double cyl_x, cyl_y, cyl_z;
+  double ell_x, ell_y;
+  double q_z;
+  double alpha, vol, cos_val;
+  double nu, cos_nu;
+  double answer;
+  //convert angle degree to radian
+  double pi = 4.0*atan(1.0);
+  double theta = pars->cyl_theta * pi/180.0;
+  double phi = pars->cyl_phi * pi/180.0;
+  double psi = pars->cyl_psi * pi/180.0;
+
+  //Cylinder orientation
+  cyl_x = sin(theta) * cos(phi);
+  cyl_y = sin(theta) * sin(phi);
+  cyl_z = cos(theta);
+
+  // q vector
+  q_z = 0;
+
+  // Compute the angle btw vector q and the
+  // axis of the cylinder
+  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
+
+  // The following test should always pass
+  if (fabs(cos_val)>1.0) {
+    printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
+    return 0;
+  }
+
+  // Note: cos(alpha) = 0 and 1 will get an
+  // undefined value from CylKernel
+  alpha = acos( cos_val );
+
+  //ellipse orientation:
+  // the elliptical corss section was transformed and projected
+  // into the detector plane already through sin(alpha)and furthermore psi remains as same
+  // on the detector plane.
+  // So, all we need is to calculate the angle (nu) of the minor axis of the ellipse wrt
+  // the wave vector q.
+
+  //x- y- component on the detector plane.
+  ell_x =  cos(psi);
+  ell_y =  sin(psi);
+
+  // calculate the axis of the ellipse wrt q-coord.
+  cos_nu = ell_x*q_x + ell_y*q_y;
+  nu = acos(cos_nu);
+
+  // The following test should always pass
+  if (fabs(cos_nu)>1.0) {
+    printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n");
+    return 0;
+  }
+
+  answer = elliptical_cylinder_kernel(pars, q, alpha,nu);
+
+  // Multiply by contrast^2
+  answer *= (pars->sldCyl - pars->sldSolv) * (pars->sldCyl - pars->sldSolv);
+
+  //normalize by cylinder volume
+  //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
+  vol = acos(-1.0) * pars->r_minor * pars->r_minor * pars->r_ratio * pars->length;
+  answer *= vol;
+
+  //convert to [cm-1]
+  answer *= 1.0e8;
+
+  //Scale
+  answer *= pars->scale;
+
+  // add in the background
+  answer += pars->background;
+
+  return answer;
+}
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @return: function value
+ */
+static double elliptical_cylinder_analytical_2DXY(EllipticalCylinderParameters *pars, double qx, double qy) {
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return elliptical_cylinder_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
 EllipticalCylinderModel :: EllipticalCylinderModel() {
-        scale      = Parameter(1.0);
-        r_minor    = Parameter(20.0, true);
-        r_minor.set_min(0.0);
-        r_ratio    = Parameter(1.5, true);
-        r_ratio.set_min(0.0);
-        length     = Parameter(400.0, true);
-        length.set_min(0.0);
-        sldCyl   = Parameter(4.e-6);
-        sldSolv   = Parameter(1.e-6);
-        background = Parameter(0.0);
-        cyl_theta  = Parameter(57.325, true);
-        cyl_phi    = Parameter(0.0, true);
-        cyl_psi    = Parameter(0.0, true);
+  scale      = Parameter(1.0);
+  r_minor    = Parameter(20.0, true);
+  r_minor.set_min(0.0);
+  r_ratio    = Parameter(1.5, true);
+  r_ratio.set_min(0.0);
+  length     = Parameter(400.0, true);
+  length.set_min(0.0);
+  sldCyl   = Parameter(4.e-6);
+  sldSolv   = Parameter(1.e-6);
+  background = Parameter(0.0);
+  cyl_theta  = Parameter(57.325, true);
+  cyl_phi    = Parameter(0.0, true);
+  cyl_psi    = Parameter(0.0, true);
 }
 
@@ -56 +196 @@
  */
 double EllipticalCylinderModel :: operator()(double q) {
-        double dp[7];
-
-        dp[0] = scale();
-        dp[1] = r_minor();
-        dp[2] = r_ratio();
-        dp[3] = length();
-        dp[4] = sldCyl();
-        dp[5] = sldSolv();
-        dp[6] = 0.0;
-
-        // Get the dispersion points for the r_minor
-        vector<WeightPoint> weights_rad;
-        r_minor.get_weights(weights_rad);
-
-        // Get the dispersion points for the r_ratio
-        vector<WeightPoint> weights_rat;
-        r_ratio.get_weights(weights_rat);
-
-        // Get the dispersion points for the length
-        vector<WeightPoint> weights_len;
-        length.get_weights(weights_len);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double vol = 0.0;
-
-        // Loop over r_minor weight points
-        for(size_t i=0; i<weights_rad.size(); i++) {
-                dp[1] = weights_rad[i].value;
-
-                // Loop over r_ratio weight points
-                for(size_t j=0; j<weights_rat.size(); j++) {
-                        dp[2] = weights_rat[j].value;
-
-                        // Loop over length weight points
-                        for(size_t k=0; k<weights_len.size(); k++) {
-                                dp[3] = weights_len[k].value;
-                                //Un-normalize  by volume
-                                sum += weights_rad[i].weight
-                                        * weights_len[k].weight
-                                        * weights_rat[j].weight
-                                        * EllipCyl20(dp, q)
-                                        * pow(weights_rad[i].value,2) * weights_rat[j].value
-                                        * weights_len[k].value;
-                                //Find average volume
-                                vol += weights_rad[i].weight
-                                        * weights_len[k].weight
-                                        * weights_rat[j].weight
-                                        * pow(weights_rad[i].value,2) * weights_rat[j].value
-                                        * weights_len[k].value;
-                                norm += weights_rad[i].weight
-                                        * weights_len[k].weight
-                                        * weights_rat[j].weight;
-                        }
-                }
-        }
-
-        if (vol != 0.0 && norm != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm);}
-
-        return sum/norm + background();
+  double dp[7];
+
+  dp[0] = scale();
+  dp[1] = r_minor();
+  dp[2] = r_ratio();
+  dp[3] = length();
+  dp[4] = sldCyl();
+  dp[5] = sldSolv();
+  dp[6] = 0.0;
+
+  // Get the dispersion points for the r_minor
+  vector<WeightPoint> weights_rad;
+  r_minor.get_weights(weights_rad);
+
+  // Get the dispersion points for the r_ratio
+  vector<WeightPoint> weights_rat;
+  r_ratio.get_weights(weights_rat);
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_len;
+  length.get_weights(weights_len);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over r_minor weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp[1] = weights_rad[i].value;
+
+    // Loop over r_ratio weight points
+    for(size_t j=0; j<weights_rat.size(); j++) {
+      dp[2] = weights_rat[j].value;
+
+      // Loop over length weight points
+      for(size_t k=0; k<weights_len.size(); k++) {
+        dp[3] = weights_len[k].value;
+        //Un-normalize  by volume
+        sum += weights_rad[i].weight
+            * weights_len[k].weight
+            * weights_rat[j].weight
+            * EllipCyl20(dp, q)
+        * pow(weights_rad[i].value,2) * weights_rat[j].value
+        * weights_len[k].value;
+        //Find average volume
+        vol += weights_rad[i].weight
+            * weights_len[k].weight
+            * weights_rat[j].weight
+            * pow(weights_rad[i].value,2) * weights_rat[j].value
+            * weights_len[k].value;
+        norm += weights_rad[i].weight
+            * weights_len[k].weight
+            * weights_rat[j].weight;
+      }
+    }
+  }
+
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -128 +268 @@
  */
 double EllipticalCylinderModel :: operator()(double qx, double qy) {
-        EllipticalCylinderParameters dp;
-        // Fill parameter array
-        dp.scale      = scale();
-        dp.r_minor    = r_minor();
-        dp.r_ratio    = r_ratio();
-        dp.length     = length();
-        dp.sldCyl   = sldCyl();
-        dp.sldSolv   = sldSolv();
-        dp.background = 0.0;
-        dp.cyl_theta  = cyl_theta();
-        dp.cyl_phi    = cyl_phi();
-        dp.cyl_psi    = cyl_psi();
-
-        // Get the dispersion points for the r_minor
-        vector<WeightPoint> weights_rad;
-        r_minor.get_weights(weights_rad);
-
-        // Get the dispersion points for the r_ratio
-        vector<WeightPoint> weights_rat;
-        r_ratio.get_weights(weights_rat);
-
-        // Get the dispersion points for the length
-        vector<WeightPoint> weights_len;
-        length.get_weights(weights_len);
-
-        // Get angular averaging for theta
-        vector<WeightPoint> weights_theta;
-        cyl_theta.get_weights(weights_theta);
-
-        // Get angular averaging for phi
-        vector<WeightPoint> weights_phi;
-        cyl_phi.get_weights(weights_phi);
-
-        // Get angular averaging for psi
-        vector<WeightPoint> weights_psi;
-        cyl_psi.get_weights(weights_psi);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double norm_vol = 0.0;
-        double vol = 0.0;
-        double pi = 4.0*atan(1.0);
-        // Loop over minor radius weight points
-        for(size_t i=0; i<weights_rad.size(); i++) {
-                dp.r_minor = weights_rad[i].value;
-
-
-                // Loop over length weight points
-                for(size_t j=0; j<weights_len.size(); j++) {
-                        dp.length = weights_len[j].value;
-
-                        // Loop over r_ration weight points
-                        for(size_t m=0; m<weights_rat.size(); m++) {
-                                dp.r_ratio = weights_rat[m].value;
-
-                        // Average over theta distribution
-                        for(size_t k=0; k<weights_theta.size(); k++) {
-                                dp.cyl_theta = weights_theta[k].value;
-
-                                // Average over phi distribution
-                                for(size_t l=0; l<weights_phi.size(); l++) {
-                                        dp.cyl_phi = weights_phi[l].value;
-
-                                // Average over phi distribution
-                                for(size_t o=0; o<weights_psi.size(); o++) {
-                                        dp.cyl_psi = weights_psi[o].value;
-                                        //Un-normalize by volume
-                                        double _ptvalue = weights_rad[i].weight
-                                                * weights_len[j].weight
-                                                * weights_rat[m].weight
-                                                * weights_theta[k].weight
-                                                * weights_phi[l].weight
-                                                * weights_psi[o].weight
-                                                * elliptical_cylinder_analytical_2DXY(&dp, qx, qy)
-                                                * pow(weights_rad[i].value,2)
-                                                * weights_len[j].value
-                                                * weights_rat[m].value;
-                                        if (weights_theta.size()>1) {
-                                                _ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
-                                        }
-                                        sum += _ptvalue;
-                                        //Find average volume
-                                        vol += weights_rad[i].weight
-                                                * weights_len[j].weight
-                                                * weights_rat[m].weight
-                                                * pow(weights_rad[i].value,2)
-                                                * weights_len[j].value
-                                                * weights_rat[m].value;
-                                        //Find norm for volume
-                                        norm_vol += weights_rad[i].weight
-                                                * weights_len[j].weight
-                                                * weights_rat[m].weight;
-
-                                        norm += weights_rad[i].weight
-                                                * weights_len[j].weight
-                                                * weights_rat[m].weight
-                                                * weights_theta[k].weight
-                                                * weights_phi[l].weight
-                                                * weights_psi[o].weight;
-
-                                }
-                                }
-                        }
-                        }
-                }
-        }
-        // Averaging in theta needs an extra normalization
-        // factor to account for the sin(theta) term in the
-        // integration (see documentation).
-        if (weights_theta.size()>1) norm = norm / asin(1.0);
-
-        if (vol != 0.0 && norm_vol != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm_vol);}
-
-        return sum/norm + background();
+  EllipticalCylinderParameters dp;
+  // Fill parameter array
+  dp.scale      = scale();
+  dp.r_minor    = r_minor();
+  dp.r_ratio    = r_ratio();
+  dp.length     = length();
+  dp.sldCyl   = sldCyl();
+  dp.sldSolv   = sldSolv();
+  dp.background = 0.0;
+  dp.cyl_theta  = cyl_theta();
+  dp.cyl_phi    = cyl_phi();
+  dp.cyl_psi    = cyl_psi();
+
+  // Get the dispersion points for the r_minor
+  vector<WeightPoint> weights_rad;
+  r_minor.get_weights(weights_rad);
+
+  // Get the dispersion points for the r_ratio
+  vector<WeightPoint> weights_rat;
+  r_ratio.get_weights(weights_rat);
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_len;
+  length.get_weights(weights_len);
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_theta;
+  cyl_theta.get_weights(weights_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_phi;
+  cyl_phi.get_weights(weights_phi);
+
+  // Get angular averaging for psi
+  vector<WeightPoint> weights_psi;
+  cyl_psi.get_weights(weights_psi);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double norm_vol = 0.0;
+  double vol = 0.0;
+  double pi = 4.0*atan(1.0);
+  // Loop over minor radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp.r_minor = weights_rad[i].value;
+
+
+    // Loop over length weight points
+    for(size_t j=0; j<weights_len.size(); j++) {
+      dp.length = weights_len[j].value;
+
+      // Loop over r_ration weight points
+      for(size_t m=0; m<weights_rat.size(); m++) {
+        dp.r_ratio = weights_rat[m].value;
+
+        // Average over theta distribution
+        for(size_t k=0; k<weights_theta.size(); k++) {
+          dp.cyl_theta = weights_theta[k].value;
+
+          // Average over phi distribution
+          for(size_t l=0; l<weights_phi.size(); l++) {
+            dp.cyl_phi = weights_phi[l].value;
+
+            // Average over phi distribution
+            for(size_t o=0; o<weights_psi.size(); o++) {
+              dp.cyl_psi = weights_psi[o].value;
+              //Un-normalize by volume
+              double _ptvalue = weights_rad[i].weight
+                  * weights_len[j].weight
+                  * weights_rat[m].weight
+                  * weights_theta[k].weight
+                  * weights_phi[l].weight
+                  * weights_psi[o].weight
+                  * elliptical_cylinder_analytical_2DXY(&dp, qx, qy)
+              * pow(weights_rad[i].value,2)
+              * weights_len[j].value
+              * weights_rat[m].value;
+              if (weights_theta.size()>1) {
+                _ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
+              }
+              sum += _ptvalue;
+              //Find average volume
+              vol += weights_rad[i].weight
+                  * weights_len[j].weight
+                  * weights_rat[m].weight
+                  * pow(weights_rad[i].value,2)
+              * weights_len[j].value
+              * weights_rat[m].value;
+              //Find norm for volume
+              norm_vol += weights_rad[i].weight
+                  * weights_len[j].weight
+                  * weights_rat[m].weight;
+
+              norm += weights_rad[i].weight
+                  * weights_len[j].weight
+                  * weights_rat[m].weight
+                  * weights_theta[k].weight
+                  * weights_phi[l].weight
+                  * weights_psi[o].weight;
+
+            }
+          }
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_theta.size()>1) norm = norm / asin(1.0);
+
+  if (vol != 0.0 && norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm_vol);}
+
+  return sum/norm + background();
 
 }
@@ -256 +396 @@
  */
 double EllipticalCylinderModel :: evaluate_rphi(double q, double phi) {
-        double qx = q*cos(phi);
-        double qy = q*sin(phi);
-        return (*this).operator()(qx, qy);
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
 }
 /**
@@ -265 +405 @@
  */
 double EllipticalCylinderModel :: calculate_ER() {
-        EllipticalCylinderParameters dp;
-        dp.r_minor    = r_minor();
-        dp.r_ratio    = r_ratio();
-        dp.length     = length();
-        double rad_out = 0.0;
-        double suf_rad = sqrt(dp.r_minor*dp.r_minor*dp.r_ratio);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-
-        // Get the dispersion points for the r_minor
-        vector<WeightPoint> weights_rad;
-        r_minor.get_weights(weights_rad);
-
-        // Get the dispersion points for the r_ratio
-        vector<WeightPoint> weights_rat;
-        r_ratio.get_weights(weights_rat);
-
-        // Get the dispersion points for the length
-        vector<WeightPoint> weights_len;
-        length.get_weights(weights_len);
-
-        // Loop over minor radius weight points
-        for(size_t i=0; i<weights_rad.size(); i++) {
-                dp.r_minor = weights_rad[i].value;
-
-                // Loop over length weight points
-                for(size_t j=0; j<weights_len.size(); j++) {
-                        dp.length = weights_len[j].value;
-
-                        // Loop over r_ration weight points
-                        for(size_t m=0; m<weights_rat.size(); m++) {
-                                dp.r_ratio = weights_rat[m].value;
-                                //Calculate surface averaged radius
-                                suf_rad = sqrt(dp.r_minor * dp.r_minor * dp.r_ratio);
-
-                                //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-                                sum +=weights_rad[i].weight *weights_len[j].weight
-                                        * weights_rat[m].weight*DiamCyl(dp.length, suf_rad)/2.0;
-                                norm += weights_rad[i].weight *weights_len[j].weight* weights_rat[m].weight;
-                        }
-                }
-        }
-        if (norm != 0){
-                //return the averaged value
-                rad_out =  sum/norm;}
-        else{
-                //return normal value
-                //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-                rad_out = DiamCyl(dp.length, suf_rad)/2.0;}
-
-        return rad_out;
-}
+  EllipticalCylinderParameters dp;
+  dp.r_minor    = r_minor();
+  dp.r_ratio    = r_ratio();
+  dp.length     = length();
+  double rad_out = 0.0;
+  double suf_rad = sqrt(dp.r_minor*dp.r_minor*dp.r_ratio);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the r_minor
+  vector<WeightPoint> weights_rad;
+  r_minor.get_weights(weights_rad);
+
+  // Get the dispersion points for the r_ratio
+  vector<WeightPoint> weights_rat;
+  r_ratio.get_weights(weights_rat);
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_len;
+  length.get_weights(weights_len);
+
+  // Loop over minor radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp.r_minor = weights_rad[i].value;
+
+    // Loop over length weight points
+    for(size_t j=0; j<weights_len.size(); j++) {
+      dp.length = weights_len[j].value;
+
+      // Loop over r_ration weight points
+      for(size_t m=0; m<weights_rat.size(); m++) {
+        dp.r_ratio = weights_rat[m].value;
+        //Calculate surface averaged radius
+        suf_rad = sqrt(dp.r_minor * dp.r_minor * dp.r_ratio);
+
+        //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+        sum +=weights_rad[i].weight *weights_len[j].weight
+            * weights_rat[m].weight*DiamCyl(dp.length, suf_rad)/2.0;
+        norm += weights_rad[i].weight *weights_len[j].weight* weights_rat[m].weight;
+      }
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+    rad_out = DiamCyl(dp.length, suf_rad)/2.0;}
+
+  return rad_out;
+}

sansmodels/src/c_models/flexcyl_ellipX.cpp

@@ -21 +21 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "flexcyl_ellipX.h"
 
 extern "C" {
-        #include "libCylinder.h"
-        #include "libStructureFactor.h"
-        #include "flexcyl_ellipX.h"
-}
+#include "libCylinder.h"
+#include "libStructureFactor.h"
+}
+
+typedef struct {
+  double scale;
+  double length;
+  double kuhn_length;
+  double radius;
+  double axis_ratio;
+  double sldCyl;
+  double sldSolv;
+  double background;
+} FlexCylEXParameters;
 
 FlexCylEllipXModel :: FlexCylEllipXModel() {
-        scale      = Parameter(1.0);
-        length     = Parameter(1000.0, true);
-        length.set_min(0.0);
-        kuhn_length = Parameter(100.0, true);
-        kuhn_length.set_min(0.0);
-        radius  = Parameter(20.0, true);
-        radius.set_min(0.0);
-        axis_ratio = Parameter(1.5);
-        axis_ratio.set_min(0.0);
-        sldCyl   = Parameter(1.0e-6);
-        sldSolv   = Parameter(6.3e-6);
-        background = Parameter(0.0001);
+  scale      = Parameter(1.0);
+  length     = Parameter(1000.0, true);
+  length.set_min(0.0);
+  kuhn_length = Parameter(100.0, true);
+  kuhn_length.set_min(0.0);
+  radius  = Parameter(20.0, true);
+  radius.set_min(0.0);
+  axis_ratio = Parameter(1.5);
+  axis_ratio.set_min(0.0);
+  sldCyl   = Parameter(1.0e-6);
+  sldSolv   = Parameter(6.3e-6);
+  background = Parameter(0.0001);
 }
 
@@ -54 +64 @@
  */
 double FlexCylEllipXModel :: operator()(double q) {
-        double dp[8];
-
-        // Fill parameter array for IGOR library
-        // Add the background after averaging
-        dp[0] = scale();
-        dp[1] = length();
-        dp[2] = kuhn_length();
-        dp[3] = radius();
-        dp[4] = axis_ratio();
-        dp[5] = sldCyl();
-        dp[6] = sldSolv();
-        dp[7] = 0.0;
-
-        // Get the dispersion points for the length
-        vector<WeightPoint> weights_len;
-        length.get_weights(weights_len);
-
-        // Get the dispersion points for the kuhn_length
-        vector<WeightPoint> weights_kuhn;
-        kuhn_length.get_weights(weights_kuhn);
-
-        // Get the dispersion points for the radius
-        vector<WeightPoint> weights_rad;
-        radius.get_weights(weights_rad);
-
-        // Get the dispersion points for the axis_ratio
-        vector<WeightPoint> weights_ratio;
-        axis_ratio.get_weights(weights_ratio);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double vol = 0.0;
-
-        // Loop over length weight points
-        for(int i=0; i< (int)weights_len.size(); i++) {
-                dp[1] = weights_len[i].value;
-
-                // Loop over kuhn_length weight points
-                for(int j=0; j< (int)weights_kuhn.size(); j++) {
-                        dp[2] = weights_kuhn[j].value;
-
-                        // Loop over radius weight points
-                        for(int k=0; k< (int)weights_rad.size(); k++) {
-                                dp[3] = weights_rad[k].value;
-                                // Loop over axis_ratio weight points
-                                for(int l=0; l< (int)weights_ratio.size(); l++) {
-                                        dp[4] = weights_ratio[l].value;
-
-                                        //Un-normalize by volume
-                                        sum += weights_len[i].weight * weights_kuhn[j].weight*weights_rad[k].weight
-                                                * weights_ratio[l].weight * FlexCyl_Ellip(dp, q)
-                                                * (pow(weights_rad[k].value,2.0) * weights_ratio[l].value * weights_len[i].value);
-                                        //Find weighted volume
-                                        vol += weights_rad[k].weight * weights_kuhn[j].weight
-                                                * weights_len[i].weight * weights_ratio[l].weight
-                                                *pow(weights_rad[k].value,2.0)* weights_ratio[l].weight*weights_len[i].value;
-                                        norm += weights_len[i].weight * weights_kuhn[j].weight
-                                                *weights_rad[k].weight* weights_ratio[l].weight;
-                                }
-                        }
-                }
-        }
-        if (vol != 0.0 && norm != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm);}
-
-        return sum/norm + background();
+  double dp[8];
+
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = length();
+  dp[2] = kuhn_length();
+  dp[3] = radius();
+  dp[4] = axis_ratio();
+  dp[5] = sldCyl();
+  dp[6] = sldSolv();
+  dp[7] = 0.0;
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_len;
+  length.get_weights(weights_len);
+
+  // Get the dispersion points for the kuhn_length
+  vector<WeightPoint> weights_kuhn;
+  kuhn_length.get_weights(weights_kuhn);
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad;
+  radius.get_weights(weights_rad);
+
+  // Get the dispersion points for the axis_ratio
+  vector<WeightPoint> weights_ratio;
+  axis_ratio.get_weights(weights_ratio);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over length weight points
+  for(int i=0; i< (int)weights_len.size(); i++) {
+    dp[1] = weights_len[i].value;
+
+    // Loop over kuhn_length weight points
+    for(int j=0; j< (int)weights_kuhn.size(); j++) {
+      dp[2] = weights_kuhn[j].value;
+
+      // Loop over radius weight points
+      for(int k=0; k< (int)weights_rad.size(); k++) {
+        dp[3] = weights_rad[k].value;
+        // Loop over axis_ratio weight points
+        for(int l=0; l< (int)weights_ratio.size(); l++) {
+          dp[4] = weights_ratio[l].value;
+
+          //Un-normalize by volume
+          sum += weights_len[i].weight * weights_kuhn[j].weight*weights_rad[k].weight
+              * weights_ratio[l].weight * FlexCyl_Ellip(dp, q)
+          * (pow(weights_rad[k].value,2.0) * weights_ratio[l].value * weights_len[i].value);
+          //Find weighted volume
+          vol += weights_rad[k].weight * weights_kuhn[j].weight
+              * weights_len[i].weight * weights_ratio[l].weight
+              *pow(weights_rad[k].value,2.0)* weights_ratio[l].weight*weights_len[i].value;
+          norm += weights_len[i].weight * weights_kuhn[j].weight
+              *weights_rad[k].weight* weights_ratio[l].weight;
+        }
+      }
+    }
+  }
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -131 +141 @@
  */
 double FlexCylEllipXModel :: operator()(double qx, double qy) {
-        double q = sqrt(qx*qx + qy*qy);
-        return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -143 +153 @@
  */
 double FlexCylEllipXModel :: evaluate_rphi(double q, double phi) {
-        //double qx = q*cos(phi);
-        //double qy = q*sin(phi);
-        return (*this).operator()(q);
+  //double qx = q*cos(phi);
+  //double qy = q*sin(phi);
+  return (*this).operator()(q);
 }
 /**
@@ -152 +162 @@
  */
 double FlexCylEllipXModel :: calculate_ER() {
-        FlexCylEXParameters dp;
-
-        dp.radius  = radius();
-        dp.length     = length();
-        dp.axis_ratio  = axis_ratio();
-
-        double rad_out = 0.0;
-        double suf_rad = sqrt(dp.radius*dp.radius*dp.axis_ratio );
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-
-        // Get the dispersion points for the total length
-        vector<WeightPoint> weights_length;
-        length.get_weights(weights_length);
-
-        // Get the dispersion points for minor radius
-        vector<WeightPoint> weights_radius ;
-        radius.get_weights(weights_radius);
-
-        // Get the dispersion points for axis ratio = major_radius/minor_radius
-        vector<WeightPoint> weights_ratio ;
-        axis_ratio.get_weights(weights_ratio);
-
-        // Loop over major shell weight points
-        for(int i=0; i< (int)weights_length.size(); i++) {
-                dp.length = weights_length[i].value;
-                for(int k=0; k< (int)weights_radius.size(); k++) {
-                        dp.radius = weights_radius[k].value;
-                        // Loop over axis_ratio weight points
-                        for(int l=0; l< (int)weights_ratio.size(); l++) {
-                                dp.axis_ratio = weights_ratio[l].value;
-                                suf_rad = sqrt(dp.radius  * dp.radius  * dp.axis_ratio);
-                                //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-                                sum +=weights_length[i].weight * weights_radius[k].weight
-                                                * weights_ratio[l].weight *DiamCyl(dp.length,suf_rad)/2.0;
-                                norm += weights_length[i].weight* weights_radius[k].weight* weights_ratio[l].weight;
-                        }
-                }
-        }
-        if (norm != 0){
-                //return the averaged value
-                rad_out =  sum/norm;}
-        else{
-                //return normal value
-                //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-                rad_out = DiamCyl(dp.length,suf_rad)/2.0;}
-
-        return rad_out;
-}
+  FlexCylEXParameters dp;
+
+  dp.radius  = radius();
+  dp.length     = length();
+  dp.axis_ratio  = axis_ratio();
+
+  double rad_out = 0.0;
+  double suf_rad = sqrt(dp.radius*dp.radius*dp.axis_ratio );
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the total length
+  vector<WeightPoint> weights_length;
+  length.get_weights(weights_length);
+
+  // Get the dispersion points for minor radius
+  vector<WeightPoint> weights_radius ;
+  radius.get_weights(weights_radius);
+
+  // Get the dispersion points for axis ratio = major_radius/minor_radius
+  vector<WeightPoint> weights_ratio ;
+  axis_ratio.get_weights(weights_ratio);
+
+  // Loop over major shell weight points
+  for(int i=0; i< (int)weights_length.size(); i++) {
+    dp.length = weights_length[i].value;
+    for(int k=0; k< (int)weights_radius.size(); k++) {
+      dp.radius = weights_radius[k].value;
+      // Loop over axis_ratio weight points
+      for(int l=0; l< (int)weights_ratio.size(); l++) {
+        dp.axis_ratio = weights_ratio[l].value;
+        suf_rad = sqrt(dp.radius  * dp.radius  * dp.axis_ratio);
+        //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+        sum +=weights_length[i].weight * weights_radius[k].weight
+            * weights_ratio[l].weight *DiamCyl(dp.length,suf_rad)/2.0;
+        norm += weights_length[i].weight* weights_radius[k].weight* weights_ratio[l].weight;
+      }
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+    rad_out = DiamCyl(dp.length,suf_rad)/2.0;}
+
+  return rad_out;
+}

sansmodels/src/c_models/flexiblecylinder.cpp

@@ -18 +18 @@
  *   sansmodels/src/libigor
  *
- *      TODO: refactor so that we pull in the old sansmodels.c_extensions
  *      TODO: add 2d
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "flexible_cylinder.h"
 
 extern "C" {
         #include "libCylinder.h"
         #include "libStructureFactor.h"
-        #include "flexible_cylinder.h"
 }
+
+typedef struct {
+  double scale;
+  double length;
+  double kuhn_length;
+  double radius;
+  double sldCyl;
+  double sldSolv;
+  double background;
+} FlexibleCylinderParameters;
 
 FlexibleCylinderModel :: FlexibleCylinderModel() {

sansmodels/src/c_models/fractal.cpp

@@ -21 +21 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "fractal.h"
 
 extern "C" {
-        #include "libTwoPhase.h"
-        #include "fractal.h"
+#include "libTwoPhase.h"
 }
 
 FractalModel :: FractalModel() {
-        scale      = Parameter(1.0);
-    radius = Parameter(5.0, true);
-    radius.set_min(0.0);
-    fractal_dim = Parameter(2.0);
-    cor_length = Parameter(100.0);
-    sldBlock = Parameter(2.0e-6);
-    sldSolv= Parameter(6.35e-6);
-        background = Parameter(0.0);
+  scale      = Parameter(1.0);
+  radius = Parameter(5.0, true);
+  radius.set_min(0.0);
+  fractal_dim = Parameter(2.0);
+  cor_length = Parameter(100.0);
+  sldBlock = Parameter(2.0e-6);
+  sldSolv= Parameter(6.35e-6);
+  background = Parameter(0.0);
 }
 
@@ -49 +48 @@
  */
 double FractalModel :: operator()(double q) {
-        double dp[7];
+  double dp[7];
 
-        // Fill parameter array for IGOR library
-        // Add the background after averaging
-        dp[0] = scale();
-        dp[1] = radius();
-        dp[2] = fractal_dim();
-        dp[3] = cor_length();
-        dp[4] = sldBlock();
-        dp[5] = sldSolv();
-        dp[6] = 0.0;
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = radius();
+  dp[2] = fractal_dim();
+  dp[3] = cor_length();
+  dp[4] = sldBlock();
+  dp[5] = sldSolv();
+  dp[6] = 0.0;
 
-        // Get the dispersion points for the radius
-        vector<WeightPoint> weights_rad;
-        radius.get_weights(weights_rad);
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_rad;
+  radius.get_weights(weights_rad);
 
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double vol = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
 
-        // Loop over radius weight points
-        for(size_t i=0; i<weights_rad.size(); i++) {
-                dp[1] = weights_rad[i].value;
+  // Loop over radius weight points
+  for(size_t i=0; i<weights_rad.size(); i++) {
+    dp[1] = weights_rad[i].value;
 
-                //Un-normalize Fractal by volume
-                sum += weights_rad[i].weight
-                        * Fractal(dp, fabs(q)) * pow(weights_rad[i].value,3);
-                //Find average volume
-                vol += weights_rad[i].weight
-                        * pow(weights_rad[i].value,3);
+    //Un-normalize Fractal by volume
+    sum += weights_rad[i].weight
+        * Fractal(dp, fabs(q)) * pow(weights_rad[i].value,3);
+    //Find average volume
+    vol += weights_rad[i].weight
+        * pow(weights_rad[i].value,3);
 
-                norm += weights_rad[i].weight;
-        }
+    norm += weights_rad[i].weight;
+  }
 
-        if (vol != 0.0 && norm != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm);}
-        return sum/norm + background();
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+  return sum/norm + background();
 }
 
@@ -97 +96 @@
  */
 double FractalModel :: operator()(double qx, double qy) {
-        double q = sqrt(qx*qx + qy*qy);
-        return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -109 +108 @@
  */
 double FractalModel :: evaluate_rphi(double q, double phi) {
-        return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 
@@ -117 +116 @@
  */
 double FractalModel :: calculate_ER() {
-        //NOT implemented yet!!! 'cause None shape Model
-        return 0.0;
+  //NOT implemented yet!!! 'cause None shape Model
+  return 0.0;
 }

sansmodels/src/c_models/hollowcylinder.cpp

@@ -18 +18 @@
  *   sansmodels/src/libigor
  *
- *      TODO: refactor so that we pull in the old sansmodels.c_extensions
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "hollow_cylinder.h"
 
 extern "C" {
-        #include "libCylinder.h"
-        #include "libStructureFactor.h"
-        #include "hollow_cylinder.h"
+#include "libCylinder.h"
+#include "libStructureFactor.h"
+}
+
+typedef struct {
+  double scale;
+  double core_radius;
+  double radius;
+  double length;
+  double sldCyl;
+  double sldSolv;
+  double background;
+  double axis_theta;
+  double axis_phi;
+} HollowCylinderParameters;
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the hollow cylinder
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double hollow_cylinder_analytical_2D_scaled(HollowCylinderParameters *pars, double q, double q_x, double q_y) {
+  double cyl_x, cyl_y, cyl_z;
+  double q_z;
+  double  alpha,vol, cos_val;
+  double answer;
+  //convert angle degree to radian
+  double pi = 4.0*atan(1.0);
+  double theta = pars->axis_theta * pi/180.0;
+  double phi = pars->axis_phi * pi/180.0;
+
+  // Cylinder orientation
+  cyl_x = sin(theta) * cos(phi);
+  cyl_y = sin(theta) * sin(phi);
+  cyl_z = cos(theta);
+
+  // q vector
+  q_z = 0;
+
+  // Compute the angle btw vector q and the
+  // axis of the cylinder
+  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
+
+  // The following test should always pass
+  if (fabs(cos_val)>1.0) {
+    printf("core_shell_cylinder_analytical_2D: Unexpected error: cos(alpha)=%g\n", cos_val);
+    return 0;
+  }
+
+  alpha = acos( cos_val );
+
+  // Call the IGOR library function to get the kernel
+  answer = HolCylKernel(q, pars->core_radius, pars->radius, pars->length, cos_val);
+
+  // Multiply by contrast^2
+  answer *= (pars->sldCyl - pars->sldSolv)*(pars->sldCyl - pars->sldSolv);
+
+  //normalize by cylinder volume
+  vol=pi*((pars->radius*pars->radius)-(pars->core_radius *pars->core_radius))
+          *(pars->length);
+  answer *= vol;
+
+  //convert to [cm-1]
+  answer *= 1.0e8;
+
+  //Scale
+  answer *= pars->scale;
+
+  // add in the background
+  answer += pars->background;
+
+  return answer;
+}
+
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the Hollow cylinder
+ * @param q: q-value
+ * @return: function value
+ */
+static double hollow_cylinder_analytical_2DXY(HollowCylinderParameters *pars, double qx, double qy) {
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return hollow_cylinder_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
 HollowCylinderModel :: HollowCylinderModel() {
-        scale      = Parameter(1.0);
-        core_radius = Parameter(20.0, true);
-        core_radius.set_min(0.0);
-        radius  = Parameter(30.0, true);
-        radius.set_min(0.0);
-        length     = Parameter(400.0, true);
-        length.set_min(0.0);
-        sldCyl  = Parameter(6.3e-6);
-        sldSolv  = Parameter(1.0e-6);
-        background = Parameter(0.0);
-        axis_theta = Parameter(0.0, true);
-        axis_phi   = Parameter(0.0, true);
+  scale      = Parameter(1.0);
+  core_radius = Parameter(20.0, true);
+  core_radius.set_min(0.0);
+  radius  = Parameter(30.0, true);
+  radius.set_min(0.0);
+  length     = Parameter(400.0, true);
+  length.set_min(0.0);
+  sldCyl  = Parameter(6.3e-6);
+  sldSolv  = Parameter(1.0e-6);
+  background = Parameter(0.0);
+  axis_theta = Parameter(0.0, true);
+  axis_phi   = Parameter(0.0, true);
 }
 
@@ -55 +140 @@
  */
 double HollowCylinderModel :: operator()(double q) {
-        double dp[7];
-
-        dp[0] = scale();
-        dp[1] = core_radius();
-        dp[2] = radius();
-        dp[3] = length();
-        dp[4] = sldCyl();
-        dp[5] = sldSolv();
-        dp[6] = 0.0;
-
-        // Get the dispersion points for the core radius
-        vector<WeightPoint> weights_core_radius;
-        core_radius.get_weights(weights_core_radius);
-
-        // Get the dispersion points for the shell radius
-        vector<WeightPoint> weights_radius;
-        radius.get_weights(weights_radius);
-
-        // Get the dispersion points for the length
-        vector<WeightPoint> weights_length;
-        length.get_weights(weights_length);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double vol = 0.0;
-
-        // Loop over core radius weight points
-        for(int i=0; i< (int)weights_core_radius.size(); i++) {
-                dp[1] = weights_core_radius[i].value;
-
-                // Loop over length weight points
-                for(int j=0; j< (int)weights_length.size(); j++) {
-                        dp[3] = weights_length[j].value;
-
-                        // Loop over shell radius weight points
-                        for(int k=0; k< (int)weights_radius.size(); k++) {
-                                dp[2] = weights_radius[k].value;
-                                //Un-normalize  by volume
-                                sum += weights_core_radius[i].weight
-                                        * weights_length[j].weight
-                                        * weights_radius[k].weight
-                                        * HollowCylinder(dp, q)
-                                        * (pow(weights_radius[k].value,2)-pow(weights_core_radius[i].value,2))
-                                        * weights_length[j].value;
-                                //Find average volume
-                                vol += weights_core_radius[i].weight
-                                        * weights_length[j].weight
-                                        * weights_radius[k].weight
-                                        * (pow(weights_radius[k].value,2)-pow(weights_core_radius[i].value,2))
-                                        * weights_length[j].value;
-
-                                norm += weights_core_radius[i].weight
-                                        * weights_length[j].weight
-                                        * weights_radius[k].weight;
-                        }
-                }
-        }
-        if (vol != 0.0 && norm != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm);}
-
-        return sum/norm + background();
+  double dp[7];
+
+  dp[0] = scale();
+  dp[1] = core_radius();
+  dp[2] = radius();
+  dp[3] = length();
+  dp[4] = sldCyl();
+  dp[5] = sldSolv();
+  dp[6] = 0.0;
+
+  // Get the dispersion points for the core radius
+  vector<WeightPoint> weights_core_radius;
+  core_radius.get_weights(weights_core_radius);
+
+  // Get the dispersion points for the shell radius
+  vector<WeightPoint> weights_radius;
+  radius.get_weights(weights_radius);
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_length;
+  length.get_weights(weights_length);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over core radius weight points
+  for(int i=0; i< (int)weights_core_radius.size(); i++) {
+    dp[1] = weights_core_radius[i].value;
+
+    // Loop over length weight points
+    for(int j=0; j< (int)weights_length.size(); j++) {
+      dp[3] = weights_length[j].value;
+
+      // Loop over shell radius weight points
+      for(int k=0; k< (int)weights_radius.size(); k++) {
+        dp[2] = weights_radius[k].value;
+        //Un-normalize  by volume
+        sum += weights_core_radius[i].weight
+            * weights_length[j].weight
+            * weights_radius[k].weight
+            * HollowCylinder(dp, q)
+        * (pow(weights_radius[k].value,2)-pow(weights_core_radius[i].value,2))
+        * weights_length[j].value;
+        //Find average volume
+        vol += weights_core_radius[i].weight
+            * weights_length[j].weight
+            * weights_radius[k].weight
+            * (pow(weights_radius[k].value,2)-pow(weights_core_radius[i].value,2))
+            * weights_length[j].value;
+
+        norm += weights_core_radius[i].weight
+            * weights_length[j].weight
+            * weights_radius[k].weight;
+      }
+    }
+  }
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -127 +212 @@
  */
 double HollowCylinderModel :: operator()(double qx, double qy) {
-        HollowCylinderParameters dp;
-        // Fill parameter array
-        dp.scale      = scale();
-        dp.core_radius     = core_radius();
-        dp.radius  = radius();
-        dp.length     = length();
-        dp.sldCyl   = sldCyl();
-        dp.sldSolv  = sldSolv();
-        dp.background = 0.0;
-        dp.axis_theta = axis_theta();
-        dp.axis_phi   = axis_phi();
-
-        // Get the dispersion points for the core radius
-        vector<WeightPoint> weights_core_radius;
-        core_radius.get_weights(weights_core_radius);
-
-        // Get the dispersion points for the shell radius
-        vector<WeightPoint> weights_radius;
-        radius.get_weights(weights_radius);
-
-        // Get the dispersion points for the length
-        vector<WeightPoint> weights_length;
-        length.get_weights(weights_length);
-
-        // Get angular averaging for theta
-        vector<WeightPoint> weights_theta;
-        axis_theta.get_weights(weights_theta);
-
-        // Get angular averaging for phi
-        vector<WeightPoint> weights_phi;
-        axis_phi.get_weights(weights_phi);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double norm_vol = 0.0;
-        double vol = 0.0;
-        double pi = 4.0*atan(1.0);
-        // Loop over core radius weight points
-        for(int i=0; i<(int)weights_core_radius.size(); i++) {
-                dp.core_radius = weights_core_radius[i].value;
-
-
-                // Loop over length weight points
-                for(int j=0; j<(int)weights_length.size(); j++) {
-                        dp.length = weights_length[j].value;
-
-                        // Loop over shell radius weight points
-                        for(int m=0; m< (int)weights_radius.size(); m++) {
-                                dp.radius = weights_radius[m].value;
-
-                        // Average over theta distribution
-                        for(int k=0; k< (int)weights_theta.size(); k++) {
-                                dp.axis_theta = weights_theta[k].value;
-
-                                // Average over phi distribution
-                                for(int l=0; l< (int)weights_phi.size(); l++) {
-                                        dp.axis_phi = weights_phi[l].value;
-                                        //Un-normalize by volume
-                                        double _ptvalue = weights_core_radius[i].weight
-                                                * weights_length[j].weight
-                                                * weights_radius[m].weight
-                                                * weights_theta[k].weight
-                                                * weights_phi[l].weight
-                                                * hollow_cylinder_analytical_2DXY(&dp, qx, qy)
-                                                / ((pow(weights_radius[m].value,2)-pow(weights_core_radius[i].value,2))
-                                                * weights_length[j].value);
-                                        if (weights_theta.size()>1) {
-                                                _ptvalue *= fabs(sin(weights_theta[k].value * pi/180.0));
-                                        }
-                                        sum += _ptvalue;
-                                        //Find average volume
-                                        vol += weights_core_radius[i].weight
-                                                * weights_length[j].weight
-                                                * weights_radius[m].weight
-                                                * (pow(weights_radius[m].value,2)-pow(weights_core_radius[i].value,2))
-                                                * weights_length[j].value;
-                                        //Find norm for volume
-                                        norm_vol += weights_core_radius[i].weight
-                                                * weights_length[j].weight
-                                                * weights_radius[m].weight;
-
-                                        norm += weights_core_radius[i].weight
-                                                * weights_length[j].weight
-                                                * weights_radius[m].weight
-                                                * weights_theta[k].weight
-                                                * weights_phi[l].weight;
-
-                                }
-                        }
-                        }
-                }
-        }
-        // Averaging in theta needs an extra normalization
-        // factor to account for the sin(theta) term in the
-        // integration (see documentation).
-        if (weights_theta.size()>1) norm = norm/asin(1.0);
-        if (vol != 0.0 || norm_vol != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum*(vol/norm_vol);}
-        return sum/norm + background();
+  HollowCylinderParameters dp;
+  // Fill parameter array
+  dp.scale      = scale();
+  dp.core_radius     = core_radius();
+  dp.radius  = radius();
+  dp.length     = length();
+  dp.sldCyl   = sldCyl();
+  dp.sldSolv  = sldSolv();
+  dp.background = 0.0;
+  dp.axis_theta = axis_theta();
+  dp.axis_phi   = axis_phi();
+
+  // Get the dispersion points for the core radius
+  vector<WeightPoint> weights_core_radius;
+  core_radius.get_weights(weights_core_radius);
+
+  // Get the dispersion points for the shell radius
+  vector<WeightPoint> weights_radius;
+  radius.get_weights(weights_radius);
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_length;
+  length.get_weights(weights_length);
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_theta;
+  axis_theta.get_weights(weights_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_phi;
+  axis_phi.get_weights(weights_phi);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double norm_vol = 0.0;
+  double vol = 0.0;
+  double pi = 4.0*atan(1.0);
+  // Loop over core radius weight points
+  for(int i=0; i<(int)weights_core_radius.size(); i++) {
+    dp.core_radius = weights_core_radius[i].value;
+
+
+    // Loop over length weight points
+    for(int j=0; j<(int)weights_length.size(); j++) {
+      dp.length = weights_length[j].value;
+
+      // Loop over shell radius weight points
+      for(int m=0; m< (int)weights_radius.size(); m++) {
+        dp.radius = weights_radius[m].value;
+
+        // Average over theta distribution
+        for(int k=0; k< (int)weights_theta.size(); k++) {
+          dp.axis_theta = weights_theta[k].value;
+
+          // Average over phi distribution
+          for(int l=0; l< (int)weights_phi.size(); l++) {
+            dp.axis_phi = weights_phi[l].value;
+            //Un-normalize by volume
+            double _ptvalue = weights_core_radius[i].weight
+                * weights_length[j].weight
+                * weights_radius[m].weight
+                * weights_theta[k].weight
+                * weights_phi[l].weight
+                * hollow_cylinder_analytical_2DXY(&dp, qx, qy)
+            / ((pow(weights_radius[m].value,2)-pow(weights_core_radius[i].value,2))
+                * weights_length[j].value);
+            if (weights_theta.size()>1) {
+              _ptvalue *= fabs(sin(weights_theta[k].value * pi/180.0));
+            }
+            sum += _ptvalue;
+            //Find average volume
+            vol += weights_core_radius[i].weight
+                * weights_length[j].weight
+                * weights_radius[m].weight
+                * (pow(weights_radius[m].value,2)-pow(weights_core_radius[i].value,2))
+                * weights_length[j].value;
+            //Find norm for volume
+            norm_vol += weights_core_radius[i].weight
+                * weights_length[j].weight
+                * weights_radius[m].weight;
+
+            norm += weights_core_radius[i].weight
+                * weights_length[j].weight
+                * weights_radius[m].weight
+                * weights_theta[k].weight
+                * weights_phi[l].weight;
+
+          }
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_theta.size()>1) norm = norm/asin(1.0);
+  if (vol != 0.0 || norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum*(vol/norm_vol);}
+  return sum/norm + background();
 }
 
@@ -238 +323 @@
  */
 double HollowCylinderModel :: evaluate_rphi(double q, double phi) {
-        double qx = q*cos(phi);
-        double qy = q*sin(phi);
-        return (*this).operator()(qx, qy);
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
 }
 /**
@@ -247 +332 @@
  */
 double HollowCylinderModel :: calculate_ER() {
-        HollowCylinderParameters dp;
-
-        dp.radius  = radius();
-        dp.length     = length();
-
-        double rad_out = 0.0;
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-
-        // Get the dispersion points for the major shell
-        vector<WeightPoint> weights_length;
-        length.get_weights(weights_length);
-
-        // Get the dispersion points for the minor shell
-        vector<WeightPoint> weights_radius ;
-        radius.get_weights(weights_radius);
-
-        // Loop over major shell weight points
-        for(int i=0; i< (int)weights_length.size(); i++) {
-                dp.length = weights_length[i].value;
-                for(int k=0; k< (int)weights_radius.size(); k++) {
-                        dp.radius = weights_radius[k].value;
-                        //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-                        sum +=weights_length[i].weight
-                                * weights_radius[k].weight*DiamCyl(dp.length,dp.radius)/2.0;
-                        norm += weights_length[i].weight* weights_radius[k].weight;
-                }
-        }
-        if (norm != 0){
-                //return the averaged value
-                rad_out =  sum/norm;}
-        else{
-                //return normal value
-                //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-                rad_out = DiamCyl(dp.length,dp.radius)/2.0;}
-
-        return rad_out;
-}
+  HollowCylinderParameters dp;
+
+  dp.radius  = radius();
+  dp.length     = length();
+
+  double rad_out = 0.0;
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the major shell
+  vector<WeightPoint> weights_length;
+  length.get_weights(weights_length);
+
+  // Get the dispersion points for the minor shell
+  vector<WeightPoint> weights_radius ;
+  radius.get_weights(weights_radius);
+
+  // Loop over major shell weight points
+  for(int i=0; i< (int)weights_length.size(); i++) {
+    dp.length = weights_length[i].value;
+    for(int k=0; k< (int)weights_radius.size(); k++) {
+      dp.radius = weights_radius[k].value;
+      //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+      sum +=weights_length[i].weight
+          * weights_radius[k].weight*DiamCyl(dp.length,dp.radius)/2.0;
+      norm += weights_length[i].weight* weights_radius[k].weight;
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+    rad_out = DiamCyl(dp.length,dp.radius)/2.0;}
+
+  return rad_out;
+}

sansmodels/src/c_models/lamellar.cpp

@@ -18 +18 @@
  *   sansmodels/src/libigor
  *
- *      TODO: refactor so that we pull in the old sansmodels.c_extensions
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "lamellar.h"
 
-extern "C" {
-//      #include "libCylinder.h"
-        #include "lamellar.h"
+/*  LamellarFFX  :  calculates the form factor of a lamellar structure - no S(q) effects included
+            -NO polydispersion included
+*/
+static double lamellar_kernel(double dp[], double q){
+  double scale,del,sld_bi,sld_sol,contr,bkg;    //local variables of coefficient wave
+  double inten, qval,Pq;
+  double Pi;
+
+
+  Pi = 4.0*atan(1.0);
+  scale = dp[0];
+  del = dp[1];
+  sld_bi = dp[2];
+  sld_sol = dp[3];
+  bkg = dp[4];
+  qval = q;
+  contr = sld_bi -sld_sol;
+
+  Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del));
+
+  inten = 2.0*Pi*scale*Pq/(qval*qval);    //this is now dimensionless...
+
+  inten /= del;     //normalize by the thickness (in A)
+
+  inten *= 1.0e8;   // 1/A to 1/cm
+
+  return(inten+bkg);
 }
 
@@ -39 +62 @@
         sld_sol    = Parameter(6.3e-6);
         background = Parameter(0.0);
-
 }
 

sansmodels/src/c_models/lamellarFF_HG.cpp

@@ -18 +18 @@
  *   sansmodels/src/libigor
  *
- *      TODO: refactor so that we pull in the old sansmodels.c_extensions
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "lamellarFF_HG.h"
 
 extern "C" {
-        #include "libCylinder.h"
-        #include "lamellarFF_HG.h"
+#include "libCylinder.h"
 }
 
 LamellarFFHGModel :: LamellarFFHGModel() {
-        scale      = Parameter(1.0);
-        t_length     = Parameter(15.0, true);
-        t_length.set_min(0.0);
-        h_thickness    = Parameter(10.0, true);
-        h_thickness.set_min(0.0);
-        sld_tail   = Parameter(4e-7);
-        sld_head  = Parameter(3e-6);
-        sld_solvent    = Parameter(6e-6);
-        background = Parameter(0.0);
+  scale      = Parameter(1.0);
+  t_length     = Parameter(15.0, true);
+  t_length.set_min(0.0);
+  h_thickness    = Parameter(10.0, true);
+  h_thickness.set_min(0.0);
+  sld_tail   = Parameter(4e-7);
+  sld_head  = Parameter(3e-6);
+  sld_solvent    = Parameter(6e-6);
+  background = Parameter(0.0);
 
 }
@@ -52 +50 @@
  */
 double LamellarFFHGModel :: operator()(double q) {
-        double dp[7];
+  double dp[7];
 
-        // Fill parameter array for IGOR library
-        // Add the background after averaging
-        dp[0] = scale();
-        dp[1] = t_length();
-        dp[2] = h_thickness();
-        dp[3] = sld_tail();
-        dp[4] = sld_head();
-        dp[5] = sld_solvent();
-        dp[6] = 0.0;
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = t_length();
+  dp[2] = h_thickness();
+  dp[3] = sld_tail();
+  dp[4] = sld_head();
+  dp[5] = sld_solvent();
+  dp[6] = 0.0;
 
-        // Get the dispersion points for the tail length
-        vector<WeightPoint> weights_t_length;
-        t_length.get_weights(weights_t_length);
+  // Get the dispersion points for the tail length
+  vector<WeightPoint> weights_t_length;
+  t_length.get_weights(weights_t_length);
 
-        // Get the dispersion points for the head thickness
-        vector<WeightPoint> weights_h_thickness;
-        h_thickness.get_weights(weights_h_thickness);
+  // Get the dispersion points for the head thickness
+  vector<WeightPoint> weights_h_thickness;
+  h_thickness.get_weights(weights_h_thickness);
 
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
 
-        // Loop over semi axis A weight points
-        for(int i=0; i< (int)weights_t_length.size(); i++) {
-                dp[1] = weights_t_length[i].value;
+  // Loop over semi axis A weight points
+  for(int i=0; i< (int)weights_t_length.size(); i++) {
+    dp[1] = weights_t_length[i].value;
 
-                for (int j=0; j< (int)weights_h_thickness.size(); j++){
-                        dp[2] = weights_h_thickness[j].value;
+    for (int j=0; j< (int)weights_h_thickness.size(); j++){
+      dp[2] = weights_h_thickness[j].value;
 
-                        sum += weights_t_length[i].weight* weights_h_thickness[j].weight*LamellarFF_HG(dp, q);
-                        norm += weights_t_length[i].weight* weights_h_thickness[j].weight;
-                }
+      sum += weights_t_length[i].weight* weights_h_thickness[j].weight*LamellarFF_HG(dp, q);
+      norm += weights_t_length[i].weight* weights_h_thickness[j].weight;
+    }
 
-        }
-        return sum/norm + background();
+  }
+  return sum/norm + background();
 }
 
@@ -99 +97 @@
 
 double LamellarFFHGModel :: operator()(double qx, double qy) {
-        double q = sqrt(qx*qx + qy*qy);
-        return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -111 +109 @@
  */
 double LamellarFFHGModel :: evaluate_rphi(double q, double phi) {
-        return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 /**
@@ -118 +116 @@
  */
 double LamellarFFHGModel :: calculate_ER() {
-//NOT implemented yet!!!
-        return 0.0;
+  //NOT implemented yet!!!
+  return 0.0;
 }

sansmodels/src/c_models/lamellarPC.cpp

@@ -20 +20 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "lamellarPC.h"
 
 extern "C" {
-        #include "libCylinder.h"
-        #include "lamellarPC.h"
+#include "libCylinder.h"
 }
 
 LamellarPCrystalModel :: LamellarPCrystalModel() {
-        scale      = Parameter(1.0);
-        thickness     = Parameter(33.0, true);
-        thickness.set_min(0.0);
-        Nlayers    = Parameter(20.0, true);
-        Nlayers.set_min(0.0);
-        spacing   = Parameter(250);
-        pd_spacing   = Parameter(0.0);
-        sld_layer  = Parameter(1.0e-6);
-        sld_solvent    = Parameter(6.34e-6);
-        background = Parameter(0.0);
+  scale      = Parameter(1.0);
+  thickness     = Parameter(33.0, true);
+  thickness.set_min(0.0);
+  Nlayers    = Parameter(20.0, true);
+  Nlayers.set_min(0.0);
+  spacing   = Parameter(250);
+  pd_spacing   = Parameter(0.0);
+  sld_layer  = Parameter(1.0e-6);
+  sld_solvent    = Parameter(6.34e-6);
+  background = Parameter(0.0);
 
 }
@@ -51 +50 @@
  */
 double LamellarPCrystalModel :: operator()(double q) {
-        double dp[8];
+  double dp[8];
 
-        // Fill parameter array for IGOR library
-        // Add the background after averaging
-        dp[0] = scale();
-        dp[1] = thickness();
-        dp[2] = Nlayers();
-        dp[3] = spacing();
-        dp[4] = pd_spacing();
-        dp[5] = sld_layer();
-        dp[6] = sld_solvent();
-        dp[7] = 0.0; // Do not apply background here.
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = thickness();
+  dp[2] = Nlayers();
+  dp[3] = spacing();
+  dp[4] = pd_spacing();
+  dp[5] = sld_layer();
+  dp[6] = sld_solvent();
+  dp[7] = 0.0; // Do not apply background here.
 
-        // Get the dispersion points for the head thickness
-        vector<WeightPoint> weights_thickness;
-        thickness.get_weights(weights_thickness);
+  // Get the dispersion points for the head thickness
+  vector<WeightPoint> weights_thickness;
+  thickness.get_weights(weights_thickness);
 
-        // Let's provide from the func which is more accurate especially for small q region.
-        // Get the dispersion points for the tail length
-        //vector<WeightPoint> weights_spacing;
-        //spacing.get_weights(weights_spacing);
+  // Let's provide from the func which is more accurate especially for small q region.
+  // Get the dispersion points for the tail length
+  //vector<WeightPoint> weights_spacing;
+  //spacing.get_weights(weights_spacing);
 
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
 
-        // Loop over thickness and spacing weight points
-        for(int i=0; i< (int)weights_thickness.size(); i++) {
-                dp[1] = weights_thickness[i].value;
-                //for (int j=0; j< (int)weights_spacing.size(); j++){
-                        //dp[3] = weights_spacing[j].value;
-                        sum += weights_thickness[i].weight*Lamellar_ParaCrystal(dp, q);
-                        norm += weights_thickness[i].weight;
-                //}
-        }
-        //apply norm and background
-        return sum/norm + background();
+  // Loop over thickness and spacing weight points
+  for(int i=0; i< (int)weights_thickness.size(); i++) {
+    dp[1] = weights_thickness[i].value;
+    //for (int j=0; j< (int)weights_spacing.size(); j++){
+    //dp[3] = weights_spacing[j].value;
+    sum += weights_thickness[i].weight*Lamellar_ParaCrystal(dp, q);
+    norm += weights_thickness[i].weight;
+    //}
+  }
+  //apply norm and background
+  return sum/norm + background();
 }
 
@@ -98 +97 @@
 
 double LamellarPCrystalModel :: operator()(double qx, double qy) {
-        double q = sqrt(qx*qx + qy*qy);
-        return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -110 +109 @@
  */
 double LamellarPCrystalModel :: evaluate_rphi(double q, double phi) {
-        return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 /**
@@ -117 +116 @@
  */
 double LamellarPCrystalModel :: calculate_ER() {
-//NOT implemented yet!!!
-        return 0.0;
+  //NOT implemented yet!!!
+  return 0.0;
 }

sansmodels/src/c_models/lamellarPS.cpp

@@ -23 +23 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
+#include <stdlib.h>
 using namespace std;
+#include "lamellarPS.h"
 
 extern "C" {
-        #include "libCylinder.h"
-        #include "lamellarPS.h"
+#include "libCylinder.h"
+}
+
+/*LamellarPS_kernel() was moved from libigor to get rid of polydipersity in del(thickness) that we provide from control panel.
+  LamellarPSX  :  calculates the form factor of a lamellar structure - with S(q) effects included
+-------
+------- resolution effects ARE NOT included, but only a CONSTANT default value, not the real q-dependent resolution!!
+
+ */
+static double LamellarPS_kernel(double dp[], double q)
+{
+  double scale,dd,del,sld_bi,sld_sol,contr,NN,Cp,bkg;   //local variables of coefficient wave
+  double inten, qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ;
+  double Pi,Euler,dQDefault,fii;
+  int ii,NNint;
+  Euler = 0.5772156649;   // Euler's constant
+  dQDefault = 0.0;    //[=] 1/A, q-resolution, default value
+  dQ = dQDefault;
+
+  Pi = 4.0*atan(1.0);
+  qval = q;
+
+  scale = dp[0];
+  dd = dp[1];
+  del = dp[2];
+  sld_bi = dp[3];
+  sld_sol = dp[4];
+  NN = trunc(dp[5]);    //be sure that NN is an integer
+  Cp = dp[6];
+  bkg = dp[7];
+
+  contr = sld_bi - sld_sol;
+
+  Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del));
+
+  NNint = (int)NN;    //cast to an integer for the loop
+  ii=0;
+  Sq = 0.0;
+  for(ii=1;ii<(NNint-1);ii+=1) {
+
+    fii = (double)ii;   //do I really need to do this?
+
+    temp = 0.0;
+    alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler);
+    t1 = 2.0*dQ*dQ*dd*dd*alpha;
+    t2 = 2.0*qval*qval*dd*dd*alpha;
+    t3 = dQ*dQ*dd*dd*ii*ii;
+
+    temp = 1.0-ii/NN;
+    temp *= cos(dd*qval*ii/(1.0+t1));
+    temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
+    temp /= sqrt(1.0+t1);
+
+    Sq += temp;
+  }
+
+  Sq *= 2.0;
+  Sq += 1.0;
+
+  inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval);
+
+  inten *= 1.0e8;   // 1/A to 1/cm
+
+  return(inten+bkg);
 }
 
 LamellarPSModel :: LamellarPSModel() {
-        scale      = Parameter(1.0);
-        spacing    = Parameter(400.0, true);
-        spacing.set_min(0.0);
-        delta     = Parameter(30.0);
-        delta.set_min(0.0);
-        sld_bi   = Parameter(6.3e-6);
-        sld_sol   = Parameter(1.0e-6);
-        n_plates     = Parameter(20.0);
-        caille = Parameter(0.1);
-        background = Parameter(0.0);
+  scale      = Parameter(1.0);
+  spacing    = Parameter(400.0, true);
+  spacing.set_min(0.0);
+  delta     = Parameter(30.0);
+  delta.set_min(0.0);
+  sld_bi   = Parameter(6.3e-6);
+  sld_sol   = Parameter(1.0e-6);
+  n_plates     = Parameter(20.0);
+  caille = Parameter(0.1);
+  background = Parameter(0.0);
 
 }
@@ -54 +117 @@
  */
 double LamellarPSModel :: operator()(double q) {
-        double dp[8];
+  double dp[8];
 
-        // Fill parameter array for IGOR library
-        // Add the background after averaging
-        dp[0] = scale();
-        dp[1] = spacing();
-        dp[2] = delta();
-        dp[3] = sld_bi();
-        dp[4] = sld_sol();
-        dp[5] = n_plates();
-        dp[6] = caille();
-        dp[7] = 0.0;
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = spacing();
+  dp[2] = delta();
+  dp[3] = sld_bi();
+  dp[4] = sld_sol();
+  dp[5] = n_plates();
+  dp[6] = caille();
+  dp[7] = 0.0;
 
 
-        // Get the dispersion points for spacing and delta (thickness)
-        vector<WeightPoint> weights_spacing;
-        spacing.get_weights(weights_spacing);
-        vector<WeightPoint> weights_delta;
-        delta.get_weights(weights_delta);
+  // Get the dispersion points for spacing and delta (thickness)
+  vector<WeightPoint> weights_spacing;
+  spacing.get_weights(weights_spacing);
+  vector<WeightPoint> weights_delta;
+  delta.get_weights(weights_delta);
 
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
 
-        // Loop over short_edgeA weight points
-        for(int i=0; i< (int)weights_spacing.size(); i++) {
-                dp[1] = weights_spacing[i].value;
-                for(int j=0; j< (int)weights_spacing.size(); j++) {
-                        dp[2] = weights_delta[i].value;
+  // Loop over short_edgeA weight points
+  for(int i=0; i< (int)weights_spacing.size(); i++) {
+    dp[1] = weights_spacing[i].value;
+    for(int j=0; j< (int)weights_spacing.size(); j++) {
+      dp[2] = weights_delta[i].value;
 
-                        sum += weights_spacing[i].weight * weights_delta[j].weight * LamellarPS_kernel(dp, q);
-                        norm += weights_spacing[i].weight * weights_delta[j].weight;
-                }
-        }
-        return sum/norm + background();
+      sum += weights_spacing[i].weight * weights_delta[j].weight * LamellarPS_kernel(dp, q);
+      norm += weights_spacing[i].weight * weights_delta[j].weight;
+    }
+  }
+  return sum/norm + background();
 }
 /**
@@ -97 +160 @@
  */
 double LamellarPSModel :: operator()(double qx, double qy) {
-        double q = sqrt(qx*qx + qy*qy);
-        return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -109 +172 @@
  */
 double LamellarPSModel :: evaluate_rphi(double q, double phi) {
-        return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 /**
@@ -116 +179 @@
  */
 double LamellarPSModel :: calculate_ER() {
-//NOT implemented yet!!!
-        return 0.0;
+  //NOT implemented yet!!!
+  return 0.0;
 }
 

sansmodels/src/c_models/lamellarPS_HG.cpp

@@ -18 +18 @@
  *   sansmodels/src/libigor
  *
- *      TODO: refactor so that we pull in the old sansmodels.c_extensions
  *      TODO: add 2D function
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "lamellarPS_HG.h"
 
 extern "C" {
         #include "libCylinder.h"
-        #include "lamellarPS_HG.h"
 }
 
@@ -137 +135 @@
         return 0.0;
 }
-
-/*
-double LamellarPSHGModel :: operator()(double qx, double qy) {
-        LamellarPSHGParameters dp;
-        // Fill parameter array
-        dp.scale      = scale();
-        dp.spacing   = spacing();
-        dp.deltaT  = deltaT();
-        dp.deltaH = deltaH();
-        dp.sld_tail   = sld_tail();
-        dp.sld_head = sld_head();
-        dp.sld_solvent   = sld_solvent();
-        dp.n_plates = n_plates();
-        dp.caille = caille();
-        dp.background    = background();
-
-        // Get the dispersion points for the deltaT
-        vector<WeightPoint> weights_deltaT;
-        deltaT.get_weights(weights_deltaT);
-
-        // Get the dispersion points for the deltaH
-        vector<WeightPoint> weights_deltaH;
-        deltaH.get_weights(weights_deltaH);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-
-        // Loop over deltaT weight points
-        for(int i=0; i< (int)weights_deltaT.size(); i++) {
-                dp.deltaT = weights_deltaT[i].value;
-
-                // Loop over deltaH weight points
-                for(int j=0; j< (int)weights_deltaH.size(); j++) {
-                        dp.deltaH = weights_deltaH[j].value;
-
-                        sum += weights_deltaT[i].weight *weights_deltaH[j].weight *lamellarPS_HG_analytical_2DXY(&dp, qx, qy);
-                        norm += weights_deltaT[i].weight * weights_deltaH[j].weight;
-                }
-        }
-        return sum/norm + background();
-}
-*/
-
-

sansmodels/src/c_models/multishell.cpp

@@ -21 +21 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "multishell.h"
 
 extern "C" {
-        #include "libSphere.h"
-        #include "multishell.h"
+#include "libSphere.h"
 }
 
+typedef struct {
+  double scale;
+  double core_radius;
+  double s_thickness;
+  double w_thickness;
+  double core_sld;
+  double shell_sld;
+  double n_pairs;
+  double background;
+
+} MultiShellParameters;
+
 MultiShellModel :: MultiShellModel() {
-        scale      = Parameter(1.0);
-        core_radius     = Parameter(60.0, true);
-        core_radius.set_min(0.0);
-        s_thickness  = Parameter(10.0, true);
-        s_thickness.set_min(0.0);
-        w_thickness   = Parameter(10.0, true);
-        w_thickness.set_min(0.0);
-        core_sld   = Parameter(6.4e-6);
-        shell_sld   = Parameter(4.0e-7);
-        n_pairs   = Parameter(2);
-        background = Parameter(0.0);
+  scale      = Parameter(1.0);
+  core_radius     = Parameter(60.0, true);
+  core_radius.set_min(0.0);
+  s_thickness  = Parameter(10.0, true);
+  s_thickness.set_min(0.0);
+  w_thickness   = Parameter(10.0, true);
+  w_thickness.set_min(0.0);
+  core_sld   = Parameter(6.4e-6);
+  shell_sld   = Parameter(4.0e-7);
+  n_pairs   = Parameter(2);
+  background = Parameter(0.0);
 }
 
@@ -52 +63 @@
  */
 double MultiShellModel :: operator()(double q) {
-        double dp[8];
+  double dp[8];
 
-        // Fill parameter array for IGOR library
-        // Add the background after averaging
-        dp[0] = scale();
-        dp[1] = core_radius();
-        dp[2] = s_thickness();
-        dp[3] = w_thickness();
-        dp[4] = core_sld();
-        dp[5] = shell_sld();
-        dp[6] = n_pairs();
-        dp[7] = 0.0;
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = core_radius();
+  dp[2] = s_thickness();
+  dp[3] = w_thickness();
+  dp[4] = core_sld();
+  dp[5] = shell_sld();
+  dp[6] = n_pairs();
+  dp[7] = 0.0;
 
-        // Get the dispersion points for the core radius
-        vector<WeightPoint> weights_core_radius;
-        core_radius.get_weights(weights_core_radius);
+  // Get the dispersion points for the core radius
+  vector<WeightPoint> weights_core_radius;
+  core_radius.get_weights(weights_core_radius);
 
-        // Get the dispersion points for the s_thickness
-        vector<WeightPoint> weights_s_thickness;
-        s_thickness.get_weights(weights_s_thickness);
+  // Get the dispersion points for the s_thickness
+  vector<WeightPoint> weights_s_thickness;
+  s_thickness.get_weights(weights_s_thickness);
 
-        // Get the dispersion points for the w_thickness
-        vector<WeightPoint> weights_w_thickness;
-        w_thickness.get_weights(weights_w_thickness);
+  // Get the dispersion points for the w_thickness
+  vector<WeightPoint> weights_w_thickness;
+  w_thickness.get_weights(weights_w_thickness);
 
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double vol = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
 
-        // Loop over radius weight points
-        for(int i=0; i< (int)weights_core_radius.size(); i++) {
-                dp[1] = weights_core_radius[i].value;
-                for(int j=0; j< (int)weights_s_thickness.size(); j++){
-                        dp[2] = weights_s_thickness[j].value;
-                        for(int k=0; k< (int)weights_w_thickness.size(); k++){
-                                dp[3] = weights_w_thickness[k].value;
-                                //Un-normalize SphereForm by volume
-                                sum += weights_core_radius[i].weight*weights_s_thickness[j].weight
-                                        *weights_w_thickness[k].weight* MultiShell(dp, q)
-                                        *pow(weights_core_radius[i].value+dp[6]*weights_s_thickness[j].value+(dp[6]-1)*weights_w_thickness[k].value,3);
-                                //Find average volume
-                                vol += weights_core_radius[i].weight*weights_s_thickness[j].weight
-                                        *weights_w_thickness[k].weight
-                                        *pow(weights_core_radius[i].value+dp[6]*weights_s_thickness[j].value+(dp[6]-1)*weights_w_thickness[k].value,3);
-                                norm += weights_core_radius[i].weight*weights_s_thickness[j].weight
-                                        *weights_w_thickness[k].weight;
-                        }
-                }
-        }
-        if (vol != 0.0 && norm != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm);}
-        return sum/norm + background();
+  // Loop over radius weight points
+  for(int i=0; i< (int)weights_core_radius.size(); i++) {
+    dp[1] = weights_core_radius[i].value;
+    for(int j=0; j< (int)weights_s_thickness.size(); j++){
+      dp[2] = weights_s_thickness[j].value;
+      for(int k=0; k< (int)weights_w_thickness.size(); k++){
+        dp[3] = weights_w_thickness[k].value;
+        //Un-normalize SphereForm by volume
+        sum += weights_core_radius[i].weight*weights_s_thickness[j].weight
+            *weights_w_thickness[k].weight* MultiShell(dp, q)
+        *pow(weights_core_radius[i].value+dp[6]*weights_s_thickness[j].value+(dp[6]-1)*weights_w_thickness[k].value,3);
+        //Find average volume
+        vol += weights_core_radius[i].weight*weights_s_thickness[j].weight
+            *weights_w_thickness[k].weight
+            *pow(weights_core_radius[i].value+dp[6]*weights_s_thickness[j].value+(dp[6]-1)*weights_w_thickness[k].value,3);
+        norm += weights_core_radius[i].weight*weights_s_thickness[j].weight
+            *weights_w_thickness[k].weight;
+      }
+    }
+  }
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+  return sum/norm + background();
 }
 
@@ -115 +126 @@
  */
 double MultiShellModel :: operator()(double qx, double qy) {
-        double q = sqrt(qx*qx + qy*qy);
-        return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -127 +138 @@
  */
 double MultiShellModel :: evaluate_rphi(double q, double phi) {
-        return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 /**
@@ -134 +145 @@
  */
 double MultiShellModel :: calculate_ER() {
-        MultiShellParameters dp;
+  MultiShellParameters dp;
 
-        dp.core_radius     = core_radius();
-        dp.s_thickness  = s_thickness();
-        dp.w_thickness  = w_thickness();
-        dp.n_pairs = n_pairs();
+  dp.core_radius     = core_radius();
+  dp.s_thickness  = s_thickness();
+  dp.w_thickness  = w_thickness();
+  dp.n_pairs = n_pairs();
 
-        double rad_out = 0.0;
+  double rad_out = 0.0;
 
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        if (dp.n_pairs <= 0.0 ){
-                dp.n_pairs = 0.0;
-        }
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  if (dp.n_pairs <= 0.0 ){
+    dp.n_pairs = 0.0;
+  }
 
-        // Get the dispersion points for the core radius
-        vector<WeightPoint> weights_core_radius;
-        core_radius.get_weights(weights_core_radius);
+  // Get the dispersion points for the core radius
+  vector<WeightPoint> weights_core_radius;
+  core_radius.get_weights(weights_core_radius);
 
-        // Get the dispersion points for the s_thickness
-        vector<WeightPoint> weights_s_thickness;
-        s_thickness.get_weights(weights_s_thickness);
+  // Get the dispersion points for the s_thickness
+  vector<WeightPoint> weights_s_thickness;
+  s_thickness.get_weights(weights_s_thickness);
 
-        // Get the dispersion points for the w_thickness
-        vector<WeightPoint> weights_w_thickness;
-        w_thickness.get_weights(weights_w_thickness);
+  // Get the dispersion points for the w_thickness
+  vector<WeightPoint> weights_w_thickness;
+  w_thickness.get_weights(weights_w_thickness);
 
-        // Loop over major shell weight points
-        for(int i=0; i< (int)weights_s_thickness.size(); i++) {
-                dp.s_thickness = weights_s_thickness[i].value;
-                for(int j=0; j< (int)weights_w_thickness.size(); j++) {
-                        dp.w_thickness = weights_w_thickness[j].value;
-                        for(int k=0; k< (int)weights_core_radius.size(); k++) {
-                                dp.core_radius = weights_core_radius[k].value;
-                                sum += weights_s_thickness[i].weight*weights_w_thickness[j].weight
-                                        * weights_core_radius[k].weight*(dp.core_radius+dp.n_pairs*dp.s_thickness+(dp.n_pairs-1.0)*dp.w_thickness);
-                                norm += weights_s_thickness[i].weight*weights_w_thickness[j].weight* weights_core_radius[k].weight;
-                        }
-                }
-        }
-        if (norm != 0){
-                //return the averaged value
-                rad_out =  sum/norm;}
-        else{
-                //return normal value
-                rad_out = (dp.core_radius+dp.n_pairs*dp.s_thickness+(dp.n_pairs-1.0)*dp.w_thickness);}
+  // Loop over major shell weight points
+  for(int i=0; i< (int)weights_s_thickness.size(); i++) {
+    dp.s_thickness = weights_s_thickness[i].value;
+    for(int j=0; j< (int)weights_w_thickness.size(); j++) {
+      dp.w_thickness = weights_w_thickness[j].value;
+      for(int k=0; k< (int)weights_core_radius.size(); k++) {
+        dp.core_radius = weights_core_radius[k].value;
+        sum += weights_s_thickness[i].weight*weights_w_thickness[j].weight
+            * weights_core_radius[k].weight*(dp.core_radius+dp.n_pairs*dp.s_thickness+(dp.n_pairs-1.0)*dp.w_thickness);
+        norm += weights_s_thickness[i].weight*weights_w_thickness[j].weight* weights_core_radius[k].weight;
+      }
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    rad_out = (dp.core_radius+dp.n_pairs*dp.s_thickness+(dp.n_pairs-1.0)*dp.w_thickness);}
 
-        return rad_out;
+  return rad_out;
 }

sansmodels/src/c_models/parallelepiped.cpp

@@ -22 +22 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>

sansmodels/src/c_models/polygausscoil.cpp

@@ -21 +21 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "polygausscoil.h"
 
 extern "C" {
         #include "libTwoPhase.h"
-        #include "polygausscoil.h"
 }
 

sansmodels/src/c_models/sc.cpp

@@ -21 +21 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>

sansmodels/src/c_models/sphere.cpp

@@ -21 +21 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>

sansmodels/src/c_models/spheroid.cpp

@@ -18 +18 @@
  *   sansmodels/src/libigor
  *
- *      TODO: refactor so that we pull in the old sansmodels.c_extensions
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
+#include <stdlib.h>
 using namespace std;
+#include "spheroid.h"
 
 extern "C" {
-        #include "libCylinder.h"
-        #include "libStructureFactor.h"
-        #include "spheroid.h"
+#include "libCylinder.h"
+#include "libStructureFactor.h"
+}
+
+typedef struct {
+  double scale;
+  double equat_core;
+  double polar_core;
+  double equat_shell;
+  double polar_shell;
+  double sld_core;
+  double sld_shell;
+  double sld_solvent;
+  double background;
+  double axis_theta;
+  double axis_phi;
+
+} SpheroidParameters;
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the prolate
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double spheroid_analytical_2D_scaled(SpheroidParameters *pars, double q, double q_x, double q_y) {
+
+  double cyl_x, cyl_y, cyl_z;
+  double q_z;
+  double alpha, vol, cos_val;
+  double answer;
+  double Pi = 4.0*atan(1.0);
+  double sldcs,sldss;
+
+  //convert angle degree to radian
+  double theta = pars->axis_theta * Pi/180.0;
+  double phi = pars->axis_phi * Pi/180.0;
+
+
+  // ellipsoid orientation, the axis of the rotation is consistent with the ploar axis.
+  cyl_x = sin(theta) * cos(phi);
+  cyl_y = sin(theta) * sin(phi);
+  cyl_z = cos(theta);
+  //del sld
+  sldcs = pars->sld_core - pars->sld_shell;
+  sldss = pars->sld_shell- pars->sld_solvent;
+
+  // q vector
+  q_z = 0;
+
+  // Compute the angle btw vector q and the
+  // axis of the cylinder
+  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
+
+  // The following test should always pass
+  if (fabs(cos_val)>1.0) {
+    printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
+    return 0;
+  }
+
+  // Note: cos(alpha) = 0 and 1 will get an
+  // undefined value from CylKernel
+  alpha = acos( cos_val );
+
+  // Call the IGOR library function to get the kernel: MUST use gfn4 not gf2 because of the def of params.
+  answer = gfn4(cos_val,pars->equat_core,pars->polar_core,pars->equat_shell,pars->polar_shell,sldcs,sldss,q);
+  //It seems that it should be normalized somehow. How???
+
+  //normalize by cylinder volume
+  //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
+  vol = 4.0*Pi/3.0*pars->equat_shell*pars->equat_shell*pars->polar_shell;
+  answer /= vol;
+
+  //convert to [cm-1]
+  answer *= 1.0e8;
+
+  //Scale
+  answer *= pars->scale;
+
+  // add in the background
+  answer += pars->background;
+
+  return answer;
 }
 
 CoreShellEllipsoidModel :: CoreShellEllipsoidModel() {
-        scale      = Parameter(1.0);
-        equat_core     = Parameter(200.0, true);
-        equat_core.set_min(0.0);
-        polar_core     = Parameter(20.0, true);
-        polar_core.set_min(0.0);
-        equat_shell   = Parameter(250.0, true);
-        equat_shell.set_min(0.0);
-        polar_shell    = Parameter(30.0, true);
-        polar_shell.set_min(0.0);
-        sld_core   = Parameter(2e-6);
-        sld_shell  = Parameter(1e-6);
-        sld_solvent = Parameter(6.3e-6);
-        background = Parameter(0.0);
-        axis_theta  = Parameter(0.0, true);
-        axis_phi    = Parameter(0.0, true);
-
+  scale      = Parameter(1.0);
+  equat_core     = Parameter(200.0, true);
+  equat_core.set_min(0.0);
+  polar_core     = Parameter(20.0, true);
+  polar_core.set_min(0.0);
+  equat_shell   = Parameter(250.0, true);
+  equat_shell.set_min(0.0);
+  polar_shell    = Parameter(30.0, true);
+  polar_shell.set_min(0.0);
+  sld_core   = Parameter(2e-6);
+  sld_shell  = Parameter(1e-6);
+  sld_solvent = Parameter(6.3e-6);
+  background = Parameter(0.0);
+  axis_theta  = Parameter(0.0, true);
+  axis_phi    = Parameter(0.0, true);
+
+}
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the prolate
+ * @param q: q-value
+ * @return: function value
+ */
+static double spheroid_analytical_2DXY(SpheroidParameters *pars, double qx, double qy) {
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return spheroid_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
@@ -59 +154 @@
  */
 double CoreShellEllipsoidModel :: operator()(double q) {
-        double dp[9];
-
-        // Fill parameter array for IGOR library
-        // Add the background after averaging
-        dp[0] = scale();
-        dp[1] = equat_core();
-        dp[2] = polar_core();
-        dp[3] = equat_shell();
-        dp[4] = polar_shell();
-        dp[5] = sld_core();
-        dp[6] = sld_shell();
-        dp[7] = sld_solvent();
-        dp[8] = 0.0;
-
-        // Get the dispersion points for the major core
-        vector<WeightPoint> weights_equat_core;
-        equat_core.get_weights(weights_equat_core);
-
-        // Get the dispersion points for the minor core
-        vector<WeightPoint> weights_polar_core;
-        polar_core.get_weights(weights_polar_core);
-
-        // Get the dispersion points for the major shell
-        vector<WeightPoint> weights_equat_shell;
-        equat_shell.get_weights(weights_equat_shell);
-
-        // Get the dispersion points for the minor_shell
-        vector<WeightPoint> weights_polar_shell;
-        polar_shell.get_weights(weights_polar_shell);
-
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double vol = 0.0;
-
-        // Loop over major core weight points
-        for(int i=0; i<(int)weights_equat_core.size(); i++) {
-                dp[1] = weights_equat_core[i].value;
-
-                // Loop over minor core weight points
-                for(int j=0; j<(int)weights_polar_core.size(); j++) {
-                        dp[2] = weights_polar_core[j].value;
-
-                        // Loop over major shell weight points
-                        for(int k=0; k<(int)weights_equat_shell.size(); k++) {
-                                dp[3] = weights_equat_shell[k].value;
-
-                                // Loop over minor shell weight points
-                                for(int l=0; l<(int)weights_polar_shell.size(); l++) {
-                                        dp[4] = weights_polar_shell[l].value;
-                                        //Un-normalize  by volume
-                                        sum += weights_equat_core[i].weight* weights_polar_core[j].weight * weights_equat_shell[k].weight
-                                                * weights_polar_shell[l].weight * OblateForm(dp, q)
-                                                * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
-                                        //Find average volume
-                                        vol += weights_equat_core[i].weight* weights_polar_core[j].weight
-                                                * weights_equat_shell[k].weight
-                                                * weights_polar_shell[l].weight
-                                                * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
-                                        norm += weights_equat_core[i].weight* weights_polar_core[j].weight * weights_equat_shell[k].weight
-                                                        * weights_polar_shell[l].weight;
-                                }
-                        }
-                }
-        }
-        if (vol != 0.0 && norm != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm);}
-        return sum/norm + background();
+  double dp[9];
+
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = equat_core();
+  dp[2] = polar_core();
+  dp[3] = equat_shell();
+  dp[4] = polar_shell();
+  dp[5] = sld_core();
+  dp[6] = sld_shell();
+  dp[7] = sld_solvent();
+  dp[8] = 0.0;
+
+  // Get the dispersion points for the major core
+  vector<WeightPoint> weights_equat_core;
+  equat_core.get_weights(weights_equat_core);
+
+  // Get the dispersion points for the minor core
+  vector<WeightPoint> weights_polar_core;
+  polar_core.get_weights(weights_polar_core);
+
+  // Get the dispersion points for the major shell
+  vector<WeightPoint> weights_equat_shell;
+  equat_shell.get_weights(weights_equat_shell);
+
+  // Get the dispersion points for the minor_shell
+  vector<WeightPoint> weights_polar_shell;
+  polar_shell.get_weights(weights_polar_shell);
+
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over major core weight points
+  for(int i=0; i<(int)weights_equat_core.size(); i++) {
+    dp[1] = weights_equat_core[i].value;
+
+    // Loop over minor core weight points
+    for(int j=0; j<(int)weights_polar_core.size(); j++) {
+      dp[2] = weights_polar_core[j].value;
+
+      // Loop over major shell weight points
+      for(int k=0; k<(int)weights_equat_shell.size(); k++) {
+        dp[3] = weights_equat_shell[k].value;
+
+        // Loop over minor shell weight points
+        for(int l=0; l<(int)weights_polar_shell.size(); l++) {
+          dp[4] = weights_polar_shell[l].value;
+          //Un-normalize  by volume
+          sum += weights_equat_core[i].weight* weights_polar_core[j].weight * weights_equat_shell[k].weight
+              * weights_polar_shell[l].weight * OblateForm(dp, q)
+          * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
+          //Find average volume
+          vol += weights_equat_core[i].weight* weights_polar_core[j].weight
+              * weights_equat_shell[k].weight
+              * weights_polar_shell[l].weight
+              * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
+          norm += weights_equat_core[i].weight* weights_polar_core[j].weight * weights_equat_shell[k].weight
+              * weights_polar_shell[l].weight;
+        }
+      }
+    }
+  }
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+  return sum/norm + background();
 }
 
@@ -143 +238 @@
         return (*this).operator()(q);
 }
-*/
+ */
 
 /**
@@ -153 +248 @@
  */
 double CoreShellEllipsoidModel :: evaluate_rphi(double q, double phi) {
-        double qx = q*cos(phi);
-        double qy = q*sin(phi);
-        return (*this).operator()(qx, qy);
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
 }
 
@@ -165 +260 @@
  */
 double CoreShellEllipsoidModel :: operator()(double qx, double qy) {
-        SpheroidParameters dp;
-        // Fill parameter array
-        dp.scale      = scale();
-        dp.equat_core = equat_core();
-        dp.polar_core = polar_core();
-        dp.equat_shell = equat_shell();
-        dp.polar_shell = polar_shell();
-        dp.sld_core = sld_core();
-        dp.sld_shell = sld_shell();
-        dp.sld_solvent = sld_solvent();
-        dp.background = 0.0;
-        dp.axis_theta = axis_theta();
-        dp.axis_phi = axis_phi();
-
-        // Get the dispersion points for the major core
-        vector<WeightPoint> weights_equat_core;
-        equat_core.get_weights(weights_equat_core);
-
-        // Get the dispersion points for the minor core
-        vector<WeightPoint> weights_polar_core;
-        polar_core.get_weights(weights_polar_core);
-
-        // Get the dispersion points for the major shell
-        vector<WeightPoint> weights_equat_shell;
-        equat_shell.get_weights(weights_equat_shell);
-
-        // Get the dispersion points for the minor shell
-        vector<WeightPoint> weights_polar_shell;
-        polar_shell.get_weights(weights_polar_shell);
-
-
-        // Get angular averaging for theta
-        vector<WeightPoint> weights_theta;
-        axis_theta.get_weights(weights_theta);
-
-        // Get angular averaging for phi
-        vector<WeightPoint> weights_phi;
-        axis_phi.get_weights(weights_phi);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double norm_vol = 0.0;
-        double vol = 0.0;
-        double pi = 4.0*atan(1.0);
-        // Loop over major core weight points
-        for(int i=0; i< (int)weights_equat_core.size(); i++) {
-                dp.equat_core = weights_equat_core[i].value;
-
-                // Loop over minor core weight points
-                for(int j=0; j< (int)weights_polar_core.size(); j++) {
-                        dp.polar_core = weights_polar_core[j].value;
-
-                        // Loop over major shell weight points
-                        for(int k=0; k< (int)weights_equat_shell.size(); k++) {
-                                dp.equat_shell = weights_equat_shell[i].value;
-
-                                // Loop over minor shell weight points
-                                for(int l=0; l< (int)weights_polar_shell.size(); l++) {
-                                        dp.polar_shell = weights_polar_shell[l].value;
-
-                                        // Average over theta distribution
-                                        for(int m=0; m< (int)weights_theta.size(); m++) {
-                                                dp.axis_theta = weights_theta[m].value;
-
-                                                // Average over phi distribution
-                                                for(int n=0; n< (int)weights_phi.size(); n++) {
-                                                        dp.axis_phi = weights_phi[n].value;
-                                                        //Un-normalize by volume
-                                                        double _ptvalue = weights_equat_core[i].weight *weights_polar_core[j].weight
-                                                                * weights_equat_shell[k].weight * weights_polar_shell[l].weight
-                                                                * weights_theta[m].weight
-                                                                * weights_phi[n].weight
-                                                                * spheroid_analytical_2DXY(&dp, qx, qy)
-                                                                * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
-                                                        if (weights_theta.size()>1) {
-                                                                _ptvalue *= fabs(sin(weights_theta[m].value*pi/180.0));
-                                                        }
-                                                        sum += _ptvalue;
-                                                        //Find average volume
-                                                        vol += weights_equat_shell[k].weight
-                                                                * weights_polar_shell[l].weight
-                                                                * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
-                                                        //Find norm for volume
-                                                        norm_vol += weights_equat_shell[k].weight
-                                                                * weights_polar_shell[l].weight;
-
-                                                        norm += weights_equat_core[i].weight *weights_polar_core[j].weight
-                                                                * weights_equat_shell[k].weight * weights_polar_shell[l].weight
-                                                                * weights_theta[m].weight * weights_phi[n].weight;
-                                                }
-                                        }
-                                }
-                        }
-                }
-        }
-        // Averaging in theta needs an extra normalization
-        // factor to account for the sin(theta) term in the
-        // integration (see documentation).
-        if (weights_theta.size()>1) norm = norm / asin(1.0);
-
-        if (vol != 0.0 && norm_vol != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm_vol);}
-
-        return sum/norm + background();
+  SpheroidParameters dp;
+  // Fill parameter array
+  dp.scale      = scale();
+  dp.equat_core = equat_core();
+  dp.polar_core = polar_core();
+  dp.equat_shell = equat_shell();
+  dp.polar_shell = polar_shell();
+  dp.sld_core = sld_core();
+  dp.sld_shell = sld_shell();
+  dp.sld_solvent = sld_solvent();
+  dp.background = 0.0;
+  dp.axis_theta = axis_theta();
+  dp.axis_phi = axis_phi();
+
+  // Get the dispersion points for the major core
+  vector<WeightPoint> weights_equat_core;
+  equat_core.get_weights(weights_equat_core);
+
+  // Get the dispersion points for the minor core
+  vector<WeightPoint> weights_polar_core;
+  polar_core.get_weights(weights_polar_core);
+
+  // Get the dispersion points for the major shell
+  vector<WeightPoint> weights_equat_shell;
+  equat_shell.get_weights(weights_equat_shell);
+
+  // Get the dispersion points for the minor shell
+  vector<WeightPoint> weights_polar_shell;
+  polar_shell.get_weights(weights_polar_shell);
+
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_theta;
+  axis_theta.get_weights(weights_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_phi;
+  axis_phi.get_weights(weights_phi);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double norm_vol = 0.0;
+  double vol = 0.0;
+  double pi = 4.0*atan(1.0);
+  // Loop over major core weight points
+  for(int i=0; i< (int)weights_equat_core.size(); i++) {
+    dp.equat_core = weights_equat_core[i].value;
+
+    // Loop over minor core weight points
+    for(int j=0; j< (int)weights_polar_core.size(); j++) {
+      dp.polar_core = weights_polar_core[j].value;
+
+      // Loop over major shell weight points
+      for(int k=0; k< (int)weights_equat_shell.size(); k++) {
+        dp.equat_shell = weights_equat_shell[i].value;
+
+        // Loop over minor shell weight points
+        for(int l=0; l< (int)weights_polar_shell.size(); l++) {
+          dp.polar_shell = weights_polar_shell[l].value;
+
+          // Average over theta distribution
+          for(int m=0; m< (int)weights_theta.size(); m++) {
+            dp.axis_theta = weights_theta[m].value;
+
+            // Average over phi distribution
+            for(int n=0; n< (int)weights_phi.size(); n++) {
+              dp.axis_phi = weights_phi[n].value;
+              //Un-normalize by volume
+              double _ptvalue = weights_equat_core[i].weight *weights_polar_core[j].weight
+                  * weights_equat_shell[k].weight * weights_polar_shell[l].weight
+                  * weights_theta[m].weight
+                  * weights_phi[n].weight
+                  * spheroid_analytical_2DXY(&dp, qx, qy)
+              * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
+              if (weights_theta.size()>1) {
+                _ptvalue *= fabs(sin(weights_theta[m].value*pi/180.0));
+              }
+              sum += _ptvalue;
+              //Find average volume
+              vol += weights_equat_shell[k].weight
+                  * weights_polar_shell[l].weight
+                  * pow(weights_equat_shell[k].value,2)*weights_polar_shell[l].value;
+              //Find norm for volume
+              norm_vol += weights_equat_shell[k].weight
+                  * weights_polar_shell[l].weight;
+
+              norm += weights_equat_core[i].weight *weights_polar_core[j].weight
+                  * weights_equat_shell[k].weight * weights_polar_shell[l].weight
+                  * weights_theta[m].weight * weights_phi[n].weight;
+            }
+          }
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_theta.size()>1) norm = norm / asin(1.0);
+
+  if (vol != 0.0 && norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm_vol);}
+
+  return sum/norm + background();
 }
 
@@ -278 +373 @@
  */
 double CoreShellEllipsoidModel :: calculate_ER() {
-        SpheroidParameters dp;
-
-        dp.equat_shell = equat_shell();
-        dp.polar_shell = polar_shell();
-
-        double rad_out = 0.0;
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-
-        // Get the dispersion points for the major shell
-        vector<WeightPoint> weights_equat_shell;
-        equat_shell.get_weights(weights_equat_shell);
-
-        // Get the dispersion points for the minor shell
-        vector<WeightPoint> weights_polar_shell;
-        polar_shell.get_weights(weights_polar_shell);
-
-        // Loop over major shell weight points
-        for(int i=0; i< (int)weights_equat_shell.size(); i++) {
-                dp.equat_shell = weights_equat_shell[i].value;
-                for(int k=0; k< (int)weights_polar_shell.size(); k++) {
-                        dp.polar_shell = weights_polar_shell[k].value;
-                        //Note: output of "DiamEllip(dp.polar_shell,dp.equat_shell)" is DIAMETER.
-                        sum +=weights_equat_shell[i].weight
-                                * weights_polar_shell[k].weight*DiamEllip(dp.polar_shell,dp.equat_shell)/2.0;
-                        norm += weights_equat_shell[i].weight* weights_polar_shell[k].weight;
-                }
-        }
-        if (norm != 0){
-                //return the averaged value
-                rad_out =  sum/norm;}
-        else{
-                //return normal value
-                //Note: output of "DiamEllip(dp.polar_shell,dp.equat_shell)" is DIAMETER.
-                rad_out = DiamEllip(dp.polar_shell,dp.equat_shell)/2.0;}
-
-        return rad_out;
-}
+  SpheroidParameters dp;
+
+  dp.equat_shell = equat_shell();
+  dp.polar_shell = polar_shell();
+
+  double rad_out = 0.0;
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the major shell
+  vector<WeightPoint> weights_equat_shell;
+  equat_shell.get_weights(weights_equat_shell);
+
+  // Get the dispersion points for the minor shell
+  vector<WeightPoint> weights_polar_shell;
+  polar_shell.get_weights(weights_polar_shell);
+
+  // Loop over major shell weight points
+  for(int i=0; i< (int)weights_equat_shell.size(); i++) {
+    dp.equat_shell = weights_equat_shell[i].value;
+    for(int k=0; k< (int)weights_polar_shell.size(); k++) {
+      dp.polar_shell = weights_polar_shell[k].value;
+      //Note: output of "DiamEllip(dp.polar_shell,dp.equat_shell)" is DIAMETER.
+      sum +=weights_equat_shell[i].weight
+          * weights_polar_shell[k].weight*DiamEllip(dp.polar_shell,dp.equat_shell)/2.0;
+      norm += weights_equat_shell[i].weight* weights_polar_shell[k].weight;
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    //Note: output of "DiamEllip(dp.polar_shell,dp.equat_shell)" is DIAMETER.
+    rad_out = DiamEllip(dp.polar_shell,dp.equat_shell)/2.0;}
+
+  return rad_out;
+}

sansmodels/src/c_models/stackeddisks.cpp

@@ -23 +23 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "stacked_disks.h"
 
 extern "C" {
-        #include "libCylinder.h"
-        #include "libStructureFactor.h"
-        #include "stacked_disks.h"
+#include "libCylinder.h"
+#include "libStructureFactor.h"
+}
+
+typedef struct {
+  double scale;
+  double radius;
+  double core_thick;
+  double layer_thick;
+  double core_sld;
+  double layer_sld;
+  double solvent_sld;
+  double n_stacking;
+  double sigma_d;
+  double background;
+  double axis_theta;
+  double axis_phi;
+} StackedDisksParameters;
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the staked disks
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double stacked_disks_analytical_2D_scaled(StackedDisksParameters *pars, double q, double q_x, double q_y) {
+  double cyl_x, cyl_y, cyl_z;
+  double q_z;
+  double alpha, vol, cos_val;
+  double d, dum, halfheight;
+  double answer;
+
+
+
+  // parallelepiped orientation
+  cyl_x = sin(pars->axis_theta) * cos(pars->axis_phi);
+  cyl_y = sin(pars->axis_theta) * sin(pars->axis_phi);
+  cyl_z = cos(pars->axis_theta);
+
+  // q vector
+  q_z = 0;
+
+  // Compute the angle btw vector q and the
+  // axis of the parallelepiped
+  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
+
+  // The following test should always pass
+  if (fabs(cos_val)>1.0) {
+    printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+    return 0;
+  }
+
+  // Note: cos(alpha) = 0 and 1 will get an
+  // undefined value from Stackdisc_kern
+  alpha = acos( cos_val );
+
+  // Call the IGOR library function to get the kernel
+  d = 2 * pars->layer_thick + pars->core_thick;
+  halfheight = pars->core_thick/2.0;
+  dum =alpha ;
+  answer = Stackdisc_kern(q, pars->radius, pars->core_sld,pars->layer_sld,
+      pars->solvent_sld, halfheight, pars->layer_thick, dum, pars->sigma_d, d, pars->n_stacking)/sin(alpha);
+
+  // Multiply by contrast^2
+  //answer *= pars->contrast*pars->contrast;
+
+  //normalize by staked disks volume
+  vol = acos(-1.0) * pars->radius * pars->radius * d * pars->n_stacking;
+  answer /= vol;
+
+  //convert to [cm-1]
+  answer *= 1.0e8;
+
+  //Scale
+  answer *= pars->scale;
+
+  // add in the background
+  answer += pars->background;
+
+  return answer;
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the staked disks
+ * @param q: q-value
+ * @return: function value
+ */
+static double stacked_disks_analytical_2DXY(StackedDisksParameters *pars, double qx, double qy) {
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return stacked_disks_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 
 StackedDisksModel :: StackedDisksModel() {
-        scale      = Parameter(1.0);
-        radius     = Parameter(3000.0, true);
-        radius.set_min(0.0);
-        core_thick  = Parameter(10.0, true);
-        core_thick.set_min(0.0);
-        layer_thick     = Parameter(15.0);
-        layer_thick.set_min(0.0);
-        core_sld = Parameter(4.0e-6);
-        layer_sld  = Parameter(-4.0e-7);
-        solvent_sld  = Parameter(5.0e-6);
-        n_stacking   = Parameter(1);
-        sigma_d   = Parameter(0);
-        background = Parameter(0.001);
-        axis_theta  = Parameter(0.0, true);
-        axis_phi    = Parameter(0.0, true);
+  scale      = Parameter(1.0);
+  radius     = Parameter(3000.0, true);
+  radius.set_min(0.0);
+  core_thick  = Parameter(10.0, true);
+  core_thick.set_min(0.0);
+  layer_thick     = Parameter(15.0);
+  layer_thick.set_min(0.0);
+  core_sld = Parameter(4.0e-6);
+  layer_sld  = Parameter(-4.0e-7);
+  solvent_sld  = Parameter(5.0e-6);
+  n_stacking   = Parameter(1);
+  sigma_d   = Parameter(0);
+  background = Parameter(0.001);
+  axis_theta  = Parameter(0.0, true);
+  axis_phi    = Parameter(0.0, true);
 }
 
@@ -59 +151 @@
  */
 double StackedDisksModel :: operator()(double q) {
-        double dp[10];
-
-        // Fill parameter array for IGOR library
-        // Add the background after averaging
-        dp[0] = scale();
-        dp[1] = radius();
-        dp[2] = core_thick();
-        dp[3] = layer_thick();
-        dp[4] = core_sld();
-        dp[5] = layer_sld();
-        dp[6] = solvent_sld();
-        dp[7] = n_stacking();
-        dp[8] = sigma_d();
-        dp[9] = 0.0;
-
-        // Get the dispersion points for the radius
-        vector<WeightPoint> weights_radius;
-        radius.get_weights(weights_radius);
-
-        // Get the dispersion points for the core_thick
-        vector<WeightPoint> weights_core_thick;
-        core_thick.get_weights(weights_core_thick);
-
-        // Get the dispersion points for the layer_thick
-        vector<WeightPoint> weights_layer_thick;
-        layer_thick.get_weights(weights_layer_thick);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double vol = 0.0;
-
-        // Loop over length weight points
-        for(int i=0; i< (int)weights_radius.size(); i++) {
-                dp[1] = weights_radius[i].value;
-
-                // Loop over radius weight points
-                for(int j=0; j< (int)weights_core_thick.size(); j++) {
-                        dp[2] = weights_core_thick[j].value;
-
-                        // Loop over thickness weight points
-                        for(int k=0; k< (int)weights_layer_thick.size(); k++) {
-                                dp[3] = weights_layer_thick[k].value;
-                                //Un-normalize by volume
-                                sum += weights_radius[i].weight
-                                        * weights_core_thick[j].weight * weights_layer_thick[k].weight* StackedDiscs(dp, q)
-                                        *pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value);
-                                //Find average volume
-                                vol += weights_radius[i].weight
-                                        * weights_core_thick[j].weight * weights_layer_thick[k].weight
-                                        *pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value);
-                                norm += weights_radius[i].weight
-                                        * weights_core_thick[j].weight* weights_layer_thick[k].weight;
-                        }
-                }
-        }
-        if (vol != 0.0 && norm != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm);}
-
-        return sum/norm + background();
+  double dp[10];
+
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = radius();
+  dp[2] = core_thick();
+  dp[3] = layer_thick();
+  dp[4] = core_sld();
+  dp[5] = layer_sld();
+  dp[6] = solvent_sld();
+  dp[7] = n_stacking();
+  dp[8] = sigma_d();
+  dp[9] = 0.0;
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_radius;
+  radius.get_weights(weights_radius);
+
+  // Get the dispersion points for the core_thick
+  vector<WeightPoint> weights_core_thick;
+  core_thick.get_weights(weights_core_thick);
+
+  // Get the dispersion points for the layer_thick
+  vector<WeightPoint> weights_layer_thick;
+  layer_thick.get_weights(weights_layer_thick);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over length weight points
+  for(int i=0; i< (int)weights_radius.size(); i++) {
+    dp[1] = weights_radius[i].value;
+
+    // Loop over radius weight points
+    for(int j=0; j< (int)weights_core_thick.size(); j++) {
+      dp[2] = weights_core_thick[j].value;
+
+      // Loop over thickness weight points
+      for(int k=0; k< (int)weights_layer_thick.size(); k++) {
+        dp[3] = weights_layer_thick[k].value;
+        //Un-normalize by volume
+        sum += weights_radius[i].weight
+            * weights_core_thick[j].weight * weights_layer_thick[k].weight* StackedDiscs(dp, q)
+        *pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value);
+        //Find average volume
+        vol += weights_radius[i].weight
+            * weights_core_thick[j].weight * weights_layer_thick[k].weight
+            *pow(weights_radius[i].value,2)*(weights_core_thick[j].value+2*weights_layer_thick[k].value);
+        norm += weights_radius[i].weight
+            * weights_core_thick[j].weight* weights_layer_thick[k].weight;
+      }
+    }
+  }
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -129 +221 @@
  */
 double StackedDisksModel :: operator()(double qx, double qy) {
-        StackedDisksParameters dp;
-        // Fill parameter array
-        dp.scale      = scale();
-        dp.core_thick    = core_thick();
-        dp.radius         = radius();
-        dp.layer_thick  = layer_thick();
-        dp.core_sld   = core_sld();
-        dp.layer_sld  = layer_sld();
-        dp.solvent_sld= solvent_sld();
-        dp.n_stacking     = n_stacking();
-        dp.sigma_d   = sigma_d();
-        dp.background = 0.0;
-        dp.axis_theta = axis_theta();
-        dp.axis_phi   = axis_phi();
-
-        // Get the dispersion points for the length
-        vector<WeightPoint> weights_core_thick;
-        core_thick.get_weights(weights_core_thick);
-
-        // Get the dispersion points for the radius
-        vector<WeightPoint> weights_radius;
-        radius.get_weights(weights_radius);
-
-        // Get the dispersion points for the thickness
-        vector<WeightPoint> weights_layer_thick;
-        layer_thick.get_weights(weights_layer_thick);
-
-        // Get angular averaging for theta
-        vector<WeightPoint> weights_theta;
-        axis_theta.get_weights(weights_theta);
-
-        // Get angular averaging for phi
-        vector<WeightPoint> weights_phi;
-        axis_phi.get_weights(weights_phi);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double norm_vol = 0.0;
-        double vol = 0.0;
-
-        // Loop over length weight points
-        for(int i=0; i< (int)weights_core_thick.size(); i++) {
-                dp.core_thick = weights_core_thick[i].value;
-
-                // Loop over radius weight points
-                for(int j=0; j< (int)weights_radius.size(); j++) {
-                        dp.radius = weights_radius[j].value;
-
-                                // Loop over thickness weight points
-                                for(int k=0; k< (int)weights_layer_thick.size(); k++) {
-                                dp.layer_thick = weights_layer_thick[k].value;
-
-                                        for(int l=0; l< (int)weights_theta.size(); l++) {
-                                        dp.axis_theta = weights_theta[l].value;
-
-                                                // Average over phi distribution
-                                                for(int m=0; m <(int)weights_phi.size(); m++) {
-                                                        dp.axis_phi = weights_phi[m].value;
-
-                                                        //Un-normalize by volume
-                                                        double _ptvalue = weights_core_thick[i].weight
-                                                                * weights_radius[j].weight
-                                                                * weights_layer_thick[k].weight
-                                                                * weights_theta[l].weight
-                                                                * weights_phi[m].weight
-                                                                * stacked_disks_analytical_2DXY(&dp, qx, qy)
-                                                                *pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value);
-                                                        if (weights_theta.size()>1) {
-                                                                _ptvalue *= fabs(sin(weights_theta[l].value));
-                                                        }
-                                                        sum += _ptvalue;
-                                                        //Find average volume
-                                                        vol += weights_radius[j].weight
-                                                                * weights_core_thick[i].weight * weights_layer_thick[k].weight
-                                                                *pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value);
-                                                        //Find norm for volume
-                                                        norm_vol += weights_radius[j].weight
-                                                                * weights_core_thick[i].weight * weights_layer_thick[k].weight;
-
-                                                        norm += weights_core_thick[i].weight
-                                                                * weights_radius[j].weight
-                                                                * weights_layer_thick[k].weight
-                                                                * weights_theta[l].weight
-                                                                * weights_phi[m].weight;
-                                                }
-                                }
-                        }
-                }
-        }
-        // Averaging in theta needs an extra normalization
-        // factor to account for the sin(theta) term in the
-        // integration (see documentation).
-        if (weights_theta.size()>1) norm = norm / asin(1.0);
-        if (vol != 0.0 && norm_vol != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm_vol);}
-        return sum/norm + background();
+  StackedDisksParameters dp;
+  // Fill parameter array
+  dp.scale      = scale();
+  dp.core_thick    = core_thick();
+  dp.radius       = radius();
+  dp.layer_thick  = layer_thick();
+  dp.core_sld   = core_sld();
+  dp.layer_sld  = layer_sld();
+  dp.solvent_sld= solvent_sld();
+  dp.n_stacking   = n_stacking();
+  dp.sigma_d   = sigma_d();
+  dp.background = 0.0;
+  dp.axis_theta = axis_theta();
+  dp.axis_phi   = axis_phi();
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_core_thick;
+  core_thick.get_weights(weights_core_thick);
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_radius;
+  radius.get_weights(weights_radius);
+
+  // Get the dispersion points for the thickness
+  vector<WeightPoint> weights_layer_thick;
+  layer_thick.get_weights(weights_layer_thick);
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_theta;
+  axis_theta.get_weights(weights_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_phi;
+  axis_phi.get_weights(weights_phi);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double norm_vol = 0.0;
+  double vol = 0.0;
+
+  // Loop over length weight points
+  for(int i=0; i< (int)weights_core_thick.size(); i++) {
+    dp.core_thick = weights_core_thick[i].value;
+
+    // Loop over radius weight points
+    for(int j=0; j< (int)weights_radius.size(); j++) {
+      dp.radius = weights_radius[j].value;
+
+      // Loop over thickness weight points
+      for(int k=0; k< (int)weights_layer_thick.size(); k++) {
+        dp.layer_thick = weights_layer_thick[k].value;
+
+        for(int l=0; l< (int)weights_theta.size(); l++) {
+          dp.axis_theta = weights_theta[l].value;
+
+          // Average over phi distribution
+          for(int m=0; m <(int)weights_phi.size(); m++) {
+            dp.axis_phi = weights_phi[m].value;
+
+            //Un-normalize by volume
+            double _ptvalue = weights_core_thick[i].weight
+                * weights_radius[j].weight
+                * weights_layer_thick[k].weight
+                * weights_theta[l].weight
+                * weights_phi[m].weight
+                * stacked_disks_analytical_2DXY(&dp, qx, qy)
+            *pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value);
+            if (weights_theta.size()>1) {
+              _ptvalue *= fabs(sin(weights_theta[l].value));
+            }
+            sum += _ptvalue;
+            //Find average volume
+            vol += weights_radius[j].weight
+                * weights_core_thick[i].weight * weights_layer_thick[k].weight
+                *pow(weights_radius[j].value,2)*(weights_core_thick[i].value+2*weights_layer_thick[k].value);
+            //Find norm for volume
+            norm_vol += weights_radius[j].weight
+                * weights_core_thick[i].weight * weights_layer_thick[k].weight;
+
+            norm += weights_core_thick[i].weight
+                * weights_radius[j].weight
+                * weights_layer_thick[k].weight
+                * weights_theta[l].weight
+                * weights_phi[m].weight;
+          }
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_theta.size()>1) norm = norm / asin(1.0);
+  if (vol != 0.0 && norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm_vol);}
+  return sum/norm + background();
 }
 
@@ -237 +329 @@
  */
 double StackedDisksModel :: evaluate_rphi(double q, double phi) {
-        double qx = q*cos(phi);
-        double qy = q*sin(phi);
-        return (*this).operator()(qx, qy);
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
 }
 /**
@@ -246 +338 @@
  */
 double StackedDisksModel :: calculate_ER() {
-        StackedDisksParameters dp;
-
-        dp.core_thick    = core_thick();
-        dp.radius         = radius();
-        dp.layer_thick  = layer_thick();
-        dp.n_stacking     = n_stacking();
-
-        double rad_out = 0.0;
-        if (dp.n_stacking <= 0.0){
-                return rad_out;
-        }
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-
-        // Get the dispersion points for the length
-        vector<WeightPoint> weights_core_thick;
-        core_thick.get_weights(weights_core_thick);
-
-        // Get the dispersion points for the radius
-        vector<WeightPoint> weights_radius;
-        radius.get_weights(weights_radius);
-
-        // Get the dispersion points for the thickness
-        vector<WeightPoint> weights_layer_thick;
-        layer_thick.get_weights(weights_layer_thick);
-
-        // Loop over major shell weight points
-        for(int i=0; i< (int)weights_core_thick.size(); i++) {
-                dp.core_thick = weights_core_thick[i].value;
-                for(int j=0; j< (int)weights_layer_thick.size(); j++) {
-                        dp.layer_thick = weights_layer_thick[j].value;
-                        for(int k=0; k< (int)weights_radius.size(); k++) {
-                                dp.radius = weights_radius[k].value;
-                                //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-                                sum +=weights_core_thick[i].weight*weights_layer_thick[j].weight
-                                        * weights_radius[k].weight*DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0;
-                                norm += weights_core_thick[i].weight*weights_layer_thick[j].weight* weights_radius[k].weight;
-                        }
-                }
-        }
-        if (norm != 0){
-                //return the averaged value
-                rad_out =  sum/norm;}
-        else{
-                //return normal value
-                //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
-                rad_out = DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0;}
-
-        return rad_out;
-}
+  StackedDisksParameters dp;
+
+  dp.core_thick    = core_thick();
+  dp.radius       = radius();
+  dp.layer_thick  = layer_thick();
+  dp.n_stacking   = n_stacking();
+
+  double rad_out = 0.0;
+  if (dp.n_stacking <= 0.0){
+    return rad_out;
+  }
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the length
+  vector<WeightPoint> weights_core_thick;
+  core_thick.get_weights(weights_core_thick);
+
+  // Get the dispersion points for the radius
+  vector<WeightPoint> weights_radius;
+  radius.get_weights(weights_radius);
+
+  // Get the dispersion points for the thickness
+  vector<WeightPoint> weights_layer_thick;
+  layer_thick.get_weights(weights_layer_thick);
+
+  // Loop over major shell weight points
+  for(int i=0; i< (int)weights_core_thick.size(); i++) {
+    dp.core_thick = weights_core_thick[i].value;
+    for(int j=0; j< (int)weights_layer_thick.size(); j++) {
+      dp.layer_thick = weights_layer_thick[j].value;
+      for(int k=0; k< (int)weights_radius.size(); k++) {
+        dp.radius = weights_radius[k].value;
+        //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+        sum +=weights_core_thick[i].weight*weights_layer_thick[j].weight
+            * weights_radius[k].weight*DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0;
+        norm += weights_core_thick[i].weight*weights_layer_thick[j].weight* weights_radius[k].weight;
+      }
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    //Note: output of "DiamCyl(dp.length,dp.radius)" is DIAMETER.
+    rad_out = DiamCyl(dp.n_stacking*(dp.layer_thick*2.0+dp.core_thick),dp.radius)/2.0;}
+
+  return rad_out;
+}

sansmodels/src/c_models/triaxialellipsoid.cpp

@@ -18 +18 @@
  *   sansmodels/src/libigor
  *
- *      TODO: refactor so that we pull in the old sansmodels.c_extensions
  */
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
+#include <stdlib.h>
 using namespace std;
+#include "triaxial_ellipsoid.h"
 
 extern "C" {
-        #include "libCylinder.h"
-        #include "libStructureFactor.h"
-        #include "triaxial_ellipsoid.h"
-}
+#include "libCylinder.h"
+#include "libStructureFactor.h"
+}
+
+typedef struct {
+  double scale;
+  double semi_axisA;
+  double semi_axisB;
+  double semi_axisC;
+  double sldEll;
+  double sldSolv;
+  double background;
+  double axis_theta;
+  double axis_phi;
+  double axis_psi;
+
+} TriaxialEllipsoidParameters;
+
+static double triaxial_ellipsoid_kernel(TriaxialEllipsoidParameters *pars, double q, double alpha, double nu) {
+  double t,a,b,c;
+  double kernel;
+
+  a = pars->semi_axisA ;
+  b = pars->semi_axisB ;
+  c = pars->semi_axisC ;
+
+  t = q * sqrt(a*a*cos(nu)*cos(nu)+b*b*sin(nu)*sin(nu)*sin(alpha)*sin(alpha)+c*c*cos(alpha)*cos(alpha));
+  if (t==0.0){
+    kernel  = 1.0;
+  }else{
+    kernel  = 3.0*(sin(t)-t*cos(t))/(t*t*t);
+  }
+  return kernel*kernel;
+}
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the triaxial ellipsoid
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double triaxial_ellipsoid_analytical_2D_scaled(TriaxialEllipsoidParameters *pars, double q, double q_x, double q_y) {
+  double cyl_x, cyl_y, cyl_z, ell_x, ell_y;
+  double q_z;
+  double cos_nu,nu;
+  double alpha, vol, cos_val;
+  double answer;
+  double pi = 4.0*atan(1.0);
+
+  //convert angle degree to radian
+  double theta = pars->axis_theta * pi/180.0;
+  double phi = pars->axis_phi * pi/180.0;
+  double psi = pars->axis_psi * pi/180.0;
+
+  // Cylinder orientation
+  cyl_x = sin(theta) * cos(phi);
+  cyl_y = sin(theta) * sin(phi);
+  cyl_z = cos(theta);
+
+  // q vector
+  q_z = 0.0;
+
+  //dx = 1.0;
+  //dy = 1.0;
+  // Compute the angle btw vector q and the
+  // axis of the cylinder
+  cos_val = cyl_x*q_x + cyl_y*q_y + cyl_z*q_z;
+
+  // The following test should always pass
+  if (fabs(cos_val)>1.0) {
+    printf("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
+    return 0;
+  }
+
+  // Note: cos(alpha) = 0 and 1 will get an
+  // undefined value from CylKernel
+  alpha = acos( cos_val );
+
+  //ellipse orientation:
+  // the elliptical corss section was transformed and projected
+  // into the detector plane already through sin(alpha)and furthermore psi remains as same
+  // on the detector plane.
+  // So, all we need is to calculate the angle (nu) of the minor axis of the ellipse wrt
+  // the wave vector q.
+
+  //x- y- component on the detector plane.
+  ell_x =  cos(psi);
+  ell_y =  sin(psi);
+
+  // calculate the axis of the ellipse wrt q-coord.
+  cos_nu = ell_x*q_x + ell_y*q_y;
+  nu = acos(cos_nu);
+
+  // Call the IGOR library function to get the kernel
+  answer = triaxial_ellipsoid_kernel(pars, q, alpha, nu);
+
+  // Multiply by contrast^2
+  answer *= (pars->sldEll- pars->sldSolv)*(pars->sldEll- pars->sldSolv);
+
+  //normalize by cylinder volume
+  //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
+  vol = 4.0* pi/3.0  * pars->semi_axisA * pars->semi_axisB * pars->semi_axisC;
+  answer *= vol;
+  //convert to [cm-1]
+  answer *= 1.0e8;
+  //Scale
+  answer *= pars->scale;
+
+  // add in the background
+  answer += pars->background;
+
+  return answer;
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the triaxial ellipsoid
+ * @param q: q-value
+ * @return: function value
+ */
+static double triaxial_ellipsoid_analytical_2DXY(TriaxialEllipsoidParameters *pars, double qx, double qy) {
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return triaxial_ellipsoid_analytical_2D_scaled(pars, q, qx/q, qy/q);
+}
+
+
 
 TriaxialEllipsoidModel :: TriaxialEllipsoidModel() {
-        scale      = Parameter(1.0);
-        semi_axisA     = Parameter(35.0, true);
-        semi_axisA.set_min(0.0);
-        semi_axisB     = Parameter(100.0, true);
-        semi_axisB.set_min(0.0);
-        semi_axisC  = Parameter(400.0, true);
-        semi_axisC.set_min(0.0);
-        sldEll   = Parameter(1.0e-6);
-        sldSolv   = Parameter(6.3e-6);
-        background = Parameter(0.0);
-        axis_theta  = Parameter(57.325, true);
-        axis_phi    = Parameter(57.325, true);
-        axis_psi    = Parameter(0.0, true);
+  scale      = Parameter(1.0);
+  semi_axisA     = Parameter(35.0, true);
+  semi_axisA.set_min(0.0);
+  semi_axisB     = Parameter(100.0, true);
+  semi_axisB.set_min(0.0);
+  semi_axisC  = Parameter(400.0, true);
+  semi_axisC.set_min(0.0);
+  sldEll   = Parameter(1.0e-6);
+  sldSolv   = Parameter(6.3e-6);
+  background = Parameter(0.0);
+  axis_theta  = Parameter(57.325, true);
+  axis_phi    = Parameter(57.325, true);
+  axis_psi    = Parameter(0.0, true);
 }
 
@@ -56 +182 @@
  */
 double TriaxialEllipsoidModel :: operator()(double q) {
-        double dp[7];
-
-        // Fill parameter array for IGOR library
-        // Add the background after averaging
-        dp[0] = scale();
-        dp[1] = semi_axisA();
-        dp[2] = semi_axisB();
-        dp[3] = semi_axisC();
-        dp[4] = sldEll();
-        dp[5] = sldSolv();
-        dp[6] = 0.0;
-
-        // Get the dispersion points for the semi axis A
-        vector<WeightPoint> weights_semi_axisA;
-        semi_axisA.get_weights(weights_semi_axisA);
-
-        // Get the dispersion points for the semi axis B
-        vector<WeightPoint> weights_semi_axisB;
-        semi_axisB.get_weights(weights_semi_axisB);
-
-        // Get the dispersion points for the semi axis C
-        vector<WeightPoint> weights_semi_axisC;
-        semi_axisC.get_weights(weights_semi_axisC);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double vol = 0.0;
-
-        // Loop over semi axis A weight points
-        for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
-                dp[1] = weights_semi_axisA[i].value;
-
-                // Loop over semi axis B weight points
-                for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
-                        dp[2] = weights_semi_axisB[j].value;
-
-                        // Loop over semi axis C weight points
-                        for(int k=0; k< (int)weights_semi_axisC.size(); k++) {
-                                dp[3] = weights_semi_axisC[k].value;
-                                //Un-normalize  by volume
-                                sum += weights_semi_axisA[i].weight
-                                        * weights_semi_axisB[j].weight * weights_semi_axisC[k].weight* TriaxialEllipsoid(dp, q)
-                                        * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
-                                //Find average volume
-                                vol += weights_semi_axisA[i].weight
-                                        * weights_semi_axisB[j].weight * weights_semi_axisC[k].weight
-                                        * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
-
-                                norm += weights_semi_axisA[i].weight
-                                        * weights_semi_axisB[j].weight * weights_semi_axisC[k].weight;
-                        }
-                }
-        }
-        if (vol != 0.0 && norm != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm);}
-
-        return sum/norm + background();
+  double dp[7];
+
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = semi_axisA();
+  dp[2] = semi_axisB();
+  dp[3] = semi_axisC();
+  dp[4] = sldEll();
+  dp[5] = sldSolv();
+  dp[6] = 0.0;
+
+  // Get the dispersion points for the semi axis A
+  vector<WeightPoint> weights_semi_axisA;
+  semi_axisA.get_weights(weights_semi_axisA);
+
+  // Get the dispersion points for the semi axis B
+  vector<WeightPoint> weights_semi_axisB;
+  semi_axisB.get_weights(weights_semi_axisB);
+
+  // Get the dispersion points for the semi axis C
+  vector<WeightPoint> weights_semi_axisC;
+  semi_axisC.get_weights(weights_semi_axisC);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
+
+  // Loop over semi axis A weight points
+  for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
+    dp[1] = weights_semi_axisA[i].value;
+
+    // Loop over semi axis B weight points
+    for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
+      dp[2] = weights_semi_axisB[j].value;
+
+      // Loop over semi axis C weight points
+      for(int k=0; k< (int)weights_semi_axisC.size(); k++) {
+        dp[3] = weights_semi_axisC[k].value;
+        //Un-normalize  by volume
+        sum += weights_semi_axisA[i].weight
+            * weights_semi_axisB[j].weight * weights_semi_axisC[k].weight* TriaxialEllipsoid(dp, q)
+        * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
+        //Find average volume
+        vol += weights_semi_axisA[i].weight
+            * weights_semi_axisB[j].weight * weights_semi_axisC[k].weight
+            * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
+
+        norm += weights_semi_axisA[i].weight
+            * weights_semi_axisB[j].weight * weights_semi_axisC[k].weight;
+      }
+    }
+  }
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
+
+  return sum/norm + background();
 }
 
@@ -124 +250 @@
  */
 double TriaxialEllipsoidModel :: operator()(double qx, double qy) {
-        TriaxialEllipsoidParameters dp;
-        // Fill parameter array
-        dp.scale      = scale();
-        dp.semi_axisA   = semi_axisA();
-        dp.semi_axisB     = semi_axisB();
-        dp.semi_axisC     = semi_axisC();
-        dp.sldEll   = sldEll();
-        dp.sldSolv   = sldSolv();
-        dp.background = 0.0;
-        dp.axis_theta  = axis_theta();
-        dp.axis_phi    = axis_phi();
-        dp.axis_psi    = axis_psi();
-
-        // Get the dispersion points for the semi_axis A
-        vector<WeightPoint> weights_semi_axisA;
-        semi_axisA.get_weights(weights_semi_axisA);
-
-        // Get the dispersion points for the semi_axis B
-        vector<WeightPoint> weights_semi_axisB;
-        semi_axisB.get_weights(weights_semi_axisB);
-
-        // Get the dispersion points for the semi_axis C
-        vector<WeightPoint> weights_semi_axisC;
-        semi_axisC.get_weights(weights_semi_axisC);
-
-        // Get angular averaging for theta
-        vector<WeightPoint> weights_theta;
-        axis_theta.get_weights(weights_theta);
-
-        // Get angular averaging for phi
-        vector<WeightPoint> weights_phi;
-        axis_phi.get_weights(weights_phi);
-
-        // Get angular averaging for psi
-        vector<WeightPoint> weights_psi;
-        axis_psi.get_weights(weights_psi);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double norm_vol = 0.0;
-        double vol = 0.0;
-        double pi = 4.0*atan(1.0);
-        // Loop over semi axis A weight points
-        for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
-                dp.semi_axisA = weights_semi_axisA[i].value;
-
-                // Loop over semi axis B weight points
-                for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
-                        dp.semi_axisB = weights_semi_axisB[j].value;
-
-                        // Loop over semi axis C weight points
-                        for(int k=0; k < (int)weights_semi_axisC.size(); k++) {
-                        dp.semi_axisC = weights_semi_axisC[k].value;
-
-                                // Average over theta distribution
-                                for(int l=0; l< (int)weights_theta.size(); l++) {
-                                        dp.axis_theta = weights_theta[l].value;
-
-                                        // Average over phi distribution
-                                        for(int m=0; m <(int)weights_phi.size(); m++) {
-                                                dp.axis_phi = weights_phi[m].value;
-                                                // Average over psi distribution
-                                                for(int n=0; n <(int)weights_psi.size(); n++) {
-                                                        dp.axis_psi = weights_psi[n].value;
-                                                        //Un-normalize  by volume
-                                                        double _ptvalue = weights_semi_axisA[i].weight
-                                                                * weights_semi_axisB[j].weight
-                                                                * weights_semi_axisC[k].weight
-                                                                * weights_theta[l].weight
-                                                                * weights_phi[m].weight
-                                                                * weights_psi[n].weight
-                                                                * triaxial_ellipsoid_analytical_2DXY(&dp, qx, qy)
-                                                                * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
-                                                        if (weights_theta.size()>1) {
-                                                                _ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
-                                                        }
-                                                        sum += _ptvalue;
-                                                        //Find average volume
-                                                        vol += weights_semi_axisA[i].weight
-                                                                * weights_semi_axisB[j].weight
-                                                                * weights_semi_axisC[k].weight
-                                                                * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
-                                                        //Find norm for volume
-                                                        norm_vol += weights_semi_axisA[i].weight
-                                                                * weights_semi_axisB[j].weight
-                                                                * weights_semi_axisC[k].weight;
-
-                                                        norm += weights_semi_axisA[i].weight
-                                                                * weights_semi_axisB[j].weight
-                                                                * weights_semi_axisC[k].weight
-                                                                * weights_theta[l].weight
-                                                                * weights_phi[m].weight
-                                                                * weights_psi[n].weight;
-                                                }
-                                        }
-
-                                }
-                        }
-                }
-        }
-        // Averaging in theta needs an extra normalization
-        // factor to account for the sin(theta) term in the
-        // integration (see documentation).
-        if (weights_theta.size()>1) norm = norm / asin(1.0);
-
-        if (vol != 0.0 && norm_vol != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm_vol);}
-
-        return sum/norm + background();
+  TriaxialEllipsoidParameters dp;
+  // Fill parameter array
+  dp.scale      = scale();
+  dp.semi_axisA   = semi_axisA();
+  dp.semi_axisB     = semi_axisB();
+  dp.semi_axisC     = semi_axisC();
+  dp.sldEll   = sldEll();
+  dp.sldSolv   = sldSolv();
+  dp.background = 0.0;
+  dp.axis_theta  = axis_theta();
+  dp.axis_phi    = axis_phi();
+  dp.axis_psi    = axis_psi();
+
+  // Get the dispersion points for the semi_axis A
+  vector<WeightPoint> weights_semi_axisA;
+  semi_axisA.get_weights(weights_semi_axisA);
+
+  // Get the dispersion points for the semi_axis B
+  vector<WeightPoint> weights_semi_axisB;
+  semi_axisB.get_weights(weights_semi_axisB);
+
+  // Get the dispersion points for the semi_axis C
+  vector<WeightPoint> weights_semi_axisC;
+  semi_axisC.get_weights(weights_semi_axisC);
+
+  // Get angular averaging for theta
+  vector<WeightPoint> weights_theta;
+  axis_theta.get_weights(weights_theta);
+
+  // Get angular averaging for phi
+  vector<WeightPoint> weights_phi;
+  axis_phi.get_weights(weights_phi);
+
+  // Get angular averaging for psi
+  vector<WeightPoint> weights_psi;
+  axis_psi.get_weights(weights_psi);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double norm_vol = 0.0;
+  double vol = 0.0;
+  double pi = 4.0*atan(1.0);
+  // Loop over semi axis A weight points
+  for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
+    dp.semi_axisA = weights_semi_axisA[i].value;
+
+    // Loop over semi axis B weight points
+    for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
+      dp.semi_axisB = weights_semi_axisB[j].value;
+
+      // Loop over semi axis C weight points
+      for(int k=0; k < (int)weights_semi_axisC.size(); k++) {
+        dp.semi_axisC = weights_semi_axisC[k].value;
+
+        // Average over theta distribution
+        for(int l=0; l< (int)weights_theta.size(); l++) {
+          dp.axis_theta = weights_theta[l].value;
+
+          // Average over phi distribution
+          for(int m=0; m <(int)weights_phi.size(); m++) {
+            dp.axis_phi = weights_phi[m].value;
+            // Average over psi distribution
+            for(int n=0; n <(int)weights_psi.size(); n++) {
+              dp.axis_psi = weights_psi[n].value;
+              //Un-normalize  by volume
+              double _ptvalue = weights_semi_axisA[i].weight
+                  * weights_semi_axisB[j].weight
+                  * weights_semi_axisC[k].weight
+                  * weights_theta[l].weight
+                  * weights_phi[m].weight
+                  * weights_psi[n].weight
+                  * triaxial_ellipsoid_analytical_2DXY(&dp, qx, qy)
+              * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
+              if (weights_theta.size()>1) {
+                _ptvalue *= fabs(sin(weights_theta[k].value*pi/180.0));
+              }
+              sum += _ptvalue;
+              //Find average volume
+              vol += weights_semi_axisA[i].weight
+                  * weights_semi_axisB[j].weight
+                  * weights_semi_axisC[k].weight
+                  * weights_semi_axisA[i].value*weights_semi_axisB[j].value*weights_semi_axisC[k].value;
+              //Find norm for volume
+              norm_vol += weights_semi_axisA[i].weight
+                  * weights_semi_axisB[j].weight
+                  * weights_semi_axisC[k].weight;
+
+              norm += weights_semi_axisA[i].weight
+                  * weights_semi_axisB[j].weight
+                  * weights_semi_axisC[k].weight
+                  * weights_theta[l].weight
+                  * weights_phi[m].weight
+                  * weights_psi[n].weight;
+            }
+          }
+
+        }
+      }
+    }
+  }
+  // Averaging in theta needs an extra normalization
+  // factor to account for the sin(theta) term in the
+  // integration (see documentation).
+  if (weights_theta.size()>1) norm = norm / asin(1.0);
+
+  if (vol != 0.0 && norm_vol != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm_vol);}
+
+  return sum/norm + background();
 }
 
@@ -245 +371 @@
  */
 double TriaxialEllipsoidModel :: evaluate_rphi(double q, double phi) {
-        double qx = q*cos(phi);
-        double qy = q*sin(phi);
-        return (*this).operator()(qx, qy);
+  double qx = q*cos(phi);
+  double qy = q*sin(phi);
+  return (*this).operator()(qx, qy);
 }
 /**
@@ -254 +380 @@
  */
 double TriaxialEllipsoidModel :: calculate_ER() {
-        TriaxialEllipsoidParameters dp;
-
-        dp.semi_axisA   = semi_axisA();
-        dp.semi_axisB     = semi_axisB();
-        //polar axis C
-        dp.semi_axisC     = semi_axisC();
-
-        double rad_out = 0.0;
-        //Surface average radius at the equat. cross section.
-        double suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB);
-
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-
-        // Get the dispersion points for the semi_axis A
-        vector<WeightPoint> weights_semi_axisA;
-        semi_axisA.get_weights(weights_semi_axisA);
-
-        // Get the dispersion points for the semi_axis B
-        vector<WeightPoint> weights_semi_axisB;
-        semi_axisB.get_weights(weights_semi_axisB);
-
-        // Get the dispersion points for the semi_axis C
-        vector<WeightPoint> weights_semi_axisC;
-        semi_axisC.get_weights(weights_semi_axisC);
-
-        // Loop over semi axis A weight points
-        for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
-                dp.semi_axisA = weights_semi_axisA[i].value;
-
-                // Loop over semi axis B weight points
-                for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
-                        dp.semi_axisB = weights_semi_axisB[j].value;
-
-                        // Loop over semi axis C weight points
-                        for(int k=0; k < (int)weights_semi_axisC.size(); k++) {
-                                dp.semi_axisC = weights_semi_axisC[k].value;
-
-                                //Calculate surface averaged radius
-                                suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB);
-
-                                //Sum
-                                sum += weights_semi_axisA[i].weight
-                                        * weights_semi_axisB[j].weight
-                                        * weights_semi_axisC[k].weight * DiamEllip(dp.semi_axisC, suf_rad)/2.0;
-                                //Norm
-                                norm += weights_semi_axisA[i].weight* weights_semi_axisB[j].weight
-                                        * weights_semi_axisC[k].weight;
-                        }
-                }
-        }
-        if (norm != 0){
-                //return the averaged value
-                rad_out =  sum/norm;}
-        else{
-                //return normal value
-                rad_out = DiamEllip(dp.semi_axisC, suf_rad)/2.0;}
-
-        return rad_out;
-}
+  TriaxialEllipsoidParameters dp;
+
+  dp.semi_axisA   = semi_axisA();
+  dp.semi_axisB     = semi_axisB();
+  //polar axis C
+  dp.semi_axisC     = semi_axisC();
+
+  double rad_out = 0.0;
+  //Surface average radius at the equat. cross section.
+  double suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB);
+
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+
+  // Get the dispersion points for the semi_axis A
+  vector<WeightPoint> weights_semi_axisA;
+  semi_axisA.get_weights(weights_semi_axisA);
+
+  // Get the dispersion points for the semi_axis B
+  vector<WeightPoint> weights_semi_axisB;
+  semi_axisB.get_weights(weights_semi_axisB);
+
+  // Get the dispersion points for the semi_axis C
+  vector<WeightPoint> weights_semi_axisC;
+  semi_axisC.get_weights(weights_semi_axisC);
+
+  // Loop over semi axis A weight points
+  for(int i=0; i< (int)weights_semi_axisA.size(); i++) {
+    dp.semi_axisA = weights_semi_axisA[i].value;
+
+    // Loop over semi axis B weight points
+    for(int j=0; j< (int)weights_semi_axisB.size(); j++) {
+      dp.semi_axisB = weights_semi_axisB[j].value;
+
+      // Loop over semi axis C weight points
+      for(int k=0; k < (int)weights_semi_axisC.size(); k++) {
+        dp.semi_axisC = weights_semi_axisC[k].value;
+
+        //Calculate surface averaged radius
+        suf_rad = sqrt(dp.semi_axisA * dp.semi_axisB);
+
+        //Sum
+        sum += weights_semi_axisA[i].weight
+            * weights_semi_axisB[j].weight
+            * weights_semi_axisC[k].weight * DiamEllip(dp.semi_axisC, suf_rad)/2.0;
+        //Norm
+        norm += weights_semi_axisA[i].weight* weights_semi_axisB[j].weight
+            * weights_semi_axisC[k].weight;
+      }
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    rad_out = DiamEllip(dp.semi_axisC, suf_rad)/2.0;}
+
+  return rad_out;
+}

sansmodels/src/c_models/vesicle.cpp

@@ -21 +21 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "vesicle.h"
 
 extern "C" {
-        #include "libSphere.h"
-        #include "vesicle.h"
+#include "libSphere.h"
 }
 
+typedef struct {
+  double scale;
+  double radius;
+  double thickness;
+  double core_sld;
+  double shell_sld;
+  double background;
+} VesicleParameters;
+
 VesicleModel :: VesicleModel() {
-        scale      = Parameter(1.0);
-        radius     = Parameter(100.0, true);
-        radius.set_min(0.0);
-        thickness  = Parameter(30.0, true);
-        thickness.set_min(0.0);
-        core_sld   = Parameter(6.36e-6);
-        shell_sld   = Parameter(5.0e-7);
-        background = Parameter(0.0);
+  scale      = Parameter(1.0);
+  radius     = Parameter(100.0, true);
+  radius.set_min(0.0);
+  thickness  = Parameter(30.0, true);
+  thickness.set_min(0.0);
+  core_sld   = Parameter(6.36e-6);
+  shell_sld   = Parameter(5.0e-7);
+  background = Parameter(0.0);
 }
 
@@ -49 +57 @@
  */
 double VesicleModel :: operator()(double q) {
-        double dp[6];
+  double dp[6];
 
-        // Fill parameter array for IGOR library
-        // Add the background after averaging
-        dp[0] = scale();
-        dp[1] = radius();
-        dp[2] = thickness();
-        dp[3] = core_sld();
-        dp[4] = shell_sld();
-        dp[5] = 0.0;
+  // Fill parameter array for IGOR library
+  // Add the background after averaging
+  dp[0] = scale();
+  dp[1] = radius();
+  dp[2] = thickness();
+  dp[3] = core_sld();
+  dp[4] = shell_sld();
+  dp[5] = 0.0;
 
 
-        // Get the dispersion points for the core radius
-        vector<WeightPoint> weights_radius;
-        radius.get_weights(weights_radius);
-        // Get the dispersion points for the thickness
-        vector<WeightPoint> weights_thickness;
-        thickness.get_weights(weights_thickness);
+  // Get the dispersion points for the core radius
+  vector<WeightPoint> weights_radius;
+  radius.get_weights(weights_radius);
+  // Get the dispersion points for the thickness
+  vector<WeightPoint> weights_thickness;
+  thickness.get_weights(weights_thickness);
 
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
-        double vol = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
+  double vol = 0.0;
 
-        // Loop over radius weight points
-        for(int i=0; i< (int)weights_radius.size(); i++) {
-                dp[1] = weights_radius[i].value;
-                for(int j=0; j< (int)weights_thickness.size(); j++) {
-                        dp[2] = weights_thickness[j].value;
-                        sum += weights_radius[i].weight
-                                * weights_thickness[j].weight * VesicleForm(dp, q)
-                                *(pow(weights_radius[i].value+weights_thickness[j].value,3)-pow(weights_radius[i].value,3));
-                        //Find average volume
-                        vol += weights_radius[i].weight * weights_thickness[j].weight
-                                *(pow(weights_radius[i].value+weights_thickness[j].value,3)-pow(weights_radius[i].value,3));
-                        norm += weights_radius[i].weight * weights_thickness[j].weight;
-                }
-        }
-        if (vol != 0.0 && norm != 0.0) {
-                //Re-normalize by avg volume
-                sum = sum/(vol/norm);}
+  // Loop over radius weight points
+  for(int i=0; i< (int)weights_radius.size(); i++) {
+    dp[1] = weights_radius[i].value;
+    for(int j=0; j< (int)weights_thickness.size(); j++) {
+      dp[2] = weights_thickness[j].value;
+      sum += weights_radius[i].weight
+          * weights_thickness[j].weight * VesicleForm(dp, q)
+      *(pow(weights_radius[i].value+weights_thickness[j].value,3)-pow(weights_radius[i].value,3));
+      //Find average volume
+      vol += weights_radius[i].weight * weights_thickness[j].weight
+          *(pow(weights_radius[i].value+weights_thickness[j].value,3)-pow(weights_radius[i].value,3));
+      norm += weights_radius[i].weight * weights_thickness[j].weight;
+    }
+  }
+  if (vol != 0.0 && norm != 0.0) {
+    //Re-normalize by avg volume
+    sum = sum/(vol/norm);}
 
-        return sum/norm + background();
+  return sum/norm + background();
 }
 
@@ -101 +109 @@
  */
 double VesicleModel :: operator()(double qx, double qy) {
-        double q = sqrt(qx*qx + qy*qy);
-        return (*this).operator()(q);
+  double q = sqrt(qx*qx + qy*qy);
+  return (*this).operator()(q);
 }
 
@@ -113 +121 @@
  */
 double VesicleModel :: evaluate_rphi(double q, double phi) {
-        return (*this).operator()(q);
+  return (*this).operator()(q);
 }
 /**
@@ -120 +128 @@
  */
 double VesicleModel :: calculate_ER() {
-        VesicleParameters dp;
+  VesicleParameters dp;
 
-        dp.radius     = radius();
-        dp.thickness  = thickness();
+  dp.radius     = radius();
+  dp.thickness  = thickness();
 
-        double rad_out = 0.0;
+  double rad_out = 0.0;
 
-        // Perform the computation, with all weight points
-        double sum = 0.0;
-        double norm = 0.0;
+  // Perform the computation, with all weight points
+  double sum = 0.0;
+  double norm = 0.0;
 
 
-        // Get the dispersion points for the major shell
-        vector<WeightPoint> weights_thickness;
-        thickness.get_weights(weights_thickness);
+  // Get the dispersion points for the major shell
+  vector<WeightPoint> weights_thickness;
+  thickness.get_weights(weights_thickness);
 
-        // Get the dispersion points for the minor shell
-        vector<WeightPoint> weights_radius ;
-        radius.get_weights(weights_radius);
+  // Get the dispersion points for the minor shell
+  vector<WeightPoint> weights_radius ;
+  radius.get_weights(weights_radius);
 
-        // Loop over major shell weight points
-        for(int j=0; j< (int)weights_thickness.size(); j++) {
-                dp.thickness = weights_thickness[j].value;
-                for(int k=0; k< (int)weights_radius.size(); k++) {
-                        dp.radius = weights_radius[k].value;
-                        sum += weights_thickness[j].weight
-                                * weights_radius[k].weight*(dp.radius+dp.thickness);
-                        norm += weights_thickness[j].weight* weights_radius[k].weight;
-                }
-        }
-        if (norm != 0){
-                //return the averaged value
-                rad_out =  sum/norm;}
-        else{
-                //return normal value
-                rad_out = (dp.radius+dp.thickness);}
+  // Loop over major shell weight points
+  for(int j=0; j< (int)weights_thickness.size(); j++) {
+    dp.thickness = weights_thickness[j].value;
+    for(int k=0; k< (int)weights_radius.size(); k++) {
+      dp.radius = weights_radius[k].value;
+      sum += weights_thickness[j].weight
+          * weights_radius[k].weight*(dp.radius+dp.thickness);
+      norm += weights_thickness[j].weight* weights_radius[k].weight;
+    }
+  }
+  if (norm != 0){
+    //return the averaged value
+    rad_out =  sum/norm;}
+  else{
+    //return normal value
+    rad_out = (dp.radius+dp.thickness);}
 
-        return rad_out;
+  return rad_out;
 }