Changeset 318b5bbb in sasview for sansmodels/src/c_models

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

Location:

sansmodels/src/c_models

Files:

• TwoYukawa.cpp (added)
• barbell.cpp (modified) (2 diffs)
• bcc.cpp (modified) (4 diffs)
• capcyl.cpp (modified) (3 diffs)
• corefourshell.cpp (modified) (6 diffs)
• coreshellbicelle.cpp (modified) (3 diffs)
• coreshellcylinder.cpp (modified) (3 diffs)
• coreshellsphere.cpp (modified) (7 diffs)
• csparallelepiped.cpp (modified) (4 diffs)
• cylinder.cpp (modified) (7 diffs)
• ellipsoid.cpp (modified) (3 diffs)
• ellipticalcylinder.cpp (modified) (6 diffs)
• fcc.cpp (modified) (4 diffs)
• hollowcylinder.cpp (modified) (3 diffs)
• libfunc.c (modified) (4 diffs)
• parallelepiped.cpp (modified) (7 diffs)
• sc.cpp (modified) (3 diffs)
• sphere.cpp (modified) (5 diffs)
• spheroid.cpp (modified) (4 diffs)
• stackeddisks.cpp (modified) (4 diffs)
• triaxialellipsoid.cpp (modified) (8 diffs)

sansmodels/src/c_models/barbell.cpp

@@ -258 +258 @@
 
                                           // Cylinder orientation
-                                    const double cyl_x = sin(_theta * pi/180.0) * cos(_phi * pi/180.0);
-                                    const double cyl_y = sin(_theta * pi/180.0) * sin(_phi * pi/180.0);
+                                    const double cyl_x = cos(_theta * pi/180.0) * cos(_phi * pi/180.0);
+                                    const double cyl_y = sin(_theta * pi/180.0);
 
                                     // Compute the angle btw vector q and the
@@ -291 +291 @@
                                                 }
                                                 if (weights_theta.size()>1) {
-                                                        _ptvalue *= fabs(sin(weights_theta[l].value*pi/180.0));
+                                                        _ptvalue *= fabs(cos(weights_theta[l].value*pi/180.0));
                                                 }
                                                 sum += _ptvalue;

sansmodels/src/c_models/bcc.cpp

@@ -54 +54 @@
  */
 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 b3_x, b3_y, b3_z, b1_x, b1_y, b2_x, b2_y;
   double q_z;
-  double alpha, cos_val_b3, cos_val_b2, cos_val_b1;
+  double cos_val_b3, cos_val_b2, cos_val_b1;
   double a1_dot_q, a2_dot_q,a3_dot_q;
   double answer;
@@ -86 +86 @@
   ///  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);
+    b3_x = cos(theta) * cos(phi);
+    b3_y = sin(theta);
+    //b3_z = -cos(theta) * sin(phi);
+    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);
+    b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
+    b1_y = sin(psi)*cos(theta);
+    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);
-
+    b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
+        b2_y = cos(theta)*cos(psi);
+    cos_val_b2 = b2_x*q_x + b2_y*q_y;
+    
+    // The following test should always pass
+    if (fabs(cos_val_b3)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b3 = 1.0;
+    }
+    if (fabs(cos_val_b2)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b2 = 1.0;
+    }
+    if (fabs(cos_val_b1)>1.0) {
+      //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+      cos_val_b1 = 1.0;
+    }
     // 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);
@@ -111 +123 @@
     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));
@@ -287 +295 @@
                                         }
                                         if (weights_theta.size()>1) {
-                                                _ptvalue *= fabs(sin(weights_theta[j].value*pi/180.0));
+                                                _ptvalue *= fabs(cos(weights_theta[j].value*pi/180.0));
                                         }
                                         sum += _ptvalue;

sansmodels/src/c_models/capcyl.cpp

@@ -136 +136 @@
  */
 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 cyl_x, cyl_y;//, cyl_z;
+  //double q_z;
   double alpha, cos_val;
   double answer;
@@ -157 +157 @@
   //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);
+    cyl_x = cos(theta) * cos(phi);
+    cyl_y = sin(theta);
+    //cyl_z = -cos(theta) * sin(phi);
 
     // q vector
-    q_z = 0;
+    //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;
+    cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
 
     // The following test should always pass
@@ -348 +348 @@
                                                 }
                                                 if (weights_theta.size()>1) {
-                                                        _ptvalue *= fabs(sin(weights_theta[l].value*pi/180.0));
+                                                        _ptvalue *= fabs(cos(weights_theta[l].value*pi/180.0));
                                                 }
                                                 sum += _ptvalue;

sansmodels/src/c_models/corefourshell.cpp

@@ -28 +28 @@
 extern "C" {
         #include "libSphere.h"
+        #include "libmultifunc/libfunc.h"
 }
 
@@ -44 +45 @@
   double sld_solv;
   double background;
+  double M0_sld_shell1;
+  double M_theta_shell1;
+  double M_phi_shell1;
+  double M0_sld_shell2;
+  double M_theta_shell2;
+  double M_phi_shell2;
+  double M0_sld_shell3;
+  double M_theta_shell3;
+  double M_phi_shell3;
+  double M0_sld_shell4;
+  double M_theta_shell4;
+  double M_phi_shell4;
+  double M0_sld_core0;
+  double M_theta_core0;
+  double M_phi_core0;
+  double M0_sld_solv;
+  double M_theta_solv;
+  double M_phi_solv;
+  double Up_frac_i;
+  double Up_frac_f;
+  double Up_theta;
 } CoreFourShellParameters;
 
@@ -65 +87 @@
         sld_solv   = Parameter(6.4e-6);
         background = Parameter(0.001);
+        M0_sld_shell1 = Parameter(0.0e-6);
+        M_theta_shell1 = Parameter(0.0);
+        M_phi_shell1 = Parameter(0.0);
+        M0_sld_shell2 = Parameter(0.0e-6);
+        M_theta_shell2 = Parameter(0.0);
+        M_phi_shell2 = Parameter(0.0);
+        M0_sld_shell3 = Parameter(0.0e-6);
+        M_theta_shell3 = Parameter(0.0);
+        M_phi_shell3 = Parameter(0.0);
+        M0_sld_shell4 = Parameter(0.0e-6);
+        M_theta_shell4 = Parameter(0.0);
+        M_phi_shell4 = Parameter(0.0);
+        M0_sld_core0 = Parameter(0.0e-6);
+        M_theta_core0 = Parameter(0.0);
+        M_phi_core0 = Parameter(0.0);
+        M0_sld_solv = Parameter(0.0e-6);
+        M_theta_solv = Parameter(0.0);
+        M_phi_solv = Parameter(0.0);
+        Up_frac_i = Parameter(0.5);
+        Up_frac_f = Parameter(0.5);
+        Up_theta = Parameter(0.0);
+
 }
 
@@ -78 +122 @@
         // Fill parameter array for IGOR library
         // Add the background after averaging
+
         dp[0] = scale();
         dp[1] = rad_core0();
@@ -158 +203 @@
  * @return: function value
  */
+
+static double corefourshell_analytical_2D_scaled(CoreFourShellParameters *pars, double q, double q_x, double q_y) {
+        double dp[13];
+
+        // Fill parameter array for IGOR library
+        // Add the background after averaging
+
+        dp[0] = pars->scale;
+        dp[1] = pars->rad_core0;
+        dp[2] = 0.0; //sld_core0;
+        dp[3] = pars->thick_shell1;
+        dp[4] = 0.0; //sld_shell1;
+        dp[5] = pars->thick_shell2;
+        dp[6] = 0.0; //sld_shell2;
+        dp[7] = pars->thick_shell3;
+        dp[8] = 0.0; //sld_shell3;
+        dp[9] = pars->thick_shell4;
+        dp[10] = 0.0; //sld_shell4;
+        dp[11] = 0.0; //sld_solv;
+        dp[12] = 0.0;
+
+    double sld_core0 = pars->sld_core0;
+    double sld_shell1 = pars->sld_shell1;
+        double sld_shell2 = pars->sld_shell2;
+        double sld_shell3 = pars->sld_shell3;
+        double sld_shell4 = pars->sld_shell4;
+        double sld_solv = pars->sld_solv;
+        double answer = 0.0;
+        double m_max0 = pars->M0_sld_core0;
+        double m_max_shell1 = pars->M0_sld_shell1;
+        double m_max_shell2 = pars->M0_sld_shell2;
+        double m_max_shell3 = pars->M0_sld_shell3;
+        double m_max_shell4 = pars->M0_sld_shell4;
+        double m_max_solv = pars->M0_sld_solv;
+
+        if (m_max0 < 1.0e-32 && m_max_solv < 1.0e-32 && m_max_shell1 < 1.0e-32 && m_max_shell2 <
+                                                1.0e-32 && m_max_shell3 < 1.0e-32 && m_max_shell4 < 1.0e-32){
+                dp[2] = sld_core0;
+                dp[4] = sld_shell1;
+                dp[6] = sld_shell2;
+                dp[8] = sld_shell3;
+                dp[10] = sld_shell4;
+                dp[11] = sld_solv;
+                answer = FourShell(dp, q);
+        }
+        else{
+                double qx = q_x;
+                double qy = q_y;
+                double s_theta = pars->Up_theta;
+                double m_phi0 = pars->M_phi_core0;
+                double m_theta0 = pars->M_theta_core0;
+                double m_phi_shell1 = pars->M_phi_shell1;
+                double m_theta_shell1 = pars->M_theta_shell1;
+                double m_phi_shell2 = pars->M_phi_shell2;
+                double m_theta_shell2 = pars->M_theta_shell2;
+                double m_phi_shell3 = pars->M_phi_shell3;
+                double m_theta_shell3 = pars->M_theta_shell3;
+                double m_phi_shell4 = pars->M_phi_shell4;
+                double m_theta_shell4 = pars->M_theta_shell4;
+                double m_phi_solv = pars->M_phi_solv;
+                double m_theta_solv = pars->M_theta_solv;
+                double in_spin = pars->Up_frac_i;
+                double out_spin = pars->Up_frac_f;
+                polar_sld p_sld_core0;
+                polar_sld p_sld_shell1;
+                polar_sld p_sld_shell2;
+                polar_sld p_sld_shell3;
+                polar_sld p_sld_shell4;
+                polar_sld p_sld_solv;
+                //Find (b+m) slds
+                p_sld_core0 = cal_msld(1, qx, qy, sld_core0, m_max0, m_theta0, m_phi0,
+                                                                in_spin, out_spin, s_theta);
+                p_sld_shell1 = cal_msld(1, qx, qy, sld_shell1, m_max_shell1, m_theta_shell1, m_phi_shell1,
+                                                                in_spin, out_spin, s_theta);
+                p_sld_shell2 = cal_msld(1, qx, qy, sld_shell2, m_max_shell2, m_theta_shell2, m_phi_shell2,
+                                                                in_spin, out_spin, s_theta);
+                p_sld_shell3 = cal_msld(1, qx, qy, sld_shell3, m_max_shell3, m_theta_shell3, m_phi_shell3,
+                                                                in_spin, out_spin, s_theta);
+                p_sld_shell4 = cal_msld(1, qx, qy, sld_shell4, m_max_shell4, m_theta_shell4, m_phi_shell4,
+                                                                in_spin, out_spin, s_theta);
+                p_sld_solv = cal_msld(1, qx, qy, sld_solv, m_max_solv, m_theta_solv, m_phi_solv,
+                                                                in_spin, out_spin, s_theta);
+                //up_up
+                if (in_spin > 0.0 && out_spin > 0.0){
+                        dp[2] = p_sld_core0.uu;
+                        dp[4] = p_sld_shell1.uu;
+                        dp[6] = p_sld_shell2.uu;
+                        dp[8] = p_sld_shell3.uu;
+                        dp[10] = p_sld_shell4.uu;
+                        dp[11] = p_sld_solv.uu;
+                        answer += FourShell(dp, q);
+                        }
+                //down_down
+                if (in_spin < 1.0 && out_spin < 1.0){
+                        dp[2] = p_sld_core0.dd;
+                        dp[4] = p_sld_shell1.dd;
+                        dp[6] = p_sld_shell2.dd;
+                        dp[8] = p_sld_shell3.dd;
+                        dp[10] = p_sld_shell4.dd;
+                        dp[11] = p_sld_solv.dd;
+                        answer += FourShell(dp, q);
+                        }
+                //up_down
+                if (in_spin > 0.0 && out_spin < 1.0){
+                        dp[2] = p_sld_core0.re_ud;
+                        dp[4] = p_sld_shell1.re_ud;
+                        dp[6] = p_sld_shell2.re_ud;
+                        dp[8] = p_sld_shell3.re_ud;
+                        dp[10] = p_sld_shell4.re_ud;
+                        dp[11] = p_sld_solv.re_ud;
+                        answer += FourShell(dp, q);
+                        dp[2] = p_sld_core0.im_ud;
+                        dp[4] = p_sld_shell1.im_ud;
+                        dp[6] = p_sld_shell2.im_ud;
+                        dp[8] = p_sld_shell3.im_ud;
+                        dp[10] = p_sld_shell4.im_ud;
+                        dp[11] = p_sld_solv.im_ud;
+                        answer += FourShell(dp, q);
+                        }
+                //down_up
+                if (in_spin < 1.0 && out_spin > 0.0){
+                        dp[2] = p_sld_core0.re_du;
+                        dp[4] = p_sld_shell1.re_du;
+                        dp[6] = p_sld_shell2.re_du;
+                        dp[8] = p_sld_shell3.re_du;
+                        dp[10] = p_sld_shell4.re_du;
+                        dp[11] = p_sld_solv.re_du;
+                        answer += FourShell(dp, q);
+                        dp[2] = p_sld_core0.im_du;
+                        dp[4] = p_sld_shell1.im_du;
+                        dp[6] = p_sld_shell2.im_du;
+                        dp[8] = p_sld_shell3.im_du;
+                        dp[10] = p_sld_shell4.im_du;
+                        dp[11] = p_sld_solv.im_du;
+                        answer += FourShell(dp, q);
+                        }
+        }
+        // Already normalized
+        // add in the background
+        answer += pars->background;
+        return answer;
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters
+ * @param q: q-value
+ * @return: function value
+ */
+static double corefourshell_analytical_2DXY(CoreFourShellParameters *pars, double qx, double qy) {
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return corefourshell_analytical_2D_scaled(pars, q, qx/q, qy/q);
+}
+
+
 double CoreFourShellModel :: operator()(double qx, double qy) {
-        double q = sqrt(qx*qx + qy*qy);
-        return (*this).operator()(q);
+        CoreFourShellParameters dp;
+        dp.scale = scale();
+        dp.rad_core0 = rad_core0();
+        dp.sld_core0 = sld_core0();
+        dp.thick_shell1 = thick_shell1();
+        dp.sld_shell1 = sld_shell1();
+        dp.thick_shell2 = thick_shell2();
+        dp.sld_shell2 = sld_shell2();
+        dp.thick_shell3 = thick_shell3();
+        dp.sld_shell3 = sld_shell3();
+        dp.thick_shell4 = thick_shell4();
+        dp.sld_shell4 = sld_shell4();
+        dp.sld_solv = sld_solv();
+        dp.background = 0.0;
+        dp.M0_sld_shell1 = M0_sld_shell1();
+        dp.M_theta_shell1 = M_theta_shell1();
+        dp.M_phi_shell1 = M_phi_shell1();
+        dp.M0_sld_shell2 = M0_sld_shell2();
+        dp.M_theta_shell2 = M_theta_shell2();
+        dp.M_phi_shell2 = M_phi_shell2();
+        dp.M0_sld_shell3 = M0_sld_shell3();
+        dp.M_theta_shell3 = M_theta_shell3();
+        dp.M_phi_shell3 = M_phi_shell3();
+        dp.M0_sld_shell4 = M0_sld_shell4();
+        dp.M_theta_shell4 = M_theta_shell4();
+        dp.M_phi_shell4 = M_phi_shell4();
+        dp.M0_sld_core0 = M0_sld_core0();
+        dp.M_theta_core0 = M_theta_core0();
+        dp.M_phi_core0 = M_phi_core0();
+        dp.M0_sld_solv = M0_sld_solv();
+        dp.M_theta_solv = M_theta_solv();
+        dp.M_phi_solv = M_phi_solv();
+        dp.Up_frac_i = Up_frac_i();
+        dp.Up_frac_f = Up_frac_f();
+        dp.Up_theta = Up_theta();
+    
+        // Get the dispersion points for the radius
+        vector<WeightPoint> weights_rad;
+        rad_core0.get_weights(weights_rad);
+
+        // Get the dispersion points for the thick 1
+        vector<WeightPoint> weights_s1;
+        thick_shell1.get_weights(weights_s1);
+
+        // Get the dispersion points for the thick 2
+        vector<WeightPoint> weights_s2;
+        thick_shell2.get_weights(weights_s2);
+
+        // Get the dispersion points for the thick 3
+        vector<WeightPoint> weights_s3;
+        thick_shell3.get_weights(weights_s3);
+
+        // Get the dispersion points for the thick 4
+        vector<WeightPoint> weights_s4;
+        thick_shell4.get_weights(weights_s4);
+
+        // 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.rad_core0 = weights_rad[i].value;
+                // Loop over radius weight points
+                for(size_t j=0; j<weights_s1.size(); j++) {
+                        dp.thick_shell1 = weights_s1[j].value;
+                        // Loop over radius weight points
+                        for(size_t k=0; k<weights_s2.size(); k++) {
+                                dp.thick_shell2 = weights_s2[k].value;
+                                // Loop over radius weight points
+                                for(size_t l=0; l<weights_s3.size(); l++) {
+                                        dp.thick_shell3 = weights_s3[l].value;
+                                        // Loop over radius weight points
+                                        for(size_t m=0; m<weights_s4.size(); m++) {
+                                                dp.thick_shell4 = weights_s4[m].value;
+                                                //Un-normalize FourShell by volume
+                                                sum += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
+                                                        * corefourshell_analytical_2DXY(&dp, qx, qy) * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value),3);
+                                                //Find average volume
+                                                vol += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight
+                                                        * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value),3);
+
+                                                norm += weights_rad[i].weight*weights_s1[j].weight*weights_s2[k].weight*weights_s3[l].weight*weights_s4[m].weight;
+                                        }
+                                }
+                        }
+                }
+        }
+
+        if (vol != 0.0 && norm != 0.0) {
+                //Re-normalize by avg volume
+                sum = sum/(vol/norm);}
+        return sum/norm + background();
 }
 
@@ -171 +464 @@
  */
 double CoreFourShellModel :: evaluate_rphi(double q, double phi) {
-        return (*this).operator()(q);
+        double qx = q*cos(phi);
+        double qy = q*sin(phi);
+        return (*this).operator()(qx, qy);
 }
 

sansmodels/src/c_models/coreshellbicelle.cpp

@@ -56 +56 @@
  */
 static double core_shell_bicelle_analytical_2D_scaled(CoreShellBicelleParameters *pars, double q, double q_x, double q_y) {
-  double cyl_x, cyl_y, cyl_z;
-  double q_z;
+  double cyl_x, cyl_y;//, cyl_z;
+  //double q_z;
   double alpha, vol, cos_val;
   double answer;
@@ -66 +66 @@
 
   // 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;
+  //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;
+  cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
 
   // The following test should always pass
@@ -315 +315 @@
 
                                 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;

sansmodels/src/c_models/coreshellcylinder.cpp

@@ -54 +54 @@
  */
 static double core_shell_cylinder_analytical_2D_scaled(CoreShellCylinderParameters *pars, double q, double q_x, double q_y) {
-  double cyl_x, cyl_y, cyl_z;
-  double q_z;
+  double cyl_x, cyl_y;//, cyl_z;
+  //double q_z;
   double alpha, vol, cos_val;
   double answer;
@@ -64 +64 @@
 
   // 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;
+  //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;
+  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;
+    cos_val = 1.0;
   }
 
   alpha = acos( cos_val );
-
+  if (alpha == 0.0){
+        alpha = 1.0e-26;
+  }
   // Call the IGOR library function to get the kernel
   answer = CoreShellCylKernel(q, pars->radius, pars->thickness,
@@ -285 +287 @@
 
             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;

sansmodels/src/c_models/coreshellsphere.cpp

@@ -28 +28 @@
 extern "C" {
 #include "libSphere.h"
+#include "libmultifunc/libfunc.h"
 }
 
@@ -38 +39 @@
   double solvent_sld;
   double background;
+  double M0_sld_shell;
+  double M_theta_shell;
+  double M_phi_shell;
+  double M0_sld_core;
+  double M_theta_core;
+  double M_phi_core;
+  double M0_sld_solv;
+  double M_theta_solv;
+  double M_phi_solv;
+  double Up_frac_i;
+  double Up_frac_f;
+  double Up_theta;
 } CoreShellParameters;
 
@@ -50 +63 @@
   solvent_sld = Parameter(3.e-6);
   background = Parameter(0.0);
+  M0_sld_shell = Parameter(0.0e-6);
+  M_theta_shell = Parameter(0.0);
+  M_phi_shell = Parameter(0.0);
+  M0_sld_core = Parameter(0.0e-6);
+  M_theta_core = Parameter(0.0);
+  M_phi_core = Parameter(0.0);
+  M0_sld_solv = Parameter(0.0e-6);
+  M_theta_solv = Parameter(0.0);
+  M_phi_solv = Parameter(0.0);
+  Up_frac_i = Parameter(0.5);
+  Up_frac_f = Parameter(0.5);
+  Up_theta = Parameter(0.0);
 }
 
@@ -72 +97 @@
   dp[6] = 0.0;
 
+  //im
+  ///dp[7] = 0.0;
+  ///dp[8] = 0.0;
+  ///dp[9] = 0.0;
 
   // Get the dispersion points for the radius
@@ -112 +141 @@
 }
 
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+
+static double coreshell_analytical_2D_scaled(CoreShellParameters *pars, double q, double q_x, double q_y) {
+        double dp[7];
+        //convert angle degree to radian
+        dp[0] = pars->scale;
+        dp[1] = pars->radius;
+        dp[2] = pars->thickness;
+        dp[3] = 0.0;
+        dp[4] = 0.0;
+        dp[5] = 0.0;
+        dp[6] = 0.0;
+        //im
+        ///dp[7] = 0.0;
+        ///dp[8] = 0.0;
+        ///dp[9] = 0.0;
+
+    double sld_core = pars->core_sld;
+    double sld_shell = pars->shell_sld;
+        double sld_solv = pars->solvent_sld;
+        double answer = 0.0;
+        double m_max = pars->M0_sld_core;
+        double m_max_shell = pars->M0_sld_shell;
+        double m_max_solv = pars->M0_sld_solv;
+
+        if (m_max < 1.0e-32 && m_max_solv < 1.0e-32 && m_max_shell < 1.0e-32){
+                dp[3] = sld_core;
+                dp[4] = sld_shell;
+                dp[5] = sld_solv;
+                answer = CoreShellForm(dp, q);
+        }
+        else{
+                double qx = q_x;
+                double qy = q_y;
+                double s_theta = pars->Up_theta;
+                double m_phi = pars->M_phi_core;
+                double m_theta = pars->M_theta_core;
+                double m_phi_shell = pars->M_phi_shell;
+                double m_theta_shell = pars->M_theta_shell;
+                double m_phi_solv = pars->M_phi_solv;
+                double m_theta_solv = pars->M_theta_solv;
+                double in_spin = pars->Up_frac_i;
+                double out_spin = pars->Up_frac_f;
+                polar_sld p_sld_core;
+                polar_sld p_sld_shell;
+                polar_sld p_sld_solv;
+                //Find (b+m) slds
+                p_sld_core = cal_msld(1, qx, qy, sld_core, m_max, m_theta, m_phi,
+                                                                in_spin, out_spin, s_theta);
+                p_sld_shell = cal_msld(1, qx, qy, sld_shell, m_max_shell, m_theta_shell, m_phi_shell,
+                                                                in_spin, out_spin, s_theta);
+                p_sld_solv = cal_msld(1, qx, qy, sld_solv, m_max_solv, m_theta_solv, m_phi_solv,
+                                                                in_spin, out_spin, s_theta);
+                //up_up
+                if (in_spin > 0.0 && out_spin > 0.0){
+                        dp[3] = p_sld_core.uu;
+                        dp[4] = p_sld_shell.uu;
+                        dp[5] = p_sld_solv.uu;
+                        answer += CoreShellForm(dp, q);
+                        }
+                //down_down
+                if (in_spin < 1.0 && out_spin < 1.0){
+                        dp[3] = p_sld_core.dd;
+                        dp[4] = p_sld_shell.dd;
+                        dp[5] = p_sld_solv.dd;
+                        answer += CoreShellForm(dp, q);
+                        }
+                //up_down
+                if (in_spin > 0.0 && out_spin < 1.0){
+                        dp[3] = p_sld_core.re_ud;
+                        dp[4] = p_sld_shell.re_ud;
+                        dp[5] = p_sld_solv.re_ud;
+                        answer += CoreShellForm(dp, q);
+                        dp[3] = p_sld_core.im_ud;
+                        dp[4] = p_sld_shell.im_ud;
+                        dp[5] = p_sld_solv.im_ud;
+                        answer += CoreShellForm(dp, q);
+                        }
+                //down_up
+                if (in_spin < 1.0 && out_spin > 0.0){
+                        dp[3] = p_sld_core.re_du;
+                        dp[4] = p_sld_shell.re_du;
+                        dp[5] = p_sld_solv.re_du;
+                        answer += CoreShellForm(dp, q);
+                        dp[3] = p_sld_core.im_du;
+                        dp[4] = p_sld_shell.im_du;
+                        dp[5] = p_sld_solv.im_du;
+                        answer += CoreShellForm(dp, q);
+                        }
+        }
+        // Already normalized
+        // add in the background
+        answer += pars->background;
+        return answer;
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters
+ * @param q: q-value
+ * @return: function value
+ */
+static double coreshell_analytical_2DXY(CoreShellParameters *pars, double qx, double qy) {
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return coreshell_analytical_2D_scaled(pars, q, qx/q, qy/q);
+}
+
+
 /**
  * Function to evaluate 2D scattering function
@@ -119 +266 @@
  */
 double CoreShellModel :: operator()(double qx, double qy) {
-  double q = sqrt(qx*qx + qy*qy);
-  return (*this).operator()(q);
+        CoreShellParameters dp;
+    dp.scale = scale();
+    dp.radius = radius();
+    dp.thickness = thickness();
+    dp.core_sld = core_sld();
+    dp.shell_sld = shell_sld();
+    dp.solvent_sld = solvent_sld();
+    dp.background = 0.0;
+    dp.M0_sld_shell = M0_sld_shell();
+    dp.M_theta_shell = M_theta_shell();
+    dp.M_phi_shell = M_phi_shell();
+    dp.M0_sld_core = M0_sld_core();
+    dp.M_theta_core = M_theta_core();
+    dp.M_phi_core = M_phi_core();
+    dp.M0_sld_solv = M0_sld_solv();
+    dp.M_theta_solv = M_theta_solv();
+    dp.M_phi_solv = M_phi_solv();
+    dp.Up_frac_i = Up_frac_i();
+    dp.Up_frac_f = Up_frac_f();
+    dp.Up_theta = Up_theta();
+
+        // 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);
+
+        // 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.radius = weights_rad[i].value;
+
+        // Loop over thickness weight points
+        for(size_t j=0; j<weights_thick.size(); j++) {
+          dp.thickness = weights_thick[j].value;
+          //Un-normalize SphereForm by volume
+          sum += weights_rad[i].weight
+                  * weights_thick[j].weight * coreshell_analytical_2DXY(&dp, qx, qy) *
+                        pow(weights_rad[i].value+weights_thick[j].value,3);
+
+          //Find average volume
+          vol += weights_rad[i].weight * weights_thick[j].weight
+                  * pow(weights_rad[i].value+weights_thick[j].value,3);
+          norm += weights_rad[i].weight
+                  * weights_thick[j].weight;
+        }
+        }
+
+        if (vol != 0.0 && norm != 0.0) {
+                //Re-normalize by avg volume
+                sum = sum/(vol/norm);}
+        return sum/norm + background();
 }
 
@@ -131 +334 @@
  */
 double CoreShellModel :: evaluate_rphi(double q, double phi) {
-  return (*this).operator()(q);
+        double qx = q*cos(phi);
+        double qy = q*sin(phi);
+        return (*this).operator()(qx, qy);
 }
 /**

sansmodels/src/c_models/csparallelepiped.cpp

@@ -143 +143 @@
 static double csparallelepiped_analytical_2D_scaled(CSParallelepipedParameters *pars, double q, double q_x, double q_y) {
   double dp[13];
-  double cparallel_x, cparallel_y, cparallel_z, bparallel_x, bparallel_y, parallel_x, parallel_y;
-  double q_z;
-  double alpha, cos_val_c, cos_val_b, cos_val_a, edgeA, edgeB, edgeC;
+  double cparallel_x, cparallel_y, bparallel_x, bparallel_y, parallel_x, parallel_y;
+  //double q_z;
+  double cos_val_c, cos_val_b, cos_val_a, edgeA, edgeB, edgeC;
 
   double answer;
@@ -176 +176 @@
 
   // parallelepiped c axis orientation
-  cparallel_x = sin(theta) * cos(phi);
-  cparallel_y = sin(theta) * sin(phi);
-  cparallel_z = cos(theta);
+  cparallel_x = cos(theta) * cos(phi);
+  cparallel_y = sin(theta);
+  //cparallel_z = -cos(theta) * sin(phi);
 
   // q vector
-  q_z = 0.0;
+  //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);
+  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);
+  parallel_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
+  parallel_y = sin(psi)*cos(theta);
 
   cos_val_a = parallel_x*q_x + parallel_y*q_y;
@@ -197 +197 @@
 
   // 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);
+  bparallel_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
+  bparallel_y = cos(theta)*cos(psi);
   // axis of the parallelepiped
-  cos_val_b = sin(acos(cos_val_a)) ;
-
-
+  cos_val_b = bparallel_x*q_x + bparallel_y*q_y;
 
   // 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;
-  }
-
+    //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+    cos_val_c = 1.0;
+  }
+  if (fabs(cos_val_a)>1.0) {
+    //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+    cos_val_a = 1.0;
+  }
+  if (fabs(cos_val_b)>1.0) {
+    //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+    cos_val_b = 1.0;
+  }
   // Call the IGOR library function to get the kernel
-  answer = cspkernel( dp,q, sin(alpha)*cos_val_a,sin(alpha)*cos_val_b,cos_val_c);
+  answer = cspkernel( dp, q, cos_val_a, cos_val_b, cos_val_c);
 
   //convert to [cm-1]
@@ -434 +439 @@
 
               if (weights_parallel_theta.size()>1) {
-                _ptvalue *= fabs(sin(weights_parallel_theta[l].value*pi/180.0));
+                _ptvalue *= fabs(cos(weights_parallel_theta[l].value*pi/180.0));
               }
               sum += _ptvalue;

sansmodels/src/c_models/cylinder.cpp

@@ -28 +28 @@
         #include "libCylinder.h"
         #include "libStructureFactor.h"
+        #include "libmultifunc/libfunc.h"
 }
 #include "cylinder.h"
@@ -41 +42 @@
     double cyl_theta;
     double cyl_phi;
+    double M0_sld_cyl;
+    double M_theta_cyl;
+    double M_phi_cyl;
+    double M0_sld_solv;
+    double M_theta_solv;
+    double M_phi_solv;
+    double Up_frac_i;
+        double Up_frac_f;
+        double Up_theta;
 } CylinderParameters;
 
@@ -54 +64 @@
         cyl_theta  = Parameter(0.0, true);
         cyl_phi    = Parameter(0.0, true);
+        M0_sld_cyl = Parameter(0.0e-6);
+        M_theta_cyl = Parameter(0.0);
+        M_phi_cyl = Parameter(0.0); 
+        M0_sld_solv = Parameter(0.0e-6);
+        M_theta_solv = Parameter(0.0);
+        M_phi_solv = Parameter(0.0); 
+        Up_frac_i = Parameter(0.5); 
+        Up_frac_f = Parameter(0.5);
+        Up_theta = Parameter(0.0);
 }
 
@@ -122 +141 @@
  */
 static double cylinder_analytical_2D_scaled(CylinderParameters *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->cyl_theta * pi/180.0;
-  double phi = pars->cyl_phi * pi/180.0;
+    double cyl_x, cyl_y;//, cyl_z;
+    //double q_z;
+    double alpha, vol, cos_val;
+    double answer = 0.0;
+    double form = 0.0;
+    //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 sld_solv = pars->sldSolv;
+    double sld_cyl = pars->sldCyl;
+    double m_max = pars->M0_sld_cyl;
+    double m_max_solv = pars->M0_sld_solv;
+    double contrast = 0.0;
 
     // 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;
 
     // 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
@@ -157 +181 @@
         }
   // Call the IGOR library function to get the kernel
-  answer = CylKernel(q, pars->radius, pars->length/2.0, alpha) / sin(alpha);
-
-  // 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
+  //answer = CylKernel(q, pars->radius, pars->length/2.0, alpha) / sin(alpha);
+
+    // Call the IGOR library function to get the kernel
+    form = CylKernel(q, pars->radius, pars->length/2.0, alpha) / sin(alpha);
+
+        if (m_max < 1.0e-32 && m_max_solv < 1.0e-32){
+                // Multiply by contrast^2
+                contrast = (pars->sldCyl - pars->sldSolv);
+                answer = contrast * contrast * form;
+        }
+        else{
+                double qx = q_x;
+                double qy = q_y;
+                double s_theta = pars->Up_theta;
+                double m_phi = pars->M_phi_cyl;
+                double m_theta = pars->M_theta_cyl;
+                double m_phi_solv = pars->M_phi_solv;
+                double m_theta_solv = pars->M_theta_solv;
+                double in_spin = pars->Up_frac_i;
+                double out_spin = pars->Up_frac_f;
+                polar_sld p_sld;
+                polar_sld p_sld_solv;
+                p_sld = cal_msld(1, qx, qy, sld_cyl, m_max, m_theta, m_phi, 
+                                                in_spin, out_spin, s_theta);
+                p_sld_solv = cal_msld(1, qx, qy, sld_solv, m_max_solv, m_theta_solv, m_phi_solv, 
+                                                in_spin, out_spin, s_theta);
+                //up_up 
+                if (in_spin > 0.0 && out_spin > 0.0){                    
+                        answer += ((p_sld.uu- p_sld_solv.uu) * (p_sld.uu- p_sld_solv.uu) * form);
+                        }
+                //down_down
+                if (in_spin < 1.0 && out_spin < 1.0){
+                        answer += ((p_sld.dd - p_sld_solv.dd) * (p_sld.dd - p_sld_solv.dd) * form);
+                        }
+                //up_down
+                if (in_spin > 0.0 && out_spin < 1.0){
+                        answer += ((p_sld.re_ud - p_sld_solv.re_ud) * (p_sld.re_ud - p_sld_solv.re_ud) * form);
+                        answer += ((p_sld.im_ud - p_sld_solv.im_ud) * (p_sld.im_ud - p_sld_solv.im_ud) * form);
+                        }
+                //down_up       
+                if (in_spin < 1.0 && out_spin > 0.0){
+                        answer += ((p_sld.re_du - p_sld_solv.re_du) * (p_sld.re_du - p_sld_solv.re_du) * form);
+                        answer += ((p_sld.im_du - p_sld_solv.im_du) * (p_sld.im_du - p_sld_solv.im_du) * form);
+                        }
+        }
+
+    //normalize by cylinder volume
+    //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
     vol = acos(-1.0) * pars->radius * pars->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;
+    answer *= vol;
+
+    //convert to [cm-1]
+    answer *= 1.0e8;
+
+    //Scale
+    answer *= pars->scale;
+
+    // add in the background
+    answer += pars->background;
+
+    return answer;
 }
 
@@ -208 +273 @@
         dp.cyl_theta  = cyl_theta();
         dp.cyl_phi    = cyl_phi();
-
+        dp.Up_theta =  Up_theta();
+        dp.M_phi_cyl =  M_phi_cyl();
+        dp.M_theta_cyl =  M_theta_cyl();
+        dp.M0_sld_cyl =  M0_sld_cyl();
+        dp.M_phi_solv =  M_phi_solv();
+        dp.M_theta_solv =  M_theta_solv();
+        dp.M0_sld_solv =  M0_sld_solv();
+        dp.Up_frac_i =  Up_frac_i();
+        dp.Up_frac_f =  Up_frac_f();
+        
         // Get the dispersion points for the radius
         vector<WeightPoint> weights_rad;
@@ -255 +329 @@
                                                 *pow(weights_rad[i].value,2)*weights_len[j].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;

sansmodels/src/c_models/ellipsoid.cpp

@@ -51 +51 @@
  */
 static double ellipsoid_analytical_2D_scaled(EllipsoidParameters *pars, double q, double q_x, double q_y) {
-  double cyl_x, cyl_y, cyl_z;
-  double q_z;
+  double cyl_x, cyl_y;//, cyl_z;
+  //double q_z;
   double alpha, vol, cos_val;
   double answer;
@@ -61 +61 @@
 
   // Ellipsoid 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;
 
   // 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
@@ -252 +252 @@
           * 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));
+            _ptvalue *= fabs(cos(weights_theta[k].value*pi/180.0));
           }
           sum += _ptvalue;

sansmodels/src/c_models/ellipticalcylinder.cpp

@@ -46 +46 @@
 
 
-static double elliptical_cylinder_kernel(EllipticalCylinderParameters *pars, double q, double alpha, double nu) {
+static double elliptical_cylinder_kernel(EllipticalCylinderParameters *pars, double q, double cos_val, double cos_nu, double cos_mu) {
   double qr;
   double qL;
@@ -55 +55 @@
   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;
+  qr = q*sqrt( r_major*r_major*cos_nu*cos_nu + pars->r_minor*pars->r_minor*cos_mu*cos_mu );
+  qL = q*pars->length*cos_val/2.0;
 
   if (qr==0){
@@ -83 +83 @@
  */
 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 cyl_x, cyl_y;//, cyl_z;
+  double ella_x, ella_y, ellb_x, ellb_y;
+  //double q_z;
+  double vol, cos_val;
+  double cos_mu, cos_nu;
   double answer;
   //convert angle degree to radian
@@ -96 +96 @@
 
   //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;
-
-  // 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;
-  }
+  //q_z = 0;
 
   // Note: cos(alpha) = 0 and 1 will get an
   // undefined value from CylKernel
-  alpha = acos( cos_val );
+  //alpha = acos( cos_val );
 
   //ellipse orientation:
@@ -125 +115 @@
 
   //x- y- component on the detector plane.
-  ell_x =  cos(psi);
-  ell_y =  sin(psi);
+  ella_x =  -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
+  ella_y =  sin(psi)*cos(theta);
+  ellb_x =  -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
+  ellb_y =  cos(theta)*cos(psi);
+  
+  // 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;
 
   // 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("cyl_ana_2D: Unexpected error: cos(alpha)>1\n");
+    cos_val = 1.0;
+  }
   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);
+    //printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n");
+    cos_nu = 1.0;
+  }
+  if (fabs(cos_mu)>1.0) {
+    //printf("cyl_ana_2D: Unexpected error: cos(nu)>1\n");
+    cos_mu = 1.0;
+  }
+  answer = elliptical_cylinder_kernel(pars, q, cos_val, cos_nu, cos_mu);
 
   // Multiply by contrast^2
@@ -347 +350 @@
               * weights_rat[m].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;

sansmodels/src/c_models/fcc.cpp

@@ -53 +53 @@
  */
 static double fc_analytical_2D_scaled(FCParameters *pars, double q, double q_x, double q_y) {
-  double b3_x, b3_y, b3_z, b1_x, b1_y;
+  double b3_x, b3_y, b3_z, b1_x, b1_y, b2_x, b2_y;
   double q_z;
-  double alpha, cos_val_b3, cos_val_b2, cos_val_b1;
+  double cos_val_b3, cos_val_b2, cos_val_b1;
   double a1_dot_q, a2_dot_q,a3_dot_q;
   double answer;
@@ -83 +83 @@
   ///  instead of against q coordinate in 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);
+  b3_x = cos(theta) * cos(phi);
+  b3_y = sin(theta);
+  //b3_z = -cos(theta) * sin(phi);
+  cos_val_b3 =  b3_x*q_x + b3_y*q_y;// + b3_z*q_z;
+  
   // b1 axis orientation
-  b1_x = sin(psi);
-  b1_y = cos(psi);
-  cos_val_b1 = (b1_x*q_x + b1_y*q_y);
+  b1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
+  b1_y = sin(psi)*cos(theta);
+  cos_val_b1 = b1_x*q_x + b1_y*q_y;
+  
   // b2 axis orientation
-  cos_val_b2 = sin(acos(cos_val_b1));
-  // alpha correction
-  cos_val_b2 *= sin(alpha);
-  cos_val_b1 *= sin(alpha);
-
+  b2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
+  b2_y = cos(theta)*cos(psi);
+  cos_val_b2 =  b2_x*q_x + b2_y*q_y;
+  
+  // The following test should always pass
+  if (fabs(cos_val_b3)>1.0) {
+    //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+    cos_val_b3 = 1.0;
+  }
+  if (fabs(cos_val_b2)>1.0) {
+    //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+    cos_val_b2 = 1.0;
+  }
+  if (fabs(cos_val_b1)>1.0) {
+    //printf("bcc_ana_2D: Unexpected error: cos()>1\n");
+    cos_val_b1 = 1.0;
+  }
   // Compute the angle btw vector q and the a3 axis
   a3_dot_q = 0.5*aa*q*(cos_val_b2+cos_val_b1);
@@ -107 +120 @@
   a2_dot_q = 0.5*aa*q*(cos_val_b3+cos_val_b1);
 
-  // The following test should always pass
-  if (fabs(cos_val_b3)>1.0) {
-    printf("fcc_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));
@@ -282 +291 @@
           }
           if (weights_theta.size()>1) {
-            _ptvalue *= fabs(sin(weights_theta[j].value*pi/180.0));
+            _ptvalue *= fabs(cos(weights_theta[j].value*pi/180.0));
           }
           sum += _ptvalue;

sansmodels/src/c_models/hollowcylinder.cpp

@@ -53 +53 @@
 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 q_z;
   double  alpha,vol, cos_val;
   double answer;
@@ -62 +62 @@
 
   // 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;
+  //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;
+  cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
 
   // The following test should always pass
@@ -280 +280 @@
                 * weights_length[j].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;

sansmodels/src/c_models/libfunc.c

@@ -103 +103 @@
     return -tmp+log(2.5066282746310005*ser/x);
 }
+
+// calculate magnetic sld and return total sld
+// bn : contrast (not just sld of the layer)
+// m0: max mag of M; mtheta: angle from x-z plane;
+// mphi: angle (anti-clock-wise)of x-z projection(M) from x axis
+// spinfraci: the fraction of UP among UP+Down (before sample)
+// spinfracf: the fraction of UP among UP+Down (after sample and before detector)
+// spintheta: angle (anti-clock-wise) between neutron spin(up) and x axis
+// Note: all angles are in degrees.
+polar_sld cal_msld(int isangle, double qx, double qy, double bn,
+                                double m01, double mtheta1, double mphi1,
+                                double spinfraci, double spinfracf, double spintheta)
+{
+        //locals
+        double q_x = qx;
+        double q_y = qy;
+        double sld = bn;
+        int is_angle = isangle;
+        double pi = 4.0*atan(1.0);
+        double s_theta = spintheta * pi/180.0;
+        double m_max = m01;
+        double m_phi = mphi1;
+        double m_theta = mtheta1;
+        double in_spin = spinfraci;
+        double out_spin = spinfracf;
+
+        double m_perp = 0.0;
+        double m_perp_z = 0.0;
+        double m_perp_y = 0.0;
+        double m_perp_x = 0.0;
+        double m_sigma_x = 0.0;
+        double m_sigma_z = 0.0;
+        double m_sigma_y = 0.0;
+        double b_m = 0.0;
+        double q_angle = 0.0;
+        double mx = 0.0;
+        double my = 0.0;
+        double mz = 0.0;
+        polar_sld p_sld;
+        p_sld.uu = sld;
+        p_sld.dd = sld;
+        p_sld.re_ud = 0.0;
+        p_sld.im_ud = 0.0;
+        p_sld.re_du = 0.0;
+        p_sld.im_du = 0.0;
+
+        //No mag means no further calculation
+        if (isangle>0){
+                if (m_max < 1.0e-32){
+                        p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu;
+                        p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd;
+                        return p_sld;
+                }
+        }
+        else{
+                if (fabs(m_max)< 1.0e-32 && fabs(m_phi)< 1.0e-32 && fabs(m_theta)< 1.0e-32){
+                        p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu;
+                        p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd;
+                        return p_sld;
+                }
+        }
+
+        //These are needed because of the precision of inputs
+        if (in_spin < 0.0) in_spin = 0.0;
+        if (in_spin > 1.0) in_spin = 1.0;
+        if (out_spin < 0.0) out_spin = 0.0;
+        if (out_spin > 1.0) out_spin = 1.0;
+
+        if (q_x == 0.0) q_angle = pi / 2.0;
+        else q_angle = atan(q_y/q_x);
+        if (q_y < 0.0 && q_x < 0.0) q_angle -= pi;
+        else if (q_y > 0.0 && q_x < 0.0) q_angle += pi;
+
+        q_angle = pi/2.0 - q_angle;
+        if (q_angle > pi) q_angle -= 2.0 * pi;
+        else if (q_angle < -pi) q_angle += 2.0 * pi;
+
+        if (fabs(q_x) < 1.0e-16 && fabs(q_y) < 1.0e-16){
+                m_perp = 0.0;
+                }
+        else {
+                m_perp = m_max;
+                }
+        if (is_angle > 0){
+                m_phi *= pi/180.0;
+                m_theta *= pi/180.0;
+                mx = m_perp * cos(m_theta) * cos(m_phi);
+                my = m_perp * sin(m_theta);
+                mz = -(m_perp * cos(m_theta) * sin(m_phi));
+        }
+        else{
+                mx = m_perp;
+                my = m_phi;
+                mz = m_theta;
+        }
+        //ToDo: simplify these steps
+        // m_perp1 -m_perp2
+        m_perp_x = (mx) *  cos(q_angle);
+        m_perp_x -= (my) * sin(q_angle);
+        m_perp_y = m_perp_x;
+        m_perp_x *= cos(-q_angle);
+        m_perp_y *= sin(-q_angle);
+        m_perp_z = mz;
+
+        m_sigma_x = (m_perp_x * cos(-s_theta) - m_perp_y * sin(-s_theta));
+        m_sigma_y = (m_perp_x * sin(-s_theta) + m_perp_y * cos(-s_theta));
+        m_sigma_z = (m_perp_z);
+
+        //Find b
+        p_sld.uu -= m_sigma_x;
+        p_sld.dd += m_sigma_x;
+        p_sld.re_ud = m_sigma_y;
+        p_sld.re_du = m_sigma_y;
+        p_sld.im_ud = m_sigma_z;
+        p_sld.im_du = -m_sigma_z;
+
+        p_sld.uu = sqrt(sqrt(in_spin * out_spin)) * p_sld.uu;
+        p_sld.dd = sqrt(sqrt((1.0 - in_spin) * (1.0 - out_spin))) * p_sld.dd;
+
+        p_sld.re_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * p_sld.re_ud;
+        p_sld.im_ud = sqrt(sqrt(in_spin * (1.0 - out_spin))) * p_sld.im_ud;
+        p_sld.re_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * p_sld.re_du;
+        p_sld.im_du = sqrt(sqrt((1.0 - in_spin) * out_spin)) * p_sld.im_du;
+
+        return p_sld;
+}
+
 
 /** Modifications below by kieranrcampbell@gmail.com
@@ -131 +258 @@
     ap = a;
     del = sum = 1.0/a;
-    
+
     for(n=1;n<=ITMAX;n++) {
       ++ap;
@@ -150 +277 @@
 
 /**
-   Implements the incomplete gamma function Q(a,x) evaluated by its continued fraction 
-   representation 
+   Implements the incomplete gamma function Q(a,x) evaluated by its continued fraction
+   representation
 **/
 
@@ -196 +323 @@
 }
 
-/** 
+/**
     Implementation of the error function, erf(x)
 **/

sansmodels/src/c_models/parallelepiped.cpp

@@ -29 +29 @@
         #include "libCylinder.h"
         #include "libStructureFactor.h"
+        #include "libmultifunc/libfunc.h"
 }
 #include "parallelepiped.h"
@@ -44 +45 @@
     double parallel_phi;
     double parallel_psi;
+    double M0_sld_pipe;
+    double M_theta_pipe;
+    double M_phi_pipe;
+    double M0_sld_solv;
+    double M_theta_solv;
+    double M_phi_solv;
+    double Up_frac_i;
+        double Up_frac_f;
+        double Up_theta;
 } ParallelepipedParameters;
 
@@ -85 +95 @@
  */
 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 cparallel_x, cparallel_y, bparallel_x, bparallel_y, parallel_x, parallel_y;
+  //double q_z;
+  double vol, cos_val_c, cos_val_b, cos_val_a, edgeA, edgeB, edgeC;
+
+  double answer = 0.0;
+  double form = 0.0;
   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;
-
+  double sld_solv = pars->sldSolv;
+  double sld_pipe = pars->sldPipe;
+  double m_max = pars->M0_sld_pipe;
+  double m_max_solv = pars->M0_sld_solv;
+  double contrast = 0.0;
+    
   edgeA = pars->short_a;
   edgeB = pars->short_b;
@@ -103 +119 @@
 
     // parallelepiped c axis orientation
-    cparallel_x = sin(theta) * cos(phi);
-    cparallel_y = sin(theta) * sin(phi);
-    cparallel_z = cos(theta);
+    cparallel_x = cos(theta) * cos(phi);
+    cparallel_y = sin(theta);
+    //cparallel_z = -cos(theta) * sin(phi);
 
     // q vector
-    q_z = 0.0;
+    //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);
+    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);
-
+    parallel_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
+    parallel_y = sin(psi)*cos(theta);
+        //parallel_z = sin(theta)*sin(phi)*sin(psi)+cos(phi)*cos(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);
+    bparallel_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
+    bparallel_y = cos(theta)*cos(psi);
+    //bparallel_z = sin(theta)*sin(phi)*cos(psi)-sin(psi)*cos(phi);
+    
     // axis of the parallelepiped
-    cos_val_b = sin(acos(cos_val_a)) ;
-
-
+    cos_val_b = bparallel_x*q_x + bparallel_y*q_y;
 
     // 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;
+      //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+      cos_val_c = 1.0;
     }
-
+    if (fabs(cos_val_a)>1.0) {
+      //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+      cos_val_a = 1.0;
+    }
+    if (fabs(cos_val_b)>1.0) {
+      //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+      cos_val_b = 1.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);
+  form = pkernel( q*edgeA, q*edgeB, q*edgeC, cos_val_a, cos_val_b, cos_val_c);
+  
+  if (m_max < 1.0e-32 && m_max_solv < 1.0e-32){
+      // Multiply by contrast^2
+      contrast = (pars->sldPipe - pars->sldSolv);
+          answer = contrast * contrast * form;
+  }
+  else{
+          double qx = q_x;
+          double qy = q_y;
+          double s_theta = pars->Up_theta;
+          double m_phi = pars->M_phi_pipe;
+          double m_theta = pars->M_theta_pipe;
+          double m_phi_solv = pars->M_phi_solv;
+          double m_theta_solv = pars->M_theta_solv;
+          double in_spin = pars->Up_frac_i;
+          double out_spin = pars->Up_frac_f;
+          polar_sld p_sld;
+          polar_sld p_sld_solv;
+          p_sld = cal_msld(1, qx, qy, sld_pipe, m_max, m_theta, m_phi, 
+                                        in_spin, out_spin, s_theta);
+          p_sld_solv = cal_msld(1, qx, qy, sld_solv, m_max_solv, m_theta_solv, m_phi_solv, 
+                                                in_spin, out_spin, s_theta);
+          //up_up       
+          if (in_spin > 0.0 && out_spin > 0.0){                  
+                  answer += ((p_sld.uu- p_sld_solv.uu) * (p_sld.uu- p_sld_solv.uu) * form);
+                  }
+          //down_down
+          if (in_spin < 1.0 && out_spin < 1.0){
+                  answer += ((p_sld.dd - p_sld_solv.dd) * (p_sld.dd - p_sld_solv.dd) * form);
+                  }
+          //up_down
+          if (in_spin > 0.0 && out_spin < 1.0){
+                  answer += ((p_sld.re_ud - p_sld_solv.re_ud) * (p_sld.re_ud - p_sld_solv.re_ud) * form);
+                  answer += ((p_sld.im_ud - p_sld_solv.im_ud) * (p_sld.im_ud - p_sld_solv.im_ud) * form);
+                  }
+          //down_up     
+          if (in_spin < 1.0 && out_spin > 0.0){
+                  answer += ((p_sld.re_du - p_sld_solv.re_du) * (p_sld.re_du - p_sld_solv.re_du) * form);
+                  answer += ((p_sld.im_du - p_sld_solv.im_du) * (p_sld.im_du - p_sld_solv.im_du) * form);
+                  }
+  }
+  
 
   //normalize by cylinder volume
   //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vparallel
-    vol = edgeA* edgeB * edgeC;
+  vol = edgeA* edgeB * edgeC;
   answer *= vol;
 
@@ -186 +247 @@
         parallel_phi    = Parameter(0.0, true);
         parallel_psi    = Parameter(0.0, true);
+        M0_sld_pipe = Parameter(0.0e-6);
+        M_theta_pipe = Parameter(0.0);
+        M_phi_pipe = Parameter(0.0); 
+        M0_sld_solv = Parameter(0.0e-6);
+        M_theta_solv = Parameter(0.0);
+        M_phi_solv = Parameter(0.0); 
+        Up_frac_i = Parameter(0.5); 
+        Up_frac_f = Parameter(0.5);
+        Up_theta = Parameter(0.0);
 }
 
@@ -279 +349 @@
         dp.parallel_phi    = parallel_phi();
         dp.parallel_psi    = parallel_psi();
+        dp.Up_theta =  Up_theta();
+        dp.M_phi_pipe =  M_phi_pipe();
+        dp.M_theta_pipe =  M_theta_pipe();
+        dp.M0_sld_pipe =  M0_sld_pipe();
+        dp.M_phi_solv =  M_phi_solv();
+        dp.M_theta_solv =  M_theta_solv();
+        dp.M0_sld_solv =  M0_sld_solv();
+        dp.Up_frac_i =  Up_frac_i();
+        dp.Up_frac_f =  Up_frac_f();
 
 
@@ -346 +425 @@
 
                                                         if (weights_parallel_theta.size()>1) {
-                                                                _ptvalue *= fabs(sin(weights_parallel_theta[l].value*pi/180.0));
+                                                                _ptvalue *= fabs(cos(weights_parallel_theta[l].value*pi/180.0));
                                                         }
                                                         sum += _ptvalue;

sansmodels/src/c_models/sc.cpp

@@ -55 +55 @@
   double a3_x, a3_y, a3_z, a2_x, a2_y, a1_x, a1_y;
   double q_z;
-  double alpha, cos_val_a3, cos_val_a2, cos_val_a1;
+  double cos_val_a3, cos_val_a2, cos_val_a1;
   double a1_dot_q, a2_dot_q,a3_dot_q;
   double answer;
@@ -80 +80 @@
   /// Angles here are respect to detector coordinate instead of against q coordinate(PRB 36, 3, 1754)
   // a3 axis orientation
-  a3_x = sin(theta) * cos(phi);//negative sign here???
-  a3_y = sin(theta) * sin(phi);
-  a3_z = cos(theta);
-
+  a3_x = cos(theta) * cos(phi);
+  a3_y = sin(theta);
+  //a3_z = -cos(theta) * sin(phi);
   // q vector
   q_z = 0.0;
 
   // Compute the angle btw vector q and the a3 axis
-  cos_val_a3 = a3_x*q_x + a3_y*q_y + a3_z*q_z;
-  alpha = acos(cos_val_a3);
-  //alpha = acos(cos_val_a3);
-  a3_dot_q = aa*q*cos_val_a3;
+  cos_val_a3 = a3_x*q_x + a3_y*q_y;// + a3_z*q_z;
+
   // a1 axis orientation
-  a1_x = sin(psi);
-  a1_y = cos(psi);
+  a1_x = -cos(phi)*sin(psi) * sin(theta)+sin(phi)*cos(psi);
+  a1_y = sin(psi)*cos(theta);
 
   cos_val_a1 = a1_x*q_x + a1_y*q_y;
-  a1_dot_q = aa*q*cos_val_a1*sin(alpha);
 
   // a2 axis orientation
-  a2_x = sqrt(1.0-sin(theta)*cos(phi))*cos(psi);
-  a2_y = sqrt(1.0-sin(theta)*cos(phi))*sin(psi);
+  a2_x = -sin(theta)*cos(psi)*cos(phi)-sin(psi)*sin(phi);
+  a2_y = cos(theta)*cos(psi);
   // a2 axis
-  cos_val_a2 =  sin(acos(cos_val_a1));//a2_x*q_x + a2_y*q_y;
-  a2_dot_q = aa*q*cos_val_a2*sin(alpha);
+  cos_val_a2 =  a2_x*q_x + a2_y*q_y;
 
   // The following test should always pass
   if (fabs(cos_val_a3)>1.0) {
-    printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
-    return 0;
-  }
+    //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+    cos_val_a3 = 1.0;
+  }
+   if (fabs(cos_val_a1)>1.0) {
+    //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+    cos_val_a1 = 1.0;
+  }
+   if (fabs(cos_val_a2)>1.0) {
+    //printf("parallel_ana_2D: Unexpected error: cos(alpha)>1\n");
+    cos_val_a3 = 1.0;
+  }
+  a3_dot_q = aa*q*cos_val_a3;
+  a1_dot_q = aa*q*cos_val_a1;//*sin(alpha);
+  a2_dot_q = aa*q*cos_val_a2;
+  
   // Call Zq=Z1*Z2*Z3
   Zq = (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a1_dot_q)+exp(-qDa_2));
-  Zq = Zq * (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a2_dot_q)+exp(-qDa_2));
-  Zq = Zq * (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a3_dot_q)+exp(-qDa_2));
+  Zq *= (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a2_dot_q)+exp(-qDa_2));
+  Zq *= (1.0-exp(-qDa_2))/(1.0-2.0*exp(-0.5*qDa_2)*cos(a3_dot_q)+exp(-qDa_2));
 
   // Use SphereForm directly from libigor
@@ -278 +285 @@
           }
           if (weights_theta.size()>1) {
-            _ptvalue *= fabs(sin(weights_theta[j].value*pi/180.0));
+            _ptvalue *= fabs(cos(weights_theta[j].value*pi/180.0));
           }
           sum += _ptvalue;

sansmodels/src/c_models/sphere.cpp

@@ -28 +28 @@
 extern "C" {
         #include "libSphere.h"
-}
+        #include "libmultifunc/libfunc.h"
+}
+// Convenience parameter structure
+typedef struct {
+    double scale;
+    double radius;
+    double sldSph;
+    double sldSolv;
+    double background;
+    double M0_sld_sph;
+    double M_theta_sph;
+    double M_phi_sph;
+    double M0_sld_solv;
+    double M_theta_solv;
+    double M_phi_solv;
+    double Up_frac_i;
+        double Up_frac_f;
+        double Up_theta;
+} SphereParameters;
 
 SphereModel :: SphereModel() {
@@ -37 +55 @@
         sldSolv   = Parameter(1.0e-6);
         background = Parameter(0.0);
+        M0_sld_sph = Parameter(0.0e-6);
+        M_theta_sph = Parameter(0.0);
+        M_phi_sph = Parameter(0.0); 
+        M0_sld_solv = Parameter(0.0e-6);
+        M_theta_solv = Parameter(0.0);
+        M_phi_solv = Parameter(0.0); 
+        Up_frac_i = Parameter(0.5); 
+        Up_frac_f = Parameter(0.5);
+        Up_theta = Parameter(0.0);
 }
 
@@ -87 +114 @@
 /**
  * Function to evaluate 2D scattering function
+ * @param pars: parameters
+ * @param q: q-value
+ * @param q_x: q_x / q
+ * @param q_y: q_y / q
+ * @return: function value
+ */
+
+static double sphere_analytical_2D_scaled(SphereParameters *pars, double q, double q_x, double q_y) {
+        double dp[5];
+        //convert angle degree to radian
+        dp[0] = 1.0;
+        dp[1] = pars->radius;
+        dp[2] = 0.0;
+        dp[3] = 0.0;
+        dp[4] = 0.0;
+
+    double sldSph = pars->sldSph;
+    double sldSolv = pars->sldSolv;
+        double answer = 0.0;
+        double m_max = pars->M0_sld_sph;
+        double m_max_solv = pars->M0_sld_solv;
+
+        if (m_max < 1.0e-32 && m_max_solv < 1.0e-32){
+                dp[2] = sldSph;
+                dp[3] = sldSolv;
+                answer = SphereForm(dp, q);
+        }
+        else{
+                double contrast = sldSph - sldSolv;
+                double qx = q_x;
+                double qy = q_y;
+                double s_theta = pars->Up_theta;
+                double m_phi = pars->M_phi_sph;
+                double m_theta = pars->M_theta_sph;
+                double m_phi_solv = pars->M_phi_solv;
+                double m_theta_solv = pars->M_theta_solv;
+                double in_spin = pars->Up_frac_i;
+                double out_spin = pars->Up_frac_f;
+                polar_sld p_sld;
+                polar_sld p_sld_solv;
+                p_sld = cal_msld(1, qx, qy, sldSph, m_max, m_theta, m_phi, 
+                                                        in_spin, out_spin, s_theta);
+                p_sld_solv = cal_msld(1, qx, qy, sldSolv, m_max_solv, m_theta_solv, m_phi_solv, 
+                                                                in_spin, out_spin, s_theta);
+                //up_up 
+                if (in_spin > 0.0 && out_spin > 0.0){                    
+                        dp[2] = p_sld.uu;
+                        dp[3] = p_sld_solv.uu;
+                        answer += SphereForm(dp, q);
+                        }
+                //down_down
+                if (in_spin < 1.0 && out_spin < 1.0){
+                        dp[2] = p_sld.dd;
+                        dp[3] = p_sld_solv.dd;
+                        answer += SphereForm(dp, q);
+                        }
+                //up_down
+                if (in_spin > 0.0 && out_spin < 1.0){
+                        dp[2] = p_sld.re_ud;
+                        dp[3] = p_sld_solv.re_ud;
+                        answer += SphereForm(dp, q);    
+                        dp[2] = p_sld.im_ud;
+                        dp[3] = p_sld_solv.im_ud;
+                        answer += SphereForm(dp, q);
+                        }
+                //down_up       
+                if (in_spin < 1.0 && out_spin > 0.0){
+                        dp[2] = p_sld.re_du;
+                        dp[3] = p_sld_solv.re_du;
+                        answer += SphereForm(dp, q);    
+                        dp[2] = p_sld.im_du;
+                        dp[3] = p_sld_solv.im_du;
+                        answer += SphereForm(dp, q);
+                        }
+        }
+        
+        // add in the background
+        answer *= pars->scale;
+        answer += pars->background;
+        return answer;
+}
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters
+ * @param q: q-value
+ * @return: function value
+ */
+static double sphere_analytical_2DXY(SphereParameters *pars, double qx, double qy) {
+  double q;
+  q = sqrt(qx*qx+qy*qy);
+  return sphere_analytical_2D_scaled(pars, q, qx/q, qy/q);
+}
+
+
+/**
+ * Function to evaluate 2D scattering function
  * @param q_x: value of Q along x
  * @param q_y: value of Q along y
@@ -92 +217 @@
  */
 double SphereModel :: operator()(double qx, double qy) {
-        double q = sqrt(qx*qx + qy*qy);
-        return (*this).operator()(q);
+        SphereParameters dp;
+        dp.scale = scale();
+        dp.radius = radius();
+        dp.sldSph = sldSph();
+        dp.sldSolv = sldSolv();
+        dp.background = 0.0;
+        dp.Up_theta =  Up_theta();
+        dp.M_phi_sph =  M_phi_sph();
+        dp.M_theta_sph =  M_theta_sph();
+        dp.M0_sld_sph =  M0_sld_sph();
+        dp.M_phi_solv =  M_phi_solv();
+        dp.M_theta_solv =  M_theta_solv();
+        dp.M0_sld_solv =  M0_sld_solv();
+        dp.Up_frac_i =  Up_frac_i();
+        dp.Up_frac_f =  Up_frac_f();
+
+        // 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;
+
+        // Loop over radius weight points
+        for(size_t i=0; i<weights_rad.size(); i++) {
+                dp.radius = weights_rad[i].value;
+
+                //Un-normalize SphereForm by volume
+                sum += weights_rad[i].weight
+                        * sphere_analytical_2DXY(&dp, qx, qy) * 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;
+        }
+
+        if (vol != 0.0 && norm != 0.0) {
+                //Re-normalize by avg volume
+                sum = sum/(vol/norm);}
+        return sum/norm + background();
 }
 
@@ -104 +270 @@
  */
 double SphereModel :: evaluate_rphi(double q, double phi) {
-        return (*this).operator()(q);
+        double qx = q*cos(phi);
+        double qy = q*sin(phi);
+        return (*this).operator()(qx, qy);
 }
 

sansmodels/src/c_models/spheroid.cpp

@@ -57 +57 @@
 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 cyl_x, cyl_y;//, cyl_z;
+  //double q_z;
   double alpha, vol, cos_val;
   double answer;
@@ -70 +70 @@
 
   // 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);
+  cyl_x = cos(theta) * cos(phi);
+  cyl_y = sin(theta);
+  //cyl_z = -cos(theta) * sin(phi);
   //del sld
   sldcs = pars->sld_core - pars->sld_shell;
@@ -78 +78 @@
 
   // q vector
-  q_z = 0;
+  //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;
+  cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
 
   // The following test should always pass
@@ -336 +336 @@
               * 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));
+                _ptvalue *= fabs(cos(weights_theta[m].value*pi/180.0));
               }
               sum += _ptvalue;

sansmodels/src/c_models/stackeddisks.cpp

@@ -58 +58 @@
  */
 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 cyl_x, cyl_y;//, cyl_z;
+  //double q_z;
   double alpha, vol, cos_val;
   double d, dum, halfheight;
@@ -68 +68 @@
 
 
-
   // parallelepiped 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);
 
   // q vector
-  q_z = 0;
+  //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;
+  cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z;
 
   // The following test should always pass
@@ -265 +263 @@
   double vol = 0.0;
   double pi = 4.0*atan(1.0);
-
+  
   // Loop over length weight points
   for(int i=0; i< (int)weights_core_thick.size(); i++) {
@@ -294 +292 @@
             *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*pi/180.0));
+              _ptvalue *= fabs(cos(weights_theta[l].value*pi/180.0));
             }
             sum += _ptvalue;

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;