Index: sasmodels/kernel_header.c
===================================================================
--- sasmodels/kernel_header.c	(revision e1d6983831831051b424a8d021b337435d87f0f3)
+++ sasmodels/kernel_header.c	(revision 11ca2ab20045d2ca8709faf56c10caaeaaa02734)
@@ -148,2 +148,23 @@
 inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; }
 
+ #define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \
+     SINCOS(phi*M_PI_180, sn, cn); \
+     q = sqrt(qx*qx + qy*qy); \
+     cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * cos(theta*M_PI_180));  \
+     sn = sqrt(1 - cn*cn); \
+ } while (0)
+
+ #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \
+    q = sqrt(qx*qx + qy*qy); \
+    const double qxhat = qx/q; \
+    const double qyhat = qy/q; \
+    double sin_theta, cos_theta; \
+    double sin_phi, cos_phi; \
+    double sin_psi, cos_psi; \
+    SINCOS(theta*M_PI_180, sin_theta, cos_theta); \
+    SINCOS(phi*M_PI_180, sin_phi, cos_phi); \
+    SINCOS(psi*M_PI_180, sin_psi, cos_psi); \
+    cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; \
+    cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; \
+    cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \
+ } while (0)
Index: sasmodels/models/barbell.c
===================================================================
--- sasmodels/models/barbell.c	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
+++ sasmodels/models/barbell.c	(revision 11ca2ab20045d2ca8709faf56c10caaeaaa02734)
@@ -39,7 +39,21 @@
     // translate dx in [-1,1] to dx in [lower,upper]
     const double integral = total*zm;
-    const double bell_Fq = 2*M_PI*cube(radius_bell)*integral;
-    return bell_Fq;
+    const double bell_fq = 2*M_PI*cube(radius_bell)*integral;
+    return bell_fq;
 }
+
+static double
+_fq(double q, double h,
+    double radius_bell, double radius, double half_length,
+    double sin_alpha, double cos_alpha)
+{
+    const double bell_fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha);
+    const double bj = sas_J1c(q*radius*sin_alpha);
+    const double si = sinc(q*half_length*cos_alpha);
+    const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si;
+    const double Aq = bell_fq + cyl_fq;
+    return Aq;
+}
+
 
 double form_volume(double radius_bell,
@@ -47,5 +61,4 @@
         double length)
 {
-
     // bell radius should never be less than radius when this is called
     const double hdist = sqrt(square(radius_bell) - square(radius));
@@ -71,10 +84,5 @@
         double sin_alpha, cos_alpha; // slots to hold sincos function output
         SINCOS(alpha, sin_alpha, cos_alpha);
-
-        const double bell_Fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha);
-        const double bj = sas_J1c(q*radius*sin_alpha);
-        const double si = sinc(q*half_length*cos_alpha);
-        const double cyl_Fq = M_PI*radius*radius*length*bj*si;
-        const double Aq = bell_Fq + cyl_Fq;
+        const double Aq = _fq(q, h, radius_bell, radius, half_length, sin_alpha, cos_alpha);
         total += Gauss76Wt[i] * Aq * Aq * sin_alpha;
     }
@@ -93,21 +101,9 @@
         double theta, double phi)
 {
-    // Compute angle alpha between q and the cylinder axis
-    double sn, cn; // slots to hold sincos function output
-    SINCOS(phi*M_PI_180, sn, cn);
-    const double q = sqrt(qx*qx + qy*qy);
-    const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q);
-    const double alpha = acos(cos_val); // rod angle relative to q
+    double q, sin_alpha, cos_alpha;
+    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
 
     const double h = -sqrt(square(radius_bell) - square(radius));
-    const double half_length = 0.5*length;
-
-    double sin_alpha, cos_alpha; // slots to hold sincos function output
-    SINCOS(alpha, sin_alpha, cos_alpha);
-    const double bell_Fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha);
-    const double bj = sas_J1c(q*radius*sin_alpha);
-    const double si = sinc(q*half_length*cos_alpha);
-    const double cyl_Fq = M_PI*radius*radius*length*bj*si;
-    const double Aq = cyl_Fq + bell_Fq;
+    const double Aq = _fq(q, h, radius_bell, radius, 0.5*length, sin_alpha, cos_alpha);
 
     // Multiply by contrast^2 and convert to cm-1
Index: sasmodels/models/bcc_paracrystal.c
===================================================================
--- sasmodels/models/bcc_paracrystal.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
+++ sasmodels/models/bcc_paracrystal.c	(revision 11ca2ab20045d2ca8709faf56c10caaeaaa02734)
@@ -31,4 +31,5 @@
 	const double temp8 = -sin_theta*cos_phi - sin_theta*sin_phi + cos_theta;
 	const double temp9 = -sin_theta*cos_phi + sin_theta*sin_phi - cos_theta;
+
 	const double temp10 = exp((-1.0/8.0)*temp1*(temp7*temp7 + temp8*temp8 + temp9*temp9));
 	result = pow(1.0-(temp10*temp10),3)*temp6
@@ -81,90 +82,31 @@
 
     return answer;
-
-
 }
 
 
-double Iqxy(double qx, double qy, double dnn,
-    double d_factor, double radius,double sld, double solvent_sld,
-    double theta, double phi, double psi){
+double Iqxy(double qx, double qy,
+    double dnn, double d_factor, double radius,
+    double sld, double solvent_sld,
+    double theta, double phi, double psi)
+{
+    double q, cos_a1, cos_a2, cos_a3;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1);
 
-  double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
-  //double q_z;
-  double cos_val_b3, cos_val_b2, cos_val_b1;
-  double a1_dot_q, a2_dot_q,a3_dot_q;
-  double answer;
-  double Zq, Fkq, Fkq_2;
+    const double a1 = +cos_a3 - cos_a1 + cos_a2;
+    const double a2 = +cos_a3 + cos_a1 - cos_a2;
+    const double a3 = -cos_a3 + cos_a1 + cos_a2;
 
-  //convert to q and make scaled values
-  double q = sqrt(qx*qx+qy*qy);
-  double q_x = qx/q;
-  double q_y = qy/q;
+    const double qd = 0.5*q*dnn;
+    const double exp_qd = exp(0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3));
+    const double sinh_qd = 0.5*exp_qd - 0.5/exp_qd;
+    const double cosh_qd = 0.5*exp_qd + 0.5/exp_qd;
 
-  //convert angle degree to radian
-  theta = theta * M_PI_180;
-  phi = phi * M_PI_180;
-  psi = psi * M_PI_180;
+    const double Zq = sinh_qd/(cosh_qd - cos(qd*a1))
+                    * sinh_qd/(cosh_qd - cos(qd*a2))
+                    * sinh_qd/(cosh_qd - cos(qd*a3));
 
-  const double Da = d_factor*dnn;
-  const double s1 = dnn/sqrt(0.75);
-
-
-  //the occupied volume of the lattice
-  const double latticescale = 2.0*sphere_volume(radius/s1);
-  // 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 = 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 = -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
-    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*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3);
-
-    // a1 axis
-    a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1);
-
-    // a2 axis
-    a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2);
-
-
-    // Get Fkq and Fkq_2
-    Fkq = exp(-0.5*pow(Da/dnn,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 = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale;
-
-  return answer;
- }
+    const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq;
+    //the occupied volume of the lattice
+    const double lattice_scale = 2.0*sphere_volume(sqrt(0.75)*radius/dnn);
+    return lattice_scale * Fq;
+}
Index: sasmodels/models/cylinder.c
===================================================================
--- sasmodels/models/cylinder.c	(revision 0d6e86566d51a1fba22933bdbac8c40d36f9bd76)
+++ sasmodels/models/cylinder.c	(revision 11ca2ab20045d2ca8709faf56c10caaeaaa02734)
@@ -33,5 +33,5 @@
         // alpha(theta,phi) the projection of the cylinder on the detector plane
         SINCOS(alpha, sn, cn);
-        total += Gauss76Wt[i] * fq(q, sn, cn, radius, length)* fq(q, sn, cn, radius, length) * sn;
+        total += Gauss76Wt[i] * square(fq(q, sn, cn, radius, length)) * sn;
     }
     // translate dx in [-1,1] to dx in [lower,upper]
@@ -58,15 +58,7 @@
     double phi)
 {
-    double sn, cn; // slots to hold sincos function output
-
-    // Compute angle alpha between q and the cylinder axis
-    SINCOS(phi*M_PI_180, sn, cn);
-    const double q = sqrt(qx*qx + qy*qy);
-    const double cos_val = (q==0. ? 1.0 : (cn*qx + sn*qy)*sin(theta*M_PI_180)/q);
-
-    const double alpha = acos(cos_val);
-
-    SINCOS(alpha, sn, cn);
+    double q, sin_alpha, cos_alpha;
+    ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha);
     const double s = (sld-solvent_sld) * form_volume(radius, length);
-    return 1.0e-4 * square(s * fq(q, sn, cn, radius, length));
+    return 1.0e-4 * square(s * fq(q, sin_alpha, cos_alpha, radius, length));
 }
Index: sasmodels/models/fcc_paracrystal.c
===================================================================
--- sasmodels/models/fcc_paracrystal.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
+++ sasmodels/models/fcc_paracrystal.c	(revision 11ca2ab20045d2ca8709faf56c10caaeaaa02734)
@@ -85,86 +85,27 @@
 }
 
+double Iqxy(double qx, double qy,
+    double dnn, double d_factor, double radius,
+    double sld, double solvent_sld,
+    double theta, double phi, double psi)
+{
+    double q, cos_a1, cos_a2, cos_a3;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1);
 
-double Iqxy(double qx, double qy, double dnn,
-    double d_factor, double radius,double sld, double solvent_sld,
-    double theta, double phi, double psi){
+    const double a1 = +cos_a3 + cos_a2;
+    const double a2 = +cos_a3 + cos_a1;
+    const double a3 = -cos_a3 + cos_a2;
+    const double qd = 0.5*q*dnn;
+    const double exp_qd = exp(0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3));
+    const double sinh_qd = 0.5*exp_qd - 0.5/exp_qd;
+    const double cosh_qd = 0.5*exp_qd + 0.5/exp_qd;
 
-  double b3_x, b3_y, b1_x, b1_y, b2_x, b2_y; //b3_z,
-  // double q_z;
-  double cos_val_b3, cos_val_b2, cos_val_b1;
-  double a1_dot_q, a2_dot_q,a3_dot_q;
-  double answer;
-  double Zq, Fkq, Fkq_2;
+    const double Zq = sinh_qd/(cosh_qd - cos(qd*a1))
+                    * sinh_qd/(cosh_qd - cos(qd*a2))
+                    * sinh_qd/(cosh_qd - cos(qd*a3));
 
-  //convert to q and make scaled values
-  double q = sqrt(qx*qx+qy*qy);
-  double q_x = qx/q;
-  double q_y = qy/q;
-
-  //convert angle degree to radian
-  theta = theta * M_PI_180;
-  phi = phi * M_PI_180;
-  psi = psi * M_PI_180;
-
-  const double Da = d_factor*dnn;
-  const double s1 = dnn/sqrt(0.75);
-
-
-  //the occupied volume of the lattice
-  const double latticescale = 2.0*sphere_volume(radius/s1);
-  // 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 = 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 = -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
-    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("FCC_ana_2D: Unexpected error: cos()>1\n");
-      cos_val_b3 = 1.0;
-    }
-    if (fabs(cos_val_b2)>1.0) {
-      //printf("FCC_ana_2D: Unexpected error: cos()>1\n");
-      cos_val_b2 = 1.0;
-    }
-    if (fabs(cos_val_b1)>1.0) {
-      //printf("FCC_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*dnn*q*(cos_val_b2+cos_val_b1-cos_val_b3);
-
-    // a1 axis
-    a1_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b2-cos_val_b1);
-
-    // a2 axis
-    a2_dot_q = 0.5*dnn*q*(cos_val_b3+cos_val_b1-cos_val_b2);
-
-
-    // Get Fkq and Fkq_2
-    Fkq = exp(-0.5*pow(Da/dnn,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 = sphere_form(q,radius,sld,solvent_sld)*Zq*latticescale;
-
-  return answer;
- }
+    const double Fq = sphere_form(q,radius,sld,solvent_sld)*Zq;
+    //the occupied volume of the lattice
+    const double lattice_scale = 2.0*sphere_volume(M_SQRT1_2*radius/dnn);
+    return lattice_scale * Fq;
+}
Index: sasmodels/models/sc_paracrystal.c
===================================================================
--- sasmodels/models/sc_paracrystal.c	(revision 0bef47b278e819df7ea6b572375dd7bc9607227e)
+++ sasmodels/models/sc_paracrystal.c	(revision 11ca2ab20045d2ca8709faf56c10caaeaaa02734)
@@ -60,6 +60,5 @@
 }
 
-static
-double sc_crystal_kernel(double q,
+double Iq(double q,
           double dnn,
           double d_factor,
@@ -103,6 +102,5 @@
 }
 
-static
-double sc_crystal_kernel_2d(double q, double q_x, double q_y,
+double Iqxy(double qx, double qy,
           double dnn,
           double d_factor,
@@ -114,115 +112,19 @@
           double psi)
 {
-    //convert angle degree to radian
-    theta = theta * M_PI_180;
-    phi = phi * M_PI_180;
-    psi = psi * M_PI_180;
+    double q, cos_a1, cos_a2, cos_a3;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1);
 
-    const double qda_2 = pow(q*d_factor*dnn,2.0);
+    const double qd = q*dnn;
+    const double exp_qd = exp(0.5*square(qd*d_factor));
+    const double sinh_qd = 0.5*exp_qd - 0.5/exp_qd;
+    const double cosh_qd = 0.5*exp_qd + 0.5/exp_qd;
 
-    double snt, cnt;
-    SINCOS(theta, snt, cnt);
+    const double Zq = sinh_qd/(cosh_qd - cos(qd*cos_a1))
+                    * sinh_qd/(cosh_qd - cos(qd*cos_a2))
+                    * sinh_qd/(cosh_qd - cos(qd*cos_a3));
 
-    double snp, cnp;
-    SINCOS(phi, snp, cnp);
-
-    double sns, cns;
-    SINCOS(psi, sns, cns);
-
-    /// Angles here are respect to detector coordinate instead of against
-    //  q coordinate(PRB 36, 3, 1754)
-    // a3 axis orientation
-
-    const double a3_x = cnt * cnp;
-    const double a3_y = snt;
-
-    // Compute the angle btw vector q and the a3 axis
-    double cos_val_a3 = a3_x*q_x + a3_y*q_y;
-
-    // a1 axis orientation
-    const double a1_x = -cnp*sns * snt+snp*cns;
-    const double a1_y = sns*cnt;
-
-    double cos_val_a1 = a1_x*q_x + a1_y*q_y;
-
-    // a2 axis orientation
-    const double a2_x = -snt*cns*cnp-sns*snp;
-    const double a2_y = cnt*cns;
-
-    // a2 axis
-    const double 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");
-        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;
-    }
-
-    const double a3_dot_q = dnn*q*cos_val_a3;
-    const double a1_dot_q = dnn*q*cos_val_a1;
-    const double a2_dot_q = dnn*q*cos_val_a2;
-
-    // Call Zq=Z1*Z2*Z3
-    double Zq = (1.0-exp(-qda_2))/(1.0-2.0*exp(-0.5*qda_2)*cos(a1_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
-    double answer = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq;
-
-    //consider scales
-    const double latticeScale = sphere_volume(radius/dnn);
-    answer *= latticeScale;
-
-    return answer;
+    const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq;
+    //the occupied volume of the lattice
+    const double lattice_scale = sphere_volume(radius/dnn);
+    return lattice_scale * Fq;
 }
-
-double Iq(double q,
-          double dnn,
-          double d_factor,
-          double radius,
-          double sphere_sld,
-          double solvent_sld)
-{
-    return sc_crystal_kernel(q,
-              dnn,
-              d_factor,
-              radius,
-              sphere_sld,
-              solvent_sld);
-}
-
-// Iqxy is never called since no orientation or magnetic parameters.
-double Iqxy(double qx, double qy,
-            double dnn,
-            double d_factor,
-            double radius,
-            double sphere_sld,
-            double solvent_sld,
-            double theta,
-            double phi,
-            double psi)
-{
-    double q = sqrt(qx*qx + qy*qy);
-
-
-    return sc_crystal_kernel_2d(q, qx/q, qy/q,
-                  dnn,
-                  d_factor,
-                  radius,
-                  sphere_sld,
-                  solvent_sld,
-                  theta,
-                  phi,
-                  psi);
-
-}
-
Index: sasmodels/models/triaxial_ellipsoid.c
===================================================================
--- sasmodels/models/triaxial_ellipsoid.c	(revision a8072067ae8a224e45d23ba71f663aadbe1ba4a9)
+++ sasmodels/models/triaxial_ellipsoid.c	(revision 11ca2ab20045d2ca8709faf56c10caaeaaa02734)
@@ -58,17 +58,7 @@
     double psi)
 {
-    double stheta, ctheta;
-    double sphi, cphi;
-    double spsi, cpsi;
+    double q, calpha, cmu, cnu;
+    ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, calpha, cmu, cnu);
 
-    const double q = sqrt(qx*qx + qy*qy);
-    const double qxhat = qx/q;
-    const double qyhat = qy/q;
-    SINCOS(theta*M_PI_180, stheta, ctheta);
-    SINCOS(phi*M_PI_180, sphi, cphi);
-    SINCOS(psi*M_PI_180, spsi, cpsi);
-    const double calpha = ctheta*cphi*qxhat + stheta*qyhat;
-    const double cnu = (-cphi*spsi*stheta + sphi*cpsi)*qxhat + spsi*ctheta*qyhat;
-    const double cmu = (-stheta*cpsi*cphi - spsi*sphi)*qxhat + ctheta*cpsi*qyhat;
     const double t = q*sqrt(radius_equat_minor*radius_equat_minor*cnu*cnu
                           + radius_equat_major*radius_equat_major*cmu*cmu
