Changeset 0c2389e in sasview for sansmodels

Timestamp:

Jan 4, 2012 6:00:52 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:

8343e18

Parents:

dbddbf5

Message:

refactored capped cylinder

Location:

sansmodels/src

Files:

• Make.mm (deleted)
• c_extensions/capcyl.c (deleted)
• c_extensions/capcyl.h (modified) (2 diffs)
• c_models/capcyl.cpp (modified) (3 diffs)
• c_models/cylinder.cpp (modified) (1 diff)
• c_models/models.hh (modified) (1 diff)
• python_wrapper/CCappedCylinderModel.cpp (modified) (2 diffs)

sansmodels/src/c_extensions/capcyl.h

@@ -1 +1 @@
 #if !defined(capcyl_h)
 #define capcyl_h
+#include "parameters.hh"
 
 /**
@@ -27 +28 @@
  //[ORIENTATION_PARAMS]= <text> phi; theta; phi.width; theta.width</text>
 
-typedef struct {
-    /// Scale factor
-    //  [DEFAULT]=scale= 1.0
-    double scale;
+class CappedCylinderModel{
+public:
+  // Model parameters
+  /// Scale factor
+  //  [DEFAULT]=scale= 1.0
+  Parameter scale;
 
-    /// rad_cyl [A]
-    //  [DEFAULT]=rad_cyl=20.0 [A]
-    double rad_cyl;
+  /// rad_cyl [A]
+  //  [DEFAULT]=rad_cyl=20.0 [A]
+  Parameter rad_cyl;
 
-    /// length of the cylinder
-    //  [DEFAULT]=len_cyl=400.0 [A]
-    double len_cyl;
+  /// length of the cylinder
+  //  [DEFAULT]=len_cyl=400.0 [A]
+  Parameter len_cyl;
 
-    /// Radius of sphere [A]
-    //  [DEFAULT]=rad_cap=40.0 [A]
-    double rad_cap;
+  /// Radius of sphere [A]
+  //  [DEFAULT]=rad_cap=40.0 [A]
+  Parameter rad_cap;
 
-    /// sld_capcyl [1/A^(2)]
-    //  [DEFAULT]=sld_capcyl= 1.0e-6 [1/A^(2)]
-    double sld_capcyl;
+  /// sld_capcyl [1/A^(2)]
+  //  [DEFAULT]=sld_capcyl= 1.0e-6 [1/A^(2)]
+  Parameter sld_capcyl;
 
-    /// sld_solv [1/A^(2)]
-    //  [DEFAULT]=sld_solv= 6.3e-6 [1/A^(2)]
-    double sld_solv;
+  /// sld_solv [1/A^(2)]
+  //  [DEFAULT]=sld_solv= 6.3e-6 [1/A^(2)]
+  Parameter sld_solv;
 
-        /// Incoherent Background [1/cm]
-        //  [DEFAULT]=background=0.0 [1/cm]
-        double background;
+  /// Incoherent Background [1/cm]
+  //  [DEFAULT]=background=0.0 [1/cm]
+  Parameter background;
 
-    /// Angle of the main axis against z-axis in detector plane [deg]
-    //  [DEFAULT]=theta=0.0 [deg]
-    double theta;
-    /// Azimuthal angle around z-axis in detector plane [deg]
-    //  [DEFAULT]=phi=0.0 [deg]
-    double phi;
+  /// Angle of the main axis against z-axis in detector plane [deg]
+  //  [DEFAULT]=theta=0.0 [deg]
+  Parameter theta;
 
-} CapCylParameters;
+  /// Azimuthal angle around z-axis in detector plane [deg]
+  //  [DEFAULT]=phi=0.0 [deg]
+  Parameter phi;
 
+  // Constructor
+  CappedCylinderModel();
 
-
-/// 1D scattering function
-double capcyl_analytical_1D(CapCylParameters *pars, double q);
-
-/// 2D scattering function
-double capcyl_analytical_2D(CapCylParameters *pars, double q, double phi);
-double capcyl_analytical_2DXY(CapCylParameters *pars, double qx, double qy);
-double capcyl_analytical_2D_scaled(CapCylParameters *pars, double q, double q_x, double q_y);
+  // 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);
+};
 
 #endif

sansmodels/src/c_models/capcyl.cpp

@@ -21 +21 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
@@ -27 +26 @@
 
 extern "C" {
+  #include "GaussWeights.h"
         #include "libCylinder.h"
-        #include "capcyl.h"
-}
+}
+#include "capcyl.h"
+
+// Convenience parameter structure
+typedef struct {
+    double scale;
+    double rad_cyl;
+    double len_cyl;
+    double rad_cap;
+    double sld_capcyl;
+    double sld_solv;
+    double background;
+    double theta;
+    double phi;
+} CapCylParameters;
 
 CappedCylinderModel :: CappedCylinderModel() {
@@ -44 +57 @@
         theta  = Parameter(0.0, true);
         phi    = Parameter(0.0, true);
+}
+
+static double capcyl2d_kernel(double dp[], double q, double alpha) {
+  int j;
+  double Pi;
+  double scale,contr,bkg,sldc,slds;
+  double len,rad,hDist,endRad;
+  int nordj=76;
+  double zi=alpha,yyy,answer;     //running tally of integration
+  double summj,vaj,vbj,zij;     //for the inner integration
+  double arg1,arg2,inner,be;
+
+
+  scale = dp[0];
+  rad = dp[1];
+  len = dp[2];
+  endRad = dp[3];
+  sldc = dp[4];
+  slds = dp[5];
+  bkg = dp[6];
+
+  hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad));   //by definition for this model
+
+  contr = sldc-slds;
+
+  Pi = 4.0*atan(1.0);
+  vaj = -1.0*hDist/endRad;
+  vbj = 1.0;    //endpoints of inner integral
+
+  summj=0.0;
+
+  for(j=0;j<nordj;j++) {
+    //20 gauss points for the inner integral
+    zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;    //the "t" dummy
+    yyy = Gauss76Wt[j] * ConvLens_kernel(dp,q,zij,zi);    //uses the same Kernel as the Dumbbell, here L>0
+    summj += yyy;
+  }
+  //now calculate the value of the inner integral
+  inner = (vbj-vaj)/2.0*summj;
+  inner *= 4.0*Pi*endRad*endRad*endRad;
+
+  //now calculate outer integrand
+  arg1 = q*len/2.0*cos(zi);
+  arg2 = q*rad*sin(zi);
+  yyy = inner;
+
+  if(arg2 == 0) {
+    be = 0.5;
+  } else {
+    be = NR_BessJ1(arg2)/arg2;
+  }
+
+  if(arg1 == 0.0) {   //limiting value of sinc(0) is 1; sinc is not defined in math.h
+    yyy += Pi*rad*rad*len*2.0*be;
+  } else {
+    yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
+  }
+  yyy *= yyy;   //sin(zi);
+  answer = yyy;
+
+
+  answer /= Pi*rad*rad*len + 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*hDist*hDist/3.0);   //divide by volume
+  answer *= 1.0e8;    //convert to cm^-1
+  answer *= contr*contr;
+  answer *= scale;
+  answer += bkg;
+
+  return answer;
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the BarBell
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double capcyl_analytical_2D_scaled(CapCylParameters *pars, double q, double q_x, double q_y) {
+  double cyl_x, cyl_y, cyl_z;
+  double q_z;
+  double alpha, cos_val;
+  double answer;
+  double dp[7];
+  //convert angle degree to radian
+  double pi = 4.0*atan(1.0);
+  double theta = pars->theta * pi/180.0;
+  double phi = pars->phi * pi/180.0;
+
+  dp[0] = pars->scale;
+  dp[1] = pars->rad_cyl;
+  dp[2] = pars->len_cyl;
+  dp[3] = pars->rad_cap;
+  dp[4] = pars->sld_capcyl;
+  dp[5] = pars->sld_solv;
+  dp[6] = pars->background;
+
+
+  //double Pi = 4.0*atan(1.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 );
+
+  // Call the IGOR library function to get the kernel
+  answer = capcyl2d_kernel(dp, q, alpha)/sin(alpha);
+
+
+  return answer;
+
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the BarBell
+ * @param q: q-value
+ * @return: function value
+ */
+static double capcyl_analytical_2DXY(CapCylParameters *pars, double qx, double qy){
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return capcyl_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 

sansmodels/src/c_models/cylinder.cpp

@@ -28 +28 @@
         #include "libCylinder.h"
         #include "libStructureFactor.h"
-        #include "cylinder.h"
-}
+}
+#include "cylinder.h"
 
 // Convenience parameter structure

sansmodels/src/c_models/models.hh

@@ -27 +27 @@
 using namespace std;
 
-
-class CappedCylinderModel{
-public:
-        // Model parameters
-        Parameter scale;
-        Parameter rad_cyl;
-        Parameter len_cyl;
-        Parameter rad_cap;
-        Parameter sld_capcyl;
-        Parameter sld_solv;
-        Parameter background;
-        Parameter theta;
-        Parameter phi;
-
-        // Constructor
-        CappedCylinderModel();
-
-        // 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);
-};
 
 

sansmodels/src/python_wrapper/CCappedCylinderModel.cpp

@@ -18 +18 @@
  *
  * WARNING: THIS FILE WAS GENERATED BY WRAPPERGENERATOR.PY
- *          DO NOT MODIFY THIS FILE, MODIFY capcyl.h
+ *          DO NOT MODIFY THIS FILE, MODIFY ../c_extensions/capcyl.h
  *          AND RE-RUN THE GENERATOR SCRIPT
  *
@@ -33 +33 @@
 #include <math.h>
 #include <time.h>
+
+}
+
 #include "capcyl.h"
-}
-
-#include "models.hh"
 #include "dispersion_visitor.hh"