Changeset 318b5bbb in sasview for sansmodels/src/c_models/triaxialellipsoid.cpp

Timestamp:

Dec 18, 2012 10:55:24 AM (12 years ago)

Author:

Jae Cho <jhjcho@…>

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:

6550b64

Parents:

0203ade

Message:

Added polarization and magnetic stuffs

File:

• sansmodels/src/c_models/triaxialellipsoid.cpp (modified) (8 diffs)

sansmodels/src/c_models/triaxialellipsoid.cpp

@@ -46 +46 @@
 } TriaxialEllipsoidParameters;
 
-static double triaxial_ellipsoid_kernel(TriaxialEllipsoidParameters *pars, double q, double alpha, double nu) {
+static double triaxial_ellipsoid_kernel(TriaxialEllipsoidParameters *pars, double q, double cos_val, double cos_nu, double cos_mu) {
   double t,a,b,c;
   double kernel;
@@ -54 +54 @@
   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));
+  t = q * sqrt(a*a*cos_nu*cos_nu+b*b*cos_mu*cos_mu+c*c*cos_val*cos_val);
   if (t==0.0){
     kernel  = 1.0;
@@ -73 +73 @@
  */
 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 cyl_x, cyl_y, ella_x, ella_y, ellb_x, ellb_y;
+  //double q_z;
+  double cos_nu, cos_mu;
+  double vol, cos_val;
   double answer;
   double pi = 4.0*atan(1.0);
@@ -86 +86 @@
 
   // Cylinder orientation
-  cyl_x = sin(theta) * cos(phi);
-  cyl_y = sin(theta) * sin(phi);
-  cyl_z = cos(theta);
+  cyl_x = cos(theta) * cos(phi);
+  cyl_y = sin(theta);
+  //cyl_z = -cos(theta) * sin(phi);
 
   // q vector
-  q_z = 0.0;
+  //q_z = 0.0;
 
   //dx = 1.0;
@@ -97 +97 @@
   // 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;
+  cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
 
   // The following test should always pass
@@ -107 +107 @@
   // Note: cos(alpha) = 0 and 1 will get an
   // undefined value from CylKernel
-  alpha = acos( cos_val );
+  //alpha = acos( cos_val );
 
   //ellipse orientation:
@@ -116 +116 @@
   // the wave vector q.
 
-  //x- y- component on the detector plane.
-  ell_x =  cos(psi);
-  ell_y =  sin(psi);
-
+  //x- y- component of a-axis on the detector plane.
+  ella_x =  -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
+  ella_y =  sin(psi)*cos(theta);
+  
+  //x- y- component of b-axis on the detector plane.
+  ellb_x =  -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
+  ellb_y =  cos(theta)*cos(psi);
+  
   // calculate the axis of the ellipse wrt q-coord.
-  cos_nu = ell_x*q_x + ell_y*q_y;
-  nu = acos(cos_nu);
-
+  cos_nu = ella_x*q_x + ella_y*q_y;
+  cos_mu = ellb_x*q_x + ellb_y*q_y;
+  
+  // The following test should always pass
+  if (fabs(cos_val)>1.0) {
+    //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+    cos_val = 1.0;
+  }
+   if (fabs(cos_nu)>1.0) {
+    //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+    cos_nu = 1.0;
+  }
+   if (fabs(cos_mu)>1.0) {
+    //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+    cos_mu = 1.0;
+  } 
   // Call the IGOR library function to get the kernel
-  answer = triaxial_ellipsoid_kernel(pars, q, alpha, nu);
+  answer = triaxial_ellipsoid_kernel(pars, q, cos_val, cos_nu, cos_mu);
 
   // Multiply by contrast^2
@@ -325 +342 @@
               * 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));
+                _ptvalue *= fabs(cos(weights_theta[k].value*pi/180.0));
               }
               sum += _ptvalue;