Changeset 0ba3b08 in sasview for sansmodels/src/c_models/bcc.cpp

Timestamp:

Jan 5, 2012 12:16:29 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:

011e0e4

Parents:

bbbed8c

Message:

refactored bunch of models

File:

• sansmodels/src/c_models/bcc.cpp (modified) (1 diff)

sansmodels/src/c_models/bcc.cpp

@@ -21 +21 @@
 
 #include <math.h>
-#include "models.hh"
 #include "parameters.hh"
 #include <stdio.h>
 using namespace std;
+#include "bcc.h"
 
 extern "C" {
         #include "libSphere.h"
-        #include "bcc.h"
+}
+
+// Convenience structure
+typedef struct {
+  double scale;
+  double dnn;
+  double d_factor;
+  double radius;
+  double sldSph;
+  double sldSolv;
+  double background;
+  double theta;
+  double phi;
+  double psi;
+
+} BCParameters;
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the BCCCrystalModel
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+static double bc_analytical_2D_scaled(BCParameters *pars, double q, double q_x, double q_y) {
+  double b3_x, b3_y, b3_z, b1_x, b1_y;
+  double q_z;
+  double alpha, cos_val_b3, cos_val_b2, cos_val_b1;
+  double a1_dot_q, a2_dot_q,a3_dot_q;
+  double answer;
+  double Pi = 4.0*atan(1.0);
+  double aa, Da, qDa_2, latticeScale, Zq, Fkq, Fkq_2;
+
+  //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;
+  double psi = pars->psi * pi/180.0;
+
+  double dp[5];
+  dp[0] = 1.0;
+  dp[1] = pars->radius;
+  dp[2] = pars->sldSph;
+  dp[3] = pars->sldSolv;
+  dp[4] = 0.0;
+
+  aa = pars->dnn;
+  Da = pars->d_factor*aa;
+  qDa_2 = pow(q*Da,2.0);
+
+  //the occupied volume of the lattice
+  latticeScale = 2.0*(4.0/3.0)*Pi*(dp[1]*dp[1]*dp[1])/pow(aa/sqrt(3.0/4.0),3.0);
+  // q vector
+  q_z = 0.0; // for SANS; assuming qz is negligible
+  /// Angles here are respect to detector coordinate
+  ///  instead of against q coordinate(PRB 36(46), 3(6), 1754(3854))
+    // b3 axis orientation
+    b3_x = sin(theta) * cos(phi);//negative sign here???
+    b3_y = sin(theta) * sin(phi);
+    b3_z = cos(theta);
+    cos_val_b3 =  b3_x*q_x + b3_y*q_y + b3_z*q_z;
+
+    alpha = acos(cos_val_b3);
+    // b1 axis orientation
+    b1_x = sin(psi);
+    b1_y = cos(psi);
+  cos_val_b1 = (b1_x*q_x + b1_y*q_y);
+    // b2 axis orientation
+  cos_val_b2 = sin(acos(cos_val_b1));
+  // alpha corrections
+  cos_val_b2 *= sin(alpha);
+  cos_val_b1 *= sin(alpha);
+
+    // Compute the angle btw vector q and the a3 axis
+    a3_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b1-cos_val_b3);
+
+    // a1 axis
+    a1_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b2-cos_val_b1);
+
+    // a2 axis
+    a2_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b1-cos_val_b2);
+
+    // The following test should always pass
+    if (fabs(cos_val_b3)>1.0) {
+      printf("bcc_ana_2D: Unexpected error: cos(alpha)>1\n");
+      return 0;
+    }
+    // Get Fkq and Fkq_2
+    Fkq = exp(-0.5*pow(Da/aa,2.0)*(a1_dot_q*a1_dot_q+a2_dot_q*a2_dot_q+a3_dot_q*a3_dot_q));
+    Fkq_2 = Fkq*Fkq;
+    // Call Zq=Z1*Z2*Z3
+    Zq = (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a1_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a2_dot_q)+Fkq_2);
+    Zq *= (1.0-Fkq_2)/(1.0-2.0*Fkq*cos(a3_dot_q)+Fkq_2);
+
+  // Use SphereForm directly from libigor
+  answer = SphereForm(dp,q)*Zq;
+
+  //consider scales
+  answer *= latticeScale * pars->scale;
+
+  // This FIXES a singualrity the kernel in libigor.
+  if ( answer == INFINITY || answer == NAN){
+    answer = 0.0;
+  }
+
+  // add background
+  answer += pars->background;
+
+  return answer;
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the BCC_ParaCrystal
+ * @param q: q-value
+ * @return: function value
+ */
+static double bc_analytical_2DXY(BCParameters *pars, double qx, double qy){
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return bc_analytical_2D_scaled(pars, q, qx/q, qy/q);
 }