Changeset 8343e18 in sasview

Timestamp:

Jan 4, 2012 6:16:15 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:

df88829

Parents:

0c2389e

Message:

refactored parallelepiped

Location:

sansmodels/src

Files:

• c_extensions/parallelepiped.c (deleted)
• c_extensions/parallelepiped.h (modified) (2 diffs)
• c_models/models.hh (modified) (1 diff)
• c_models/parallelepiped.cpp (modified) (2 diffs)
• python_wrapper/CParallelepipedModel.cpp (modified) (2 diffs)

sansmodels/src/c_extensions/parallelepiped.h

@@ -1 @@
-/*
-        TODO: Add 2D model
-*/
-
 #if !defined(parallelepiped_h)
 #define parallelepiped_h
+#include "parameters.hh"
+
 /** Structure definition for Parallelepiped parameters
  * [PYTHONCLASS] = ParallelepipedModel
@@ -23 +21 @@
 
  **/
-typedef struct {
-    /// Scale factor
-    //  [DEFAULT]=scale=1.0
-    double scale;
-    ///  Length of short edge of the parallelepiped [A]
-    //  [DEFAULT]=short_a=35 [A]
-    double short_a;
-        /// Length of short edge edge of the parallelepiped [A]
-    //  [DEFAULT]=short_b=75 [A]
-    double short_b;
-        /// Length of long edge of the parallelepiped [A]
-    //  [DEFAULT]=long_c=400 [A]
-    double long_c;
-    /// SLD_Pipe [1/A^(2)]
-    //  [DEFAULT]=sldPipe=6.3e-6 [1/A^(2)]
-    double sldPipe;
-    /// sldSolv [1/A^(2)]
-    //  [DEFAULT]=sldSolv=1.0e-6 [1/A^(2)]
-    double sldSolv;
-        /// Incoherent Background [1/cm]
-        //  [DEFAULT]=background=0.0 [1/cm]
-        double background;
-    /// Orientation of the parallelepiped axis w/respect incoming beam [deg]
-    //  [DEFAULT]=parallel_theta=0.0 [deg]
-    double parallel_theta;
-    /// Orientation of the longitudinal axis of the parallelepiped in the plane of the detector [deg]
-    //  [DEFAULT]=parallel_phi=0.0 [deg]
-    double parallel_phi;
-    /// Orientation of the cross-sectional minor axis of the parallelepiped in the plane of the detector [deg]
-    //  [DEFAULT]=parallel_psi=0.0 [deg]
-    double parallel_psi;
 
+class ParallelepipedModel{
+public:
+  // Model parameters
+  /// Scale factor
+  //  [DEFAULT]=scale=1.0
+  Parameter scale;
+  ///  Length of short edge of the parallelepiped [A]
+  //  [DEFAULT]=short_a=35 [A]
+  Parameter short_a;
+  /// Length of short edge edge of the parallelepiped [A]
+  //  [DEFAULT]=short_b=75 [A]
+  Parameter short_b;
+  /// Length of long edge of the parallelepiped [A]
+  //  [DEFAULT]=long_c=400 [A]
+  Parameter long_c;
+  /// SLD_Pipe [1/A^(2)]
+  //  [DEFAULT]=sldPipe=6.3e-6 [1/A^(2)]
+  Parameter sldPipe;
+  /// sldSolv [1/A^(2)]
+  //  [DEFAULT]=sldSolv=1.0e-6 [1/A^(2)]
+  Parameter sldSolv;
+  /// Incoherent Background [1/cm]
+  //  [DEFAULT]=background=0.0 [1/cm]
+  Parameter background;
+  /// Orientation of the parallelepiped axis w/respect incoming beam [deg]
+  //  [DEFAULT]=parallel_theta=0.0 [deg]
+  Parameter parallel_theta;
+  /// Orientation of the longitudinal axis of the parallelepiped in the plane of the detector [deg]
+  //  [DEFAULT]=parallel_phi=0.0 [deg]
+  Parameter parallel_phi;
+  /// Orientation of the cross-sectional minor axis of the parallelepiped in the plane of the detector [deg]
+  //  [DEFAULT]=parallel_psi=0.0 [deg]
+  Parameter parallel_psi;
 
-} ParallelepipedParameters;
+  // Constructor
+  ParallelepipedModel();
 
-
-
-/// 1D scattering function
-double parallelepiped_analytical_1D(ParallelepipedParameters *pars, double q);
-
-/// 2D scattering function
-double parallelepiped_analytical_2D(ParallelepipedParameters *pars, double q, double phi);
-double parallelepiped_analytical_2DXY(ParallelepipedParameters *pars, double qx, double qy);
-double parallelepiped_analytical_2D_scaled(ParallelepipedParameters *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/models.hh

@@ -26 +26 @@
 
 using namespace std;
-
-
-
-class ParallelepipedModel{
-public:
-        // Model parameters
-        Parameter scale;
-        Parameter short_a;
-        Parameter short_b;
-        Parameter long_c;
-        Parameter sldPipe;
-        Parameter sldSolv;
-        Parameter background;
-        Parameter parallel_theta;
-        Parameter parallel_phi;
-        Parameter parallel_psi;
-
-        // Constructor
-        ParallelepipedModel();
-
-        // 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);
-};
-
 
 class CSParallelepipedModel{

sansmodels/src/c_models/parallelepiped.cpp

@@ -18 +18 @@
  *   sansmodels/src/libigor
  *
- *      TODO: refactor so that we pull in the old sansmodels.c_extensions
  *      TODO: add 2D function
  */
@@ -31 +30 @@
         #include "libCylinder.h"
         #include "libStructureFactor.h"
-        #include "parallelepiped.h"
+}
+#include "parallelepiped.h"
+
+// Convenience parameter structure
+typedef struct {
+    double scale;
+    double short_a;
+    double short_b;
+    double long_c;
+    double sldPipe;
+    double sldSolv;
+    double background;
+    double parallel_theta;
+    double parallel_phi;
+    double parallel_psi;
+} ParallelepipedParameters;
+
+
+static double pkernel(double a, double b,double c, double ala, double alb, double alc){
+    // mu passed in is really mu*sqrt(1-sig^2)
+    double argA,argB,argC,tmp1,tmp2,tmp3;     //local variables
+
+    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
+    argA = a*ala/2.0;
+    argB = b*alb/2.0;
+    argC = c*alc/2.0;
+    if(argA==0.0) {
+    tmp1 = 1.0;
+  } else {
+    tmp1 = sin(argA)*sin(argA)/argA/argA;
+    }
+    if (argB==0.0) {
+    tmp2 = 1.0;
+  } else {
+    tmp2 = sin(argB)*sin(argB)/argB/argB;
+    }
+
+    if (argC==0.0) {
+    tmp3 = 1.0;
+  } else {
+    tmp3 = sin(argC)*sin(argC)/argC/argC;
+    }
+
+    return (tmp1*tmp2*tmp3);
+
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the parallelepiped
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double parallelepiped_analytical_2D_scaled(ParallelepipedParameters *pars, double q, double q_x, double q_y) {
+  double cparallel_x, cparallel_y, cparallel_z, bparallel_x, bparallel_y, parallel_x, parallel_y;
+  double q_z;
+  double alpha, vol, cos_val_c, cos_val_b, cos_val_a, edgeA, edgeB, edgeC;
+
+  double answer;
+  double pi = 4.0*atan(1.0);
+
+  //convert angle degree to radian
+  double theta = pars->parallel_theta * pi/180.0;
+  double phi = pars->parallel_phi * pi/180.0;
+  double psi = pars->parallel_psi * pi/180.0;
+
+  edgeA = pars->short_a;
+  edgeB = pars->short_b;
+  edgeC = pars->long_c;
+
+
+    // parallelepiped c axis orientation
+    cparallel_x = sin(theta) * cos(phi);
+    cparallel_y = sin(theta) * sin(phi);
+    cparallel_z = cos(theta);
+
+    // q vector
+    q_z = 0.0;
+
+    // Compute the angle btw vector q and the
+    // axis of the parallelepiped
+    cos_val_c = cparallel_x*q_x + cparallel_y*q_y + cparallel_z*q_z;
+    alpha = acos(cos_val_c);
+
+    // parallelepiped a axis orientation
+    parallel_x = sin(psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)*sin(pars->parallel_psi);
+    parallel_y = cos(psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)*cos(pars->parallel_psi);
+
+    cos_val_a = parallel_x*q_x + parallel_y*q_y;
+
+
+
+    // parallelepiped b axis orientation
+    bparallel_x = sqrt(1.0-sin(theta)*cos(phi))*cos(psi);//cos(pars->parallel_theta) * cos(pars->parallel_phi)* cos(pars->parallel_psi);
+    bparallel_y = sqrt(1.0-sin(theta)*cos(phi))*sin(psi);//cos(pars->parallel_theta) * sin(pars->parallel_phi)* sin(pars->parallel_psi);
+    // axis of the parallelepiped
+    cos_val_b = sin(acos(cos_val_a)) ;
+
+
+
+    // The following test should always pass
+    if (fabs(cos_val_c)>1.0) {
+      printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+      return 0;
+    }
+
+  // Call the IGOR library function to get the kernel
+  answer = pkernel( q*edgeA, q*edgeB, q*edgeC, sin(alpha)*cos_val_a,sin(alpha)*cos_val_b,cos_val_c);
+
+  // Multiply by contrast^2
+  answer *= (pars->sldPipe - pars->sldSolv) * (pars->sldPipe - pars->sldSolv);
+
+  //normalize by cylinder volume
+  //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vparallel
+    vol = edgeA* edgeB * edgeC;
+  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 parallelepiped
+ * @param q: q-value
+ * @return: function value
+ */
+static double parallelepiped_analytical_2DXY(ParallelepipedParameters *pars, double qx, double qy) {
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+    return parallelepiped_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }
 

sansmodels/src/python_wrapper/CParallelepipedModel.cpp

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