Index: sasmodels/compare.py
===================================================================
--- sasmodels/compare.py	(revision 8b25ee110ba3c3b5c878ec0f1b620daff4a34ef1)
+++ sasmodels/compare.py	(revision 82c299fa51ab091ee8ead0a235d4177ccfeecb99)
@@ -123,4 +123,15 @@
         rg_max = np.sqrt(90*np.log(10) + 3*np.log(pars['scale']))/q_max
         pars['rg'] = min(pars['rg'],rg_max)
+
+    if name == 'rpa':
+        # Make sure phi sums to 1.0
+        if pars['case_num'] < 2:
+            pars['Phia'] = 0.
+            pars['Phib'] = 0.
+        elif pars['case_num'] < 5:
+            pars['Phia'] = 0.
+        total = sum(pars['Phi'+c] for c in 'abcd')
+        for c in 'abcd':
+            pars['Phi'+c] /= total
 
 def parlist(pars):
Index: sasmodels/convert.py
===================================================================
--- sasmodels/convert.py	(revision 55a83169d3093c778de6a1637b5c8cb54039e3a0)
+++ sasmodels/convert.py	(revision 55a83169d3093c778de6a1637b5c8cb54039e3a0)
@@ -112,4 +112,8 @@
         _remove_pd(oldpars, 'num_pearls', name)
         _remove_pd(oldpars, 'thick_string', name)
+    elif name == 'rpa':
+        # convert scattering lengths from femtometers to centimeters
+        for p in "La", "Lb", "Lc", "Ld":
+            if p in oldpars: oldpars[p] *= 1e-13
 
     return oldname, oldpars
@@ -128,4 +132,6 @@
         pars['string_thickness_pd_n'] = 0
         pars['number_of_pearls_pd_n'] = 0
+    elif name == 'rpa':
+        pars['case_num'] = int(pars['case_num'])
     elif getattr(model_definition, 'category', None) == 'structure-factor':
         pars['scale'], pars['background'] = 1, 0
Index: sasmodels/generate.py
===================================================================
--- sasmodels/generate.py	(revision e97170c97d135b16697eb2bd7399dfa83eec6cf6)
+++ sasmodels/generate.py	(revision e66c9f9bddfbf0bb75d5d64484195eaf6802ca82)
@@ -265,6 +265,6 @@
 """
 
-def format_units(par):
-    return RST_UNITS.get(par, par)
+def format_units(units):
+    return "string" if isinstance(units, list) else RST_UNITS.get(units, units)
 
 def make_partable(pars):
@@ -284,5 +284,6 @@
     lines = [
         sep,
-        " ".join("%-*s" % (w, h) for w, h in zip(column_widths, PARTABLE_HEADERS)),
+        " ".join("%-*s" % (w, h)
+                 for w, h in zip(column_widths, PARTABLE_HEADERS)),
         sep,
         ]
Index: sasmodels/models/rpa.c
===================================================================
--- sasmodels/models/rpa.c	(revision 82c299fa51ab091ee8ead0a235d4177ccfeecb99)
+++ sasmodels/models/rpa.c	(revision 82c299fa51ab091ee8ead0a235d4177ccfeecb99)
@@ -0,0 +1,343 @@
+double Iq(double q, double case_num,
+    double Na, double Phia, double va, double a_sld, double ba,
+    double Nb, double Phib, double vb, double b_sld, double bb,
+    double Nc, double Phic, double vc, double c_sld, double bc,
+    double Nd, double Phid, double vd, double d_sld, double bd,
+    double Kab, double Kac, double Kad,
+    double Kbc, double Kbd, double Kcd
+    );
+
+double Iqxy(double qx, double qy, double case_num,
+    double Na, double Phia, double va, double a_sld, double ba,
+    double Nb, double Phib, double vb, double b_sld, double bb,
+    double Nc, double Phic, double vc, double c_sld, double bc,
+    double Nd, double Phid, double vd, double d_sld, double bd,
+    double Kab, double Kac, double Kad,
+    double Kbc, double Kbd, double Kcd
+    );
+
+double form_volume(void);
+
+double form_volume(void)
+{
+    return 1.0;
+}
+
+double Iq(double q, double case_num,
+    double Na, double Phia, double va, double La, double ba,
+    double Nb, double Phib, double vb, double Lb, double bb,
+    double Nc, double Phic, double vc, double Lc, double bc,
+    double Nd, double Phid, double vd, double Ld, double bd,
+    double Kab, double Kac, double Kad,
+    double Kbc, double Kbd, double Kcd
+    ) {
+  int icase = (int)case_num;
+
+#if 0  // Sasview defaults
+  if (icase <= 1) {
+    Na=Nb=1000.0;
+    Phia=Phib=0.0000001;
+    Kab=Kac=Kad=Kbc=Kbd=-0.0004;
+    La=Lb=1e-12;
+    va=vb=100.0;
+    ba=bb=5.0;
+  } else if (icase <= 4) {
+    Phia=0.0000001;
+    Kab=Kac=Kad=-0.0004;
+    La=1e-12;
+    va=100.0;
+    ba=5.0;
+  }
+#else
+  if (icase <= 1) {
+    Na=Nb=0.0;
+    Phia=Phib=0.0;
+    Kab=Kac=Kad=Kbc=Kbd=0.0;
+    La=Lb=Ld;
+    va=vb=vd;
+    ba=bb=0.0;
+  } else if (icase <= 4) {
+    Na = 0.0;
+    Phia=0.0;
+    Kab=Kac=Kad=0.0;
+    La=Ld;
+    va=vd;
+    ba=0.0;
+  }
+#endif
+
+  const double Xa = q*q*ba*ba*Na/6.0;
+  const double Xb = q*q*bb*bb*Nb/6.0;
+  const double Xc = q*q*bc*bc*Nc/6.0;
+  const double Xd = q*q*bd*bd*Nd/6.0;
+
+  // limit as Xa goes to 0 is 1
+  const double Pa = Xa==0 ? 1.0 : -expm1(-Xa)/Xa;
+  const double Pb = Xb==0 ? 1.0 : -expm1(-Xb)/Xb;
+  const double Pc = Xc==0 ? 1.0 : -expm1(-Xc)/Xc;
+  const double Pd = Xd==0 ? 1.0 : -expm1(-Xd)/Xd;
+
+  // limit as Xa goes to 0 is 1
+  const double Paa = Xa==0 ? 1.0 : 2.0*(1.0-Pa)/Xa;
+  const double Pbb = Xb==0 ? 1.0 : 2.0*(1.0-Pb)/Xb;
+  const double Pcc = Xc==0 ? 1.0 : 2.0*(1.0-Pc)/Xc;
+  const double Pdd = Xd==0 ? 1.0 : 2.0*(1.0-Pd)/Xd;
+
+
+  // Note: S0ij only defined for copolymers; otherwise set to zero
+  // 0: C/D     binary mixture
+  // 1: C-D     diblock copolymer
+  // 2: B/C/D   ternery mixture
+  // 3: B/C-D   binary mixture,1 homopolymer, 1 diblock copolymer
+  // 4: B-C-D   triblock copolymer
+  // 5: A/B/C/D quaternary mixture
+  // 6: A/B/C-D ternery mixture, 2 homopolymer, 1 diblock copolymer
+  // 7: A/B-C-D binary mixture, 1 homopolymer, 1 triblock copolymer
+  // 8: A-B/C-D binary mixture, 2 diblock copolymer
+  // 9: A-B-C-D tetra-block copolymer
+#if 0
+  const double S0aa = icase<5
+                      ? 1.0 : Na*Phia*va*Paa;
+  const double S0bb = icase<2
+                      ? 1.0 : Nb*Phib*vb*Pbb;
+  const double S0cc = Nc*Phic*vc*Pcc;
+  const double S0dd = Nd*Phid*vd*Pdd;
+  const double S0ab = icase<8
+                      ? 0.0 : sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb;
+  const double S0ac = icase<9
+                      ? 0.0 : sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);
+  const double S0ad = icase<9
+                      ? 0.0 : sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);
+  const double S0bc = (icase!=4 && icase!=7 && icase!= 9)
+                      ? 0.0 : sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc;
+  const double S0bd = (icase!=4 && icase!=7 && icase!= 9)
+                      ? 0.0 : sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);
+  const double S0cd = (icase==0 || icase==2 || icase==5)
+                      ? 0.0 : sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd;
+#else  // sasview equivalent
+//printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,Nc,Phic,vc,Pcc);
+  double S0aa = Na*Phia*va*Paa;
+  double S0bb = Nb*Phib*vb*Pbb;
+  double S0cc = Nc*Phic*vc*Pcc;
+  double S0dd = Nd*Phid*vd*Pdd;
+  double S0ab = sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb;
+  double S0ac = sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb);
+  double S0ad = sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc);
+  double S0bc = sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc;
+  double S0bd = sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc);
+  double S0cd = sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd;
+switch(icase){
+  case 0:
+    S0aa=0.000001;
+    S0ab=0.000002;
+    S0ac=0.000003;
+    S0ad=0.000004;
+    S0bb=0.000005;
+    S0bc=0.000006;
+    S0bd=0.000007;
+    S0cd=0.000008;
+    break;
+  case 1:
+    S0aa=0.000001;
+    S0ab=0.000002;
+    S0ac=0.000003;
+    S0ad=0.000004;
+    S0bb=0.000005;
+    S0bc=0.000006;
+    S0bd=0.000007;
+    break;
+  case 2:
+    S0aa=0.000001;
+    S0ab=0.000002;
+    S0ac=0.000003;
+    S0ad=0.000004;
+    S0bc=0.000005;
+    S0bd=0.000006;
+    S0cd=0.000007;
+    break;
+  case 3:
+    S0aa=0.000001;
+    S0ab=0.000002;
+    S0ac=0.000003;
+    S0ad=0.000004;
+    S0bc=0.000005;
+    S0bd=0.000006;
+    break;
+  case 4:
+    S0aa=0.000001;
+    S0ab=0.000002;
+    S0ac=0.000003;
+    S0ad=0.000004;
+    break;
+  case 5:
+    S0ab=0.000001;
+    S0ac=0.000002;
+    S0ad=0.000003;
+    S0bc=0.000004;
+    S0bd=0.000005;
+    S0cd=0.000006;
+    break;
+  case 6:
+    S0ab=0.000001;
+    S0ac=0.000002;
+    S0ad=0.000003;
+    S0bc=0.000004;
+    S0bd=0.000005;
+    break;
+  case 7:
+    S0ab=0.000001;
+    S0ac=0.000002;
+    S0ad=0.000003;
+    break;
+  case 8:
+    S0ac=0.000001;
+    S0ad=0.000002;
+    S0bc=0.000003;
+    S0bd=0.000004;
+    break;
+  default : //case 9:
+    break;
+  }
+#endif
+
+  // eq 12a: \kappa_{ij}^F = \chi_{ij}^F - \chi_{i0}^F - \chi_{j0}^F
+  const double Kaa = 0.0;
+  const double Kbb = 0.0;
+  const double Kcc = 0.0;
+  const double Kdd = 0.0;
+  const double Zaa = Kaa - Kad - Kad;
+  const double Zab = Kab - Kad - Kbd;
+  const double Zac = Kac - Kad - Kcd;
+  const double Zbb = Kbb - Kbd - Kbd;
+  const double Zbc = Kbc - Kbd - Kcd;
+  const double Zcc = Kcc - Kcd - Kcd;
+//printf("Za: %10.5g %10.5g %10.5g\n", Zaa, Zab, Zac);
+//printf("Zb: %10.5g %10.5g %10.5g\n", Zab, Zbb, Zbc);
+//printf("Zc: %10.5g %10.5g %10.5g\n", Zac, Zbc, Zcc);
+
+  // T = inv(S0)
+  const double DenT = (- S0ac*S0bb*S0ac + S0ab*S0bc*S0ac + S0ac*S0ab*S0bc
+                       - S0aa*S0bc*S0bc - S0ab*S0ab*S0cc + S0aa*S0bb*S0cc);
+  const double T11 = (-S0bc*S0bc + S0bb*S0cc)/DenT;
+  const double T12 = ( S0ac*S0bc - S0ab*S0cc)/DenT;
+  const double T13 = (-S0ac*S0bb + S0ab*S0bc)/DenT;
+  const double T22 = (-S0ac*S0ac + S0aa*S0cc)/DenT;
+  const double T23 = ( S0ac*S0ab - S0aa*S0bc)/DenT;
+  const double T33 = (-S0ab*S0ab + S0aa*S0bb)/DenT;
+
+//printf("T1: %10.5g %10.5g %10.5g\n", T11, T12, T13);
+//printf("T2: %10.5g %10.5g %10.5g\n", T12, T22, T23);
+//printf("T3: %10.5g %10.5g %10.5g\n", T13, T23, T33);
+
+  // eq 18e: m = 1/(S0_{dd} - s0^T inv(S0) s0)
+  const double ZZ = S0ad*(T11*S0ad + T12*S0bd + T13*S0cd)
+                  + S0bd*(T12*S0ad + T22*S0bd + T23*S0cd)
+                  + S0cd*(T13*S0ad + T23*S0bd + T33*S0cd);
+
+  const double m=1.0/(S0dd-ZZ);
+
+  // eq 18d: Y = inv(S0)s0 + e
+  const double Y1 = T11*S0ad + T12*S0bd + T13*S0cd + 1.0;
+  const double Y2 = T12*S0ad + T22*S0bd + T23*S0cd + 1.0;
+  const double Y3 = T13*S0ad + T23*S0bd + T33*S0cd + 1.0;
+
+  // N = mYY^T + \kappa^F
+  const double N11 = m*Y1*Y1 + Zaa;
+  const double N12 = m*Y1*Y2 + Zab;
+  const double N13 = m*Y1*Y3 + Zac;
+  const double N22 = m*Y2*Y2 + Zbb;
+  const double N23 = m*Y2*Y3 + Zbc;
+  const double N33 = m*Y3*Y3 + Zcc;
+
+//printf("N1: %10.5g %10.5g %10.5g\n", N11, N12, N13);
+//printf("N2: %10.5g %10.5g %10.5g\n", N12, N22, N23);
+//printf("N3: %10.5g %10.5g %10.5g\n", N13, N23, N33);
+//printf("S0a: %10.5g %10.5g %10.5g\n", S0aa, S0ab, S0ac);
+//printf("S0b: %10.5g %10.5g %10.5g\n", S0ab, S0bb, S0bc);
+//printf("S0c: %10.5g %10.5g %10.5g\n", S0ac, S0bc, S0cc);
+
+  // M = I + S0 N
+  const double Maa = N11*S0aa + N12*S0ab + N13*S0ac + 1.0;
+  const double Mab = N11*S0ab + N12*S0bb + N13*S0bc;
+  const double Mac = N11*S0ac + N12*S0bc + N13*S0cc;
+  const double Mba = N12*S0aa + N22*S0ab + N23*S0ac;
+  const double Mbb = N12*S0ab + N22*S0bb + N23*S0bc + 1.0;
+  const double Mbc = N12*S0ac + N22*S0bc + N23*S0cc;
+  const double Mca = N13*S0aa + N23*S0ab + N33*S0ac;
+  const double Mcb = N13*S0ab + N23*S0bb + N33*S0bc;
+  const double Mcc = N13*S0ac + N23*S0bc + N33*S0cc + 1.0;
+//printf("M1: %10.5g %10.5g %10.5g\n", Maa, Mab, Mac);
+//printf("M2: %10.5g %10.5g %10.5g\n", Mba, Mbb, Mbc);
+//printf("M3: %10.5g %10.5g %10.5g\n", Mca, Mcb, Mcc);
+
+  // Q = inv(M) = inv(I + S0 N)
+  const double DenQ = (+ Maa*Mbb*Mcc - Maa*Mbc*Mcb - Mab*Mba*Mcc
+                       + Mab*Mbc*Mca + Mac*Mba*Mcb - Mac*Mbb*Mca);
+
+  const double Q11 = ( Mbb*Mcc - Mbc*Mcb)/DenQ;
+  const double Q12 = (-Mab*Mcc + Mac*Mcb)/DenQ;
+  const double Q13 = ( Mab*Mbc - Mac*Mbb)/DenQ;
+  const double Q21 = (-Mba*Mcc + Mbc*Mca)/DenQ;
+  const double Q22 = ( Maa*Mcc - Mac*Mca)/DenQ;
+  const double Q23 = (-Maa*Mbc + Mac*Mba)/DenQ;
+  const double Q31 = ( Mba*Mcb - Mbb*Mca)/DenQ;
+  const double Q32 = (-Maa*Mcb + Mab*Mca)/DenQ;
+  const double Q33 = ( Maa*Mbb - Mab*Mba)/DenQ;
+
+//printf("Q1: %10.5g %10.5g %10.5g\n", Q11, Q12, Q13);
+//printf("Q2: %10.5g %10.5g %10.5g\n", Q21, Q22, Q23);
+//printf("Q3: %10.5g %10.5g %10.5g\n", Q31, Q32, Q33);
+  // eq 18c: inv(S) = inv(S0) + mYY^T + \kappa^F
+  // eq A1 in the appendix
+  // To solve for S, use:
+  //      S = inv(inv(S^0) + N) inv(S^0) S^0
+  //        = inv(S^0 inv(S^0) + N) S^0
+  //        = inv(I + S^0 N) S^0
+  //        = Q S^0
+  const double S11 = Q11*S0aa + Q12*S0ab + Q13*S0ac;
+  const double S12 = Q12*S0aa + Q22*S0ab + Q23*S0ac;
+  const double S13 = Q13*S0aa + Q23*S0ab + Q33*S0ac;
+  const double S22 = Q12*S0ab + Q22*S0bb + Q23*S0bc;
+  const double S23 = Q13*S0ab + Q23*S0bb + Q33*S0bc;
+  const double S33 = Q13*S0ac + Q23*S0bc + Q33*S0cc;
+  // If the full S is needed...it isn't since Ldd = (rho_d - rho_d) = 0 below
+  //const double S14=-S11-S12-S13;
+  //const double S24=-S12-S22-S23;
+  //const double S34=-S13-S23-S33;
+  //const double S44=S11+S22+S33 + 2.0*(S12+S13+S23);
+
+  // eq 12 of Akcasu, 1990: I(q) = L^T S L
+  // Note: eliminate cases without A and B polymers by setting Lij to 0
+  // Note: 1e-13 to convert from fm to cm for scattering length
+  const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13;
+  const double Lad = icase<5 ? 0.0 : (La/va - Ld/vd)*sqrt_Nav;
+  const double Lbd = icase<2 ? 0.0 : (Lb/vb - Ld/vd)*sqrt_Nav;
+  const double Lcd = (Lc/vc - Ld/vd)*sqrt_Nav;
+
+  const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33
+                    + 2.0*(Lad*Lbd*S12 + Lbd*Lcd*S23 + Lad*Lcd*S13);
+
+  return result;
+
+}
+
+double Iqxy(double qx, double qy,
+    double case_num,
+    double Na, double Phia, double va, double a_sld, double ba,
+    double Nb, double Phib, double vb, double b_sld, double bb,
+    double Nc, double Phic, double vc, double c_sld, double bc,
+    double Nd, double Phid, double vd, double d_sld, double bd,
+    double Kab, double Kac, double Kad,
+    double Kbc, double Kbd, double Kcd
+    )
+{
+    double q = sqrt(qx*qx + qy*qy);
+    return Iq(q,
+        case_num,
+        Na, Phia, va, a_sld, ba,
+        Nb, Phib, vb, b_sld, bb,
+        Nc, Phic, vc, c_sld, bc,
+        Nd, Phid, vd, d_sld, bd,
+        Kab, Kac, Kad,
+        Kbc, Kbd, Kcd);
+}
Index: sasmodels/models/rpa.py
===================================================================
--- sasmodels/models/rpa.py	(revision 82c299fa51ab091ee8ead0a235d4177ccfeecb99)
+++ sasmodels/models/rpa.py	(revision 82c299fa51ab091ee8ead0a235d4177ccfeecb99)
@@ -0,0 +1,144 @@
+r"""
+Calculates the macroscopic scattering intensity for a multi-component
+homogeneous mixture of polymers using the Random Phase Approximation.
+This general formalism contains 10 specific cases
+
+Case 0: C/D binary mixture of homopolymers
+
+Case 1: C-D diblock copolymer
+
+Case 2: B/C/D ternary mixture of homopolymers
+
+Case 3: C/C-D mixture of a homopolymer B and a diblock copolymer C-D
+
+Case 4: B-C-D triblock copolymer
+
+Case 5: A/B/C/D quaternary mixture of homopolymers
+
+Case 6: A/B/C-D mixture of two homopolymers A/B and a diblock C-D
+
+Case 7: A/B-C-D mixture of a homopolymer A and a triblock B-C-D
+
+Case 8: A-B/C-D mixture of two diblock copolymers A-B and C-D
+
+Case 9: A-B-C-D tetra-block copolymer
+
+**NB: these case numbers are different from those in the NIST SANS package!**
+
+Only one case can be used at any one time.
+
+The RPA (mean field) formalism only applies only when the multicomponent
+polymer mixture is in the homogeneous mixed-phase region.
+
+**Component D is assumed to be the "background" component (ie, all contrasts
+are calculated with respect to component D).** So the scattering contrast
+for a C/D blend = [SLD(component C) - SLD(component D)]\ :sup:`2`.
+
+Depending on which case is being used, the number of fitting parameters - the
+segment lengths (ba, bb, etc) and $\chi$ parameters (Kab, Kac, etc) - vary.
+The *scale* parameter should be held equal to unity.
+
+The input parameters are the degrees of polymerization, the volume fractions,
+the specific volumes, and the neutron scattering length densities for each
+component.
+
+.. figure:: img/image215.jpg
+
+    1D plot using the default values (w/500 data points).
+
+References
+----------
+
+A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136
+"""
+
+from numpy import inf
+
+name = "rpa"
+title = "Random Phase Approximation"
+description = """
+This formalism applies to multicomponent polymer mixtures in the
+homogeneous (mixed) phase region only.
+Case 0: C/D binary mixture of homopolymers
+Case 1: C-D diblock copolymer
+Case 2: B/C/D ternary mixture of homopolymers
+Case 3: B/C-D mixture of homopolymer b and diblock copolymer C-D
+Case 4: B-C-D triblock copolymer
+Case 5: A/B/C/D quaternary mixture of homopolymers
+Case 6: A/B/C-D mixture of two homopolymers A/B and a diblock C-D
+Case 7: A/B-C-D mixture of a homopolymer A and a triblock B-C-D
+Case 8: A-B/C-D mixture of two diblock copolymers A-B and C-D
+Case 9: A-B-C-D four-block copolymer
+See details in the model function help
+"""
+category = ""
+
+CASES = [
+    "C+D binary mixture",
+    "C:D diblock copolymer",
+    "B+C+D ternary mixture",
+    "B+C:D binary mixture",
+    "B:C:D triblock copolymer",
+    "A+B+C+D quaternary mixture",
+    "A+B+C:D ternary mixture",
+    "A+B:C:D binary mixture",
+    "A:B+C:D binary mixture",
+    "A:B:C:D quadblock copolymer",
+]
+
+#   ["name", "units", default, [lower, upper], "type","description"],
+parameters = [
+    ["case_num", CASES, 0, [0, 10], "", "Component organization" ],
+
+    ["Na", "", 1000.0, [1, inf], "", "Degree of polymerization"],
+    ["Phia", "", 0.25, [0, 1], "", "volume fraction"],
+    ["va", "mL/mol", 100.0, [0, inf], "", "specific volume"],
+    ["La", "fm", 10.0, [-inf, inf], "", "scattering length"],
+    ["ba", "Ang", 5.0, [0, inf], "", "segment length"],
+
+    ["Nb", "", 1000.0, [1, inf], "", "Degree of polymerization"],
+    ["Phib", "", 0.25, [0, 1], "", "volume fraction"],
+    ["vb", "mL/mol", 100.0, [0, inf], "", "specific volume"],
+    ["Lb", "fm", 10.0, [-inf, inf], "", "scattering length"],
+    ["bb", "Ang", 5.0, [0, inf], "", "segment length"],
+
+    ["Nc", "", 1000.0, [1, inf], "", "Degree of polymerization"],
+    ["Phic", "", 0.25, [0, 1], "", "volume fraction"],
+    ["vc", "mL/mol", 100.0, [0, inf], "", "specific volume"],
+    ["Lc", "fm", 10.0, [-inf, inf], "", "scattering length"],
+    ["bc", "Ang", 5.0, [0, inf], "", "segment length"],
+
+    ["Nd", "", 1000.0, [1, inf], "", "Degree of polymerization"],
+    ["Phid", "", 0.25, [0, 1], "", "volume fraction"],
+    ["vd", "mL/mol", 100.0, [0, inf], "", "specific volume"],
+    ["Ld", "fm", 10.0, [-inf, inf], "", "scattering length"],
+    ["bd", "Ang", 5.0, [0, inf], "", "segment length"],
+
+    ["Kab", "", -0.0004, [-inf, inf], "", "Interaction parameter"],
+    ["Kac", "", -0.0004, [-inf, inf], "", "Interaction parameter"],
+    ["Kad", "", -0.0004, [-inf, inf], "", "Interaction parameter"],
+    ["Kbc", "", -0.0004, [-inf, inf], "", "Interaction parameter"],
+    ["Kbd", "", -0.0004, [-inf, inf], "", "Interaction parameter"],
+    ["Kcd", "", -0.0004, [-inf, inf], "", "Interaction parameter"],
+]
+
+category = "shape-independent"
+
+source = ["rpa.c"]
+
+HIDE_NONE = set()
+HIDE_A = set("Na Phia va La Kab Kac Kad".split())
+HIDE_AB = set("Nb Phib vb Lb Kbc Kbd".split()).union(HIDE_A)
+def hidden(pars):
+    case_num = pars.get("case_num", parameters[0][2])
+    if case_num < 2:
+        return HIDE_A
+    elif case_num < 5:
+        return HIDE_AB
+    else:
+        return HIDE_NONE
+
+oldname = 'RPAModel'
+oldpars = dict(
+    case_num="lcase_n",
+)
Index: test/rpa.py
===================================================================
--- test/rpa.py	(revision 82c299fa51ab091ee8ead0a235d4177ccfeecb99)
+++ test/rpa.py	(revision 82c299fa51ab091ee8ead0a235d4177ccfeecb99)
@@ -0,0 +1,202 @@
+from collections import namedtuple
+import numpy as np
+from numpy import sqrt, exp, expm1
+
+AVOGADRO = 6.022e23
+
+Polymer = namedtuple("Polymer", "n phi v a b".split())
+def E(q, poly):
+    qrsq = (q*Rg(poly))**2
+    retval = exp(-qrsq)
+    return retval
+
+def F(q, poly):
+    qrsq = (q*Rg(poly))**2
+    retval = -expm1(-qrsq)/qrsq
+    return retval
+
+def P_ii(q, poly):
+    qrsq = (q*Rg(poly))**2
+    retval = 2 * (expm1(-qrsq) + qrsq)/qrsq**2
+    return retval
+
+def P_ij(q, poly_list):
+    i, j = poly_list[0], poly_list[-1]
+    retval = F(q, i) * np.prod([E(q,p) for p in poly_list[1:-1]]) * F(q, j)
+    return retval
+
+def Rg(poly):
+    return sqrt(poly.n/6.)*poly.a
+
+def S0_ii(q, poly):
+    retval = poly.n*poly.phi*poly.v*P_ii(q, poly)
+    return retval
+
+def S0_ij(q, poly_list):
+    i,j = poly_list[0], poly_list[-1]
+    retval = sqrt(i.n*i.phi*i.v*j.n*j.phi*j.v) * P_ij(q, poly_list)
+    return retval
+
+def drho(poly, base):
+    return (poly.b/poly.v -  base.b/base.v)*1e-13*sqrt(AVOGADRO)
+
+def binary_blend(q, C, D, Kcd):
+    """
+    de Gennes, 1979
+    """
+    S0cc = S0_ii(q, C)
+    S0dd = S0_ii(q, D)
+    vcc = 1/S0dd - 2*Kcd #/v0
+    Scc = S0cc/(1 + vcc*S0cc)
+    rhocd = drho(C,D)
+    Iq = rhocd**2 * Scc
+    return Iq
+
+def ternary_blend(q, B, C, D, Kbc, Kbd, Kcd):
+    S0bb = S0_ii(q, B)
+    S0cc = S0_ii(q, C)
+    S0dd = S0_ii(q, D)
+    vbb = 1/S0dd - 2*Kbd
+    vcc = 1/S0dd - 2*Kcd
+    vbc = 1/S0dd + Kbd - Kbc - Kcd
+    rhobd = drho(B,D)
+    rhocd = drho(C,D)
+    det = (1+vbb*S0bb)*(1+vcc*S0cc) - vbc**2*S0bb*S0cc
+    Sbb = S0bb*(1+vcc*S0cc)/det
+    Scc = S0cc*(1+vbb*S0bb)/det
+    Sbc = -S0bb*vbc*S0cc/det
+    Iq = rhobd**2*Sbb + rhocd**2*Scc + 2*rhobd*rhocd*Sbc
+    return Iq
+
+def diblock_copolymer(q, C, D, Kcd):
+    """
+    Leibler, 1980
+    """
+    S0cc = S0_ii(q, C)
+    S0dd = S0_ii(q, D)
+    S0cd = S0_ij(q, [C, D])
+    Scc = (S0cc*S0dd - S0cd**2)/((S0cc+S0dd + 2*S0cd)-2*Kcd*(S0cc+S0dd-2*S0cd))
+    rhocd = drho(C,D)
+    Iq = rhocd**2 * Scc
+    return Iq
+
+def matrix_form(q, case_num, polys,
+                Kab=0., Kac=0., Kad=0., Kbc=0., Kbd=0., Kcd=0.):
+    if case_num < 2:
+        C, D = polys
+        A = B = D
+    elif case_num < 5:
+        B, C, D = polys
+        A = D
+    else:
+        A, B, C, D = polys
+
+    rho = np.matrix([[drho(p, D)] for p in (A,B,C)])
+    S0aa = S0_ii(q, A)
+    S0bb = S0_ii(q, B)
+    S0cc = S0_ii(q, C)
+    S0ab = S0_ij(q, [A, B])
+    S0ac = S0_ij(q, [A, B, C])
+    S0bc = S0_ij(q, [B, C])
+    if case_num == 4: # No a-c interaction in triblock copolymer
+        S0ac *= 0.0
+    elif case_num == 9: # No a-c or a-d interaction in tetrablock copolymer
+        S0ac *= 0.0
+
+    S0 = np.array([[S0aa, S0ab, S0ac], [S0ab, S0bb, S0bc], [S0ac, S0bc, S0cc]])
+    S0dd = S0_ii(q, D)
+    vaa = 1./S0dd - 2*Kad
+    vbb = 1./S0dd - 2*Kbd
+    vcc = 1./S0dd - 2*Kcd
+    vab = 1./S0dd + Kab - Kad - Kbd
+    vac = 1./S0dd + Kac - Kad - Kcd
+    vbc = 1./S0dd + Kbc - Kbd - Kcd
+    v = np.array([[vaa, vab, vac], [vab, vbb, vbc], [vac, vbc, vcc]])
+
+    Iq = np.empty_like(q)
+    for k, qk in enumerate(q):
+        S0_k = S0[:,:,k].M
+        v_k = v[:,:,k].M
+        S_k = np.linalg.inv(1 + S0_k*v_k)*S0_k
+        Iq[k] = rho.T * S_k * rho
+
+def build_pars(case_num, polys, **interactions):
+    def set_poly(x, poly):
+        pars["N"+x] = poly.n
+        pars["Phi"+x] = poly.phi
+        pars["v"+x] = poly.v
+        pars["b"+x] = poly.a
+        pars["L"+x] = poly.b
+    pars = interactions.copy()
+    pars["case_num"] = case_num
+    polys = list(reversed(polys))
+    if len(polys) >= 4: set_poly("a",polys[3])
+    if len(polys) >= 3: set_poly("b",polys[2])
+    if len(polys) >= 2: set_poly("c",polys[1])
+    if len(polys) >= 1: set_poly("d",polys[0])
+    return pars
+
+def sasmodels_rpa(q, pars):
+    from sasmodels.models import rpa
+    from sasmodels.core import load_model
+    from sasmodels.direct_model import DirectModel
+    from sasmodels.data import empty_data1D
+    data = empty_data1D(q, resolution=0.0)
+    model = load_model(rpa, dtype="double", platform="dll")
+    #model = load_model(rpa, dtype="single", platform="ocl")
+    M = DirectModel(data, model)
+    return M(**pars)
+
+def sasview_rpa(q, pars):
+    from sasmodels.models import rpa
+    from sasmodels.compare import eval_sasview
+    from sasmodels.data import empty_data1D
+    data = empty_data1D(q, resolution=0.0)
+    M = eval_sasview(rpa, data)
+    return M(**pars)
+
+def demo():
+    import sys
+    case_num = 0 if len(sys.argv) < 2 else int(sys.argv[1])
+
+    B = Polymer(n=525,phi=0.57,v=97.5,b=-4.99,a=8)
+    C = Polymer(n=525,phi=0.57,v=97.5,b=-4.99,a=8)
+    D = Polymer(n=1105,phi=0.43,v=81.9,b=53.1,a=2)
+    q = np.logspace(-4,-1,400)
+    #q = np.array([0.1])
+    Kab=Kac=Kad=0.0
+    Kcd = 0.0106*0.0035 - 1.84e-5
+    Kbd = Kcd + 2e-4
+    Kbc = (Kcd + 1e-4)*0.5
+    K = dict(Kab=Kab,Kac=Kac,Kad=Kad,Kbc=Kbc,Kbd=Kbd,Kcd=Kcd)
+    if case_num == 0:
+        Iq = binary_blend(q, C, D, Kcd)
+    elif case_num == 1:
+        Iq = diblock_copolymer(q, C, D, Kcd)
+    elif case_num == 2:
+        Iq = ternary_blend(q, B, C, D, Kbc, Kbd, Kcd)
+    else:
+        raise ValueError("Case %d not implmented"%case_num)
+
+    pars = build_pars(case_num, [B, C, D], **K)
+    print "eval sasmodels"
+    Iq_sasmodels = sasmodels_rpa(q, pars)
+    print "eval sasview"
+    Iq_sasview = sasview_rpa(q, pars)
+    print 1./Iq[0], 1./Iq_sasmodels[0], 1./Iq_sasview[0]
+
+    #return
+    import pylab
+    pylab.subplot(121)
+    pylab.loglog(q, Iq, label='direct')
+    pylab.loglog(q, Iq_sasmodels, label='sasmodels')
+    pylab.loglog(q, Iq_sasview, label='sasview')
+    pylab.legend()
+    pylab.subplot(122)
+    pylab.loglog(q, abs(Iq_sasview - Iq_sasmodels)/Iq_sasmodels, label='sasview-sasmodels')
+    pylab.loglog(q, abs(Iq_sasmodels - Iq)/Iq, label='sasmodels-direct')
+    pylab.legend()
+    pylab.show()
+
+if __name__ == "__main__":
+    demo()
