Index: explore/angles.py
===================================================================
--- explore/angles.py	(revision ef8e68c4318261011d61de9f7eed793910372509)
+++ explore/angles.py	(revision f728001e783d4f2cad1f52bfa5374c457648c268)
@@ -38,5 +38,6 @@
 
 # include unicode symbols in output, even if piping to a pager
-sys.stdout = codecs.getwriter('utf8')(sys.stdout)
+if sys.version_info[0] < 3:
+    sys.stdout = codecs.getwriter('utf8')(sys.stdout)
 sp.init_printing(use_unicode=True)
 
Index: sasmodels/models/bcc_paracrystal.c
===================================================================
--- sasmodels/models/bcc_paracrystal.c	(revision becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/bcc_paracrystal.c	(revision f728001e783d4f2cad1f52bfa5374c457648c268)
@@ -2,21 +2,19 @@
 bcc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
 {
-#if 0  // Equations as written in Matsuoka
-    const double a1 = (+qa + qb + qc)/2.0;
-    const double a2 = (-qa - qb + qc)/2.0;
-    const double a3 = (-qa + qb - qc)/2.0;
-#else
-    const double a1 = (+qa + qb - qc)/2.0;
+    // Equations from Matsuoka 26-27-28, multiplied by |q|
+    const double a1 = (-qa + qb + qc)/2.0;
     const double a2 = (+qa - qb + qc)/2.0;
-    const double a3 = (-qa + qb + qc)/2.0;
-#endif
+    const double a3 = (+qa + qb - qc)/2.0;
 
-#if 1
-    // Numerator: (1 - exp(a)^2)^3
-    //         => (-(exp(2a) - 1))^3
-    //         => -expm1(2a)^3
-    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a))
-    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1)
-    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1)
+#if 0
+    // Matsuoka 29-30-31
+    //     Z_k numerator: 1 - exp(a)^2
+    //     Z_k denominator: 1 - 2 cos(d a_k) exp(a) + exp(2a)
+    // Rewriting numerator
+    //         => -(exp(2a) - 1)
+    //         => -expm1(2a)
+    // Rewriting denominator
+    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1)
+    //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1
     const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
     const double exp_arg = exp(arg);
@@ -25,7 +23,18 @@
           * ((exp_arg - 2.0*cos(dnn*a2))*exp_arg + 1.0)
           * ((exp_arg - 2.0*cos(dnn*a3))*exp_arg + 1.0));
-#else
-    // Alternate form, which perhaps is more approachable
-    const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+
+#elif 0
+    // ** Alternate form, which perhaps is more approachable
+    //     Z_k numerator   => -[(exp(2a) - 1) / 2.exp(a)] 2.exp(a)
+    //                     => -[sinh(a)] exp(a)
+    //     Z_k denominator => [(exp(2a) + 1) / 2.exp(a) - cos(d a_k)] 2.exp(a)
+    //                     => [cosh(a) - cos(d a_k)] 2.exp(a)
+    //     => Z_k = -sinh(a) / [cosh(a) - cos(d a_k)]
+    //            = sinh(-a) / [cosh(-a) - cos(d a_k)]
+    //
+    // One more step leads to the form in sasview 3.x for 2d models
+    //            = tanh(-a) / [1 - cos(d a_k)/cosh(-a)]
+    //
+    const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
     const double sinh_qd = sinh(arg);
     const double cosh_qd = cosh(arg);
@@ -33,4 +42,11 @@
                     * sinh_qd/(cosh_qd - cos(dnn*a2))
                     * sinh_qd/(cosh_qd - cos(dnn*a3));
+#else
+    const double arg = 0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
+    const double tanh_qd = tanh(arg);
+    const double cosh_qd = cosh(arg);
+    const double Zq = tanh_qd/(1.0 - cos(dnn*a1)/cosh_qd)
+                    * tanh_qd/(1.0 - cos(dnn*a2)/cosh_qd)
+                    * tanh_qd/(1.0 - cos(dnn*a3)/cosh_qd);
 #endif
 
Index: sasmodels/models/fcc_paracrystal.c
===================================================================
--- sasmodels/models/fcc_paracrystal.c	(revision becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/fcc_paracrystal.c	(revision f728001e783d4f2cad1f52bfa5374c457648c268)
@@ -2,20 +2,18 @@
 fcc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
 {
-#if 0  // Equations as written in Matsuoka
-    const double a1 = ( qa + qb)/2.0;
-    const double a2 = (-qa + qc)/2.0;
-    const double a3 = (-qa + qb)/2.0;
-#else
+    // Equations from Matsuoka 17-18-19, multiplied by |q|
     const double a1 = ( qa + qb)/2.0;
     const double a2 = ( qa + qc)/2.0;
     const double a3 = ( qb + qc)/2.0;
-#endif
 
-    // Numerator: (1 - exp(a)^2)^3
-    //         => (-(exp(2a) - 1))^3
-    //         => -expm1(2a)^3
-    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a))
-    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1)
-    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1)
+    // Matsuoka 23-24-25
+    //     Z_k numerator: 1 - exp(a)^2
+    //     Z_k denominator: 1 - 2 cos(d a_k) exp(a) + exp(2a)
+    // Rewriting numerator
+    //         => -(exp(2a) - 1)
+    //         => -expm1(2a)
+    // Rewriting denominator
+    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1)
+    //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1
     const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
     const double exp_arg = exp(arg);
Index: sasmodels/models/sc_paracrystal.c
===================================================================
--- sasmodels/models/sc_paracrystal.c	(revision becded3d803d05553cd978e361eeb35b6f6736bf)
+++ sasmodels/models/sc_paracrystal.c	(revision f728001e783d4f2cad1f52bfa5374c457648c268)
@@ -2,14 +2,18 @@
 sc_Zq(double qa, double qb, double qc, double dnn, double d_factor)
 {
+    // Equations from Matsuoka 9-10-11, multiplied by |q|
     const double a1 = qa;
     const double a2 = qb;
     const double a3 = qc;
 
-    // Numerator: (1 - exp(a)^2)^3
-    //         => (-(exp(2a) - 1))^3
-    //         => -expm1(2a)^3
-    // Denominator: prod(1 - 2 cos(d ak) exp(a) + exp(2a))
-    //         => prod(exp(a)^2 - 2 cos(d ak) exp(a) + 1)
-    //         => prod((exp(a) - 2 cos(d ak)) * exp(a) + 1)
+    // Matsuoka 13-14-15
+    //     Z_k numerator: 1 - exp(a)^2
+    //     Z_k denominator: 1 - 2 cos(d a_k) exp(a) + exp(2a)
+    // Rewriting numerator
+    //         => -(exp(2a) - 1)
+    //         => -expm1(2a)
+    // Rewriting denominator
+    //         => exp(a)^2 - 2 cos(d ak) exp(a) + 1)
+    //         => (exp(a) - 2 cos(d ak)) * exp(a) + 1
     const double arg = -0.5*square(dnn*d_factor)*(a1*a1 + a2*a2 + a3*a3);
     const double exp_arg = exp(arg);
