Index: example/fit.py
===================================================================
--- example/fit.py	(revision 1182da5952c5f1fedd003e65e6b154df780d0afb)
+++ example/fit.py	(revision f4b36fa94da2129ef86fb36f21c6458873815e0f)
@@ -24,5 +24,6 @@
             % section)
 data = radial_data if section != "tangential" else tan_data
-phi = 0 if section != "tangential" else 90
+theta = 89.9 if section != "tangential" else 0
+phi = 90
 kernel = load_model(name, dtype="single")
 cutoff = 1e-3
@@ -30,34 +31,32 @@
 if name == "ellipsoid":
     model = Model(kernel,
-        scale=0.08,
-        r_polar=15, r_equatorial=800,
+        scale=0.08, background=35,
+        radius_polar=15, radius_equatorial=800,
         sld=.291, sld_solvent=7.105,
-        background=0,
-        theta=90, phi=phi,
-        theta_pd=15, theta_pd_n=40, theta_pd_nsigma=3,
-        r_polar_pd=0.222296, r_polar_pd_n=1, r_polar_pd_nsigma=0,
-        r_equatorial_pd=.000128, r_equatorial_pd_n=1, r_equatorial_pd_nsigma=0,
+        theta=theta, phi=phi,
+        theta_pd=0, theta_pd_n=0, theta_pd_nsigma=3,
         phi_pd=0, phi_pd_n=20, phi_pd_nsigma=3,
+        radius_polar_pd=0.222296, radius_polar_pd_n=1, radius_polar_pd_nsigma=0,
+        radius_equatorial_pd=.000128, radius_equatorial_pd_n=1, radius_equatorial_pd_nsigma=0,
         )
 
-
     # SET THE FITTING PARAMETERS
-    model.r_polar.range(15, 1000)
-    model.r_equatorial.range(15, 1000)
-    model.theta_pd.range(0, 360)
+    model.radius_polar.range(15, 1000)
+    model.radius_equatorial.range(15, 1000)
+    #model.theta.range(0, 90)
+    #model.theta_pd.range(0,10)
+    model.phi_pd.range(0,20)
+    model.phi.range(0, 180)
     model.background.range(0,1000)
     model.scale.range(0, 10)
 
 
-
 elif name == "lamellar":
     model = Model(kernel,
-        scale=0.08,
+        scale=0.08, background=0.003,
         thickness=19.2946,
         sld=5.38,sld_sol=7.105,
-        background=0.003,
         thickness_pd= 0.37765, thickness_pd_n=10, thickness_pd_nsigma=3,
         )
-
 
     # SET THE FITTING PARAMETERS
@@ -77,15 +76,15 @@
         radius_pd=.0084, radius_pd_n=10, radius_pd_nsigma=3,
         length_pd=0.493, length_pd_n=10, length_pd_nsigma=3,
-        phi_pd=0, phi_pd_n=5, phi_pd_nsigma=3,)
+        phi_pd=0, phi_pd_n=5 phi_pd_nsigma=3,)
         """
     pars = dict(
         scale=.01, background=35,
         sld=.291, sld_solvent=5.77,
-        radius=250, length=178, 
-        theta=90, phi=phi,
+        radius=250, length=178,
         radius_pd=0.1, radius_pd_n=5, radius_pd_nsigma=3,
         length_pd=0.1,length_pd_n=5, length_pd_nsigma=3,
-        theta_pd=10, theta_pd_n=50, theta_pd_nsigma=3,
-        phi_pd=0, phi_pd_n=10, phi_pd_nsigma=3)
+        theta=theta, phi=phi,
+        theta_pd=0, theta_pd_n=0, theta_pd_nsigma=3,
+        phi_pd=10, phi_pd_n=20, phi_pd_nsigma=3)
     model = Model(kernel, **pars)
 
@@ -93,9 +92,9 @@
     model.radius.range(1, 500)
     model.length.range(1, 5000)
-    model.theta.range(-90,100)
-    model.theta_pd.range(0, 30)
-    model.theta_pd_n = model.theta_pd + 5
+    #model.theta.range(0, 90)
+    model.phi.range(0, 180)
+    model.phi_pd.range(0, 30)
     model.radius_pd.range(0, 1)
-    model.length_pd.range(0, 2)
+    model.length_pd.range(0, 1)
     model.scale.range(0, 10)
     model.background.range(0, 100)
@@ -104,21 +103,22 @@
 elif name == "core_shell_cylinder":
     model = Model(kernel,
-        scale= .031, radius=19.5, thickness=30, length=22,
-        sld_core=7.105, sld_shell=.291, sdl_solvent=7.105,
-        background=0, theta=0, phi=phi,
-
+        scale= .031, background=0,
+        radius=19.5, thickness=30, length=22,
+        sld_core=7.105, sld_shell=.291, sld_solvent=7.105,
         radius_pd=0.26, radius_pd_n=10, radius_pd_nsigma=3,
         length_pd=0.26, length_pd_n=10, length_pd_nsigma=3,
         thickness_pd=1, thickness_pd_n=1, thickness_pd_nsigma=1,
-        theta_pd=1, theta_pd_n=10, theta_pd_nsigma=3,
-        phi_pd=0.1, phi_pd_n=1, phi_pd_nsigma=1,
+        theta=theta, phi=phi,
+        theta_pd=1, theta_pd_n=1, theta_pd_nsigma=3,
+        phi_pd=0, phi_pd_n=20, phi_pd_nsigma=3,
         )
 
     # SET THE FITTING PARAMETERS
-    #model.radius.range(115, 1000)
-    #model.length.range(0, 2500)
+    model.radius.range(115, 1000)
+    model.length.range(0, 2500)
     #model.thickness.range(18, 38)
     #model.thickness_pd.range(0, 1)
     #model.phi.range(0, 90)
+    model.phi_pd.range(0,20)
     #model.radius_pd.range(0, 1)
     #model.length_pd.range(0, 1)
@@ -131,14 +131,25 @@
 elif name == "capped_cylinder":
     model = Model(kernel,
-        scale=.08, radius=20, cap_radius=40, length=400,
+        scale=.08, background=35,
+        radius=20, cap_radius=40, length=400,
         sld=1, sld_solvent=6.3,
-        background=0, theta=0, phi=phi,
         radius_pd=.1, radius_pd_n=5, radius_pd_nsigma=3,
         cap_radius_pd=.1, cap_radius_pd_n=5, cap_radius_pd_nsigma=3,
         length_pd=.1, length_pd_n=1, length_pd_nsigma=0,
-        theta_pd=.1, theta_pd_n=1, theta_pd_nsigma=0,
-        phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0,
+        theta=theta, phi=phi,
+        theta_pd=0, theta_pd_n=1, theta_pd_nsigma=0,
+        phi_pd=10, phi_pd_n=20, phi_pd_nsigma=0,
         )
 
+    model.radius.range(115, 1000)
+    model.length.range(0, 2500)
+    #model.thickness.range(18, 38)
+    #model.thickness_pd.range(0, 1)
+    #model.phi.range(0, 90)
+    model.phi_pd.range(0,20)
+    #model.radius_pd.range(0, 1)
+    #model.length_pd.range(0, 1)
+    #model.theta_pd.range(0, 360)
+    #model.background.range(0,5)
     model.scale.range(0, 1)
 
@@ -146,19 +157,20 @@
 elif name == "triaxial_ellipsoid":
     model = Model(kernel,
-        scale=0.08, req_minor=15, req_major=20, rpolar=500,
+        scale=0.08, background=35,
+        radius_equat_minor=15, radius_equat_major=20, radius_polar=500,
         sld=7.105, solvent_sld=.291,
-        background=5, theta=0, phi=phi, psi=0,
+        radius_equat_minor_pd=.1, radius_equat_minor_pd_n=1, radius_equat_minor_pd_nsigma=0,
+        radius_equat_major_pd=.1, radius_equat_major_pd_n=1, radius_equat_major_pd_nsigma=0,
+        radius_polar_pd=.1, radius_polar_pd_n=1, radius_polar_pd_nsigma=0,
+        theta=theta, phi=phi, psi=0,
         theta_pd=20, theta_pd_n=40, theta_pd_nsigma=3,
         phi_pd=.1, phi_pd_n=1, phi_pd_nsigma=0,
         psi_pd=30, psi_pd_n=1, psi_pd_nsigma=0,
-        req_minor_pd=.1, req_minor_pd_n=1, req_minor_pd_nsigma=0,
-        req_major_pd=.1, req_major_pd_n=1, req_major_pd_nsigma=0,
-        rpolar_pd=.1, rpolar_pd_n=1, rpolar_pd_nsigma=0,
         )
 
     # SET THE FITTING PARAMETERS
-    model.req_minor.range(15, 1000)
-    model.req_major.range(15, 1000)
-    #model.rpolar.range(15, 1000)
+    model.radius_equat_minor.range(15, 1000)
+    model.radius_equat_major.range(15, 1000)
+    #model.radius_polar.range(15, 1000)
     #model.background.range(0,1000)
     #model.theta_pd.range(0, 360)
@@ -173,5 +185,5 @@
 M = Experiment(data=data, model=model)
 if section == "both":
-   tan_model = Model(model.kernel, **model.parameters())
+   tan_model = Model(model.sasmodel, **model.parameters())
    tan_model.phi = model.phi - 90
    tan_model.cutoff = cutoff
Index: sasmodels/models/rpa.c
===================================================================
--- sasmodels/models/rpa.c	(revision fa4a99420607d94c8c029525a1cb0814491c18aa)
+++ sasmodels/models/rpa.c	(revision 20c856a138b7cfbe5d8eb03e1655371ef25b4b06)
@@ -39,5 +39,6 @@
   double S31,S32,S33,S34,S41,S42,S43,S44;
   double Lad,Lbd,Lcd,Nav,Intg;
-
+  
+  // Set values for non existent parameters (eg. no A or B in case 0 and 1 etc)
   //icase was shifted to N-1 from the original code
   if (icase <= 1){
@@ -57,5 +58,9 @@
     Kab = Kac = Kad = -0.0004;
   }
-
+ 
+  // Set volume fraction of component D based on constraint that sum of vol frac =1
+  Phi[3]=1.0-Phi[0]-Phi[1]-Phi[2];
+
+  //set up values for cross terms in case of block copolymers (1,3,4,6,7,8,9)
   Nab=sqrt(N[0]*N[1]);
   Nac=sqrt(N[0]*N[2]);
@@ -79,4 +84,5 @@
   Phicd=sqrt(Phi[2]*Phi[3]);
 
+  // Calculate Q^2 * Rg^2 for each homopolymer assuming random walk 
   Xa=q*q*b[0]*b[0]*N[0]/6.0;
   Xb=q*q*b[1]*b[1]*N[1]/6.0;
@@ -84,26 +90,27 @@
   Xd=q*q*b[3]*b[3]*N[3]/6.0;
 
-  Paa=2.0*(exp(-Xa)-1.0+Xa)/(Xa*Xa);
-  S0aa=N[0]*Phi[0]*v[0]*Paa;
-  Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb);
+  //calculate all partial structure factors Pij and normalize n^2
+  Paa=2.0*(exp(-Xa)-1.0+Xa)/(Xa*Xa); // free A chain form factor
+  S0aa=N[0]*Phi[0]*v[0]*Paa; // Phi * Vp * P(Q)= I(Q0)/delRho^2
+  Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb); //AB diblock (anchored Paa * anchored Pbb) partial form factor
   S0ab=(Phiab*vab*Nab)*Pab;
-  Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc);
+  Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); //ABC triblock AC partial form factor
   S0ac=(Phiac*vac*Nac)*Pac;
-  Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd);
+  Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); //ABCD four block
   S0ad=(Phiad*vad*Nad)*Pad;
 
   S0ba=S0ab;
-  Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb);
+  Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); // free B chain
   S0bb=N[1]*Phi[1]*v[1]*Pbb;
-  Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc);
+  Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); // BC diblock
   S0bc=(Phibc*vbc*Nbc)*Pbc;
-  Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd);
+  Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); // BCD triblock
   S0bd=(Phibd*vbd*Nbd)*Pbd;
 
   S0ca=S0ac;
   S0cb=S0bc;
-  Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc);
+  Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); // Free C chain
   S0cc=N[2]*Phi[2]*v[2]*Pcc;
-  Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd);
+  Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); // CD diblock
   S0cd=(Phicd*vcd*Ncd)*Pcd;
 
@@ -111,6 +118,8 @@
   S0db=S0bd;
   S0dc=S0cd;
-  Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd);
+  Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); // free D chain
   S0dd=N[3]*Phi[3]*v[3]*Pdd;
+
+  // Reset all unused partial structure factors to 0 (depends on case)
   //icase was shifted to N-1 from the original code
   switch(icase){
@@ -193,4 +202,5 @@
   S0dc=S0cd;
 
+  // self chi parameter is 0 ... of course
   Kaa=0.0;
   Kbb=0.0;
@@ -243,4 +253,5 @@
   ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd);
 
+  // D is considered the matrix or background component so enters here
   m=1.0/(S0dd-ZZ);
 
@@ -297,4 +308,8 @@
   S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23;
 
+  //calculate contrast where L[i] is the scattering length of i and D is the matrix
+  //Note that should multiply by Nav to get units of SLD which will become
+  // Nav*2 in the next line (SLD^2) but then normalization in that line would
+  //need to divide by Nav leaving only Nav or sqrt(Nav) before squaring. 
   Nav=6.022045e+23;
   Lad=(L[0]/v[0]-L[3]/v[3])*sqrt(Nav);
@@ -303,4 +318,7 @@
 
   Intg=Lad*Lad*S11+Lbd*Lbd*S22+Lcd*Lcd*S33+2.0*Lad*Lbd*S12+2.0*Lbd*Lcd*S23+2.0*Lad*Lcd*S13;
+
+  //rescale for units of Lij^2 (fm^2 to cm^2)
+  Intg *= 1.0e-26;    
 
   return Intg;
Index: sasmodels/models/rpa.py
===================================================================
--- sasmodels/models/rpa.py	(revision fa4a99420607d94c8c029525a1cb0814491c18aa)
+++ sasmodels/models/rpa.py	(revision 20c856a138b7cfbe5d8eb03e1655371ef25b4b06)
@@ -1,3 +1,6 @@
 r"""
+Definition
+----------
+
 Calculates the macroscopic scattering intensity for a multi-component
 homogeneous mixture of polymers using the Random Phase Approximation.
@@ -24,22 +27,21 @@
 Case 9: A-B-C-D tetra-block copolymer
 
-**NB: these case numbers are different from those in the NIST SANS package!**
+.. note::
+    These case numbers are different from those in the NIST SANS package!
 
-Only one case can be used at any one time.
+USAGE NOTES:
 
-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.
+* 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 can 
+  vary.  Note that in general the degrees of polymerization, the volume
+  fractions, the molar volumes, and the neutron scattering lengths for each
+  component are obtained from other methods and held fixed while the segment
+  lengths (b\ :sub:`a`, b\ :sub:`b`, etc) and $\chi$ parameters (K\ :sub:`ab`,
+  K\ :sub:`ac`, etc). The *scale* parameter should be held equal to unity.
 
 
@@ -47,5 +49,5 @@
 ----------
 
-A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136
+.. [#] A Z Akcasu, R Klein and B Hammouda, *Macromolecules*, 26 (1993) 4136
 """
 
@@ -53,5 +55,5 @@
 
 name = "rpa"
-title = "Random Phase Approximation - unfinished work in progress"
+title = "Random Phase Approximation"
 description = """
 This formalism applies to multicomponent polymer mixtures in the
@@ -90,5 +92,5 @@
     ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"],
     ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"],
-    ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"],
+    ["v[4]", "mL/mol", 100.0, [0, inf], "", "molar volume"],
     ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"],
     ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"],
@@ -107,6 +109,6 @@
 
 control = "case_num"
-HIDE_NONE = set()
-HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split())
+HIDE_ALL = set("Phi4".split())
+HIDE_A = set("N1 Phi1 v1 L1 b1 K12 K13 K14".split()).union(HIDE_ALL)
 HIDE_AB = set("N2 Phi2 v2 L2 b2 K23 K24".split()).union(HIDE_A)
 def hidden(case_num):
@@ -120,4 +122,4 @@
         return HIDE_A
     else:
-        return HIDE_NONE
+        return HIDE_ALL
 
