Index: src/sans/models/c_extension/c_models/RectangularHollowPrism.cpp
===================================================================
--- src/sans/models/c_extension/c_models/RectangularHollowPrism.cpp	(revision 201af9f481a69dc31baac2847a064fb4036d1223)
+++ src/sans/models/c_extension/c_models/RectangularHollowPrism.cpp	(revision 201af9f481a69dc31baac2847a064fb4036d1223)
@@ -0,0 +1,199 @@
+/**
+    This software was developed by the University of Tennessee as part of the
+    Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+    project funded by the US National Science Foundation.
+
+    If you use DANSE applications to do scientific research that leads to
+    publication, we ask that you acknowledge the use of the software with the
+    following sentence:
+
+    "This work benefited from DANSE software developed under NSF award DMR-0520547."
+
+    copyright 2008, University of Tennessee
+ */
+
+/**
+ * Scattering model classes
+ * The classes use the IGOR library found in
+ *   sansmodels/src/libigor
+ *
+ *    TODO: add 2D function
+ */
+
+#include <math.h>
+#include "parameters.hh"
+#include <stdio.h>
+#include <iostream>
+using namespace std;
+
+extern "C" {
+    #include "libCylinder.h"
+    #include "libStructureFactor.h"
+    #include "libmultifunc/libfunc.h"
+}
+#include "RectangularHollowPrism.h"
+
+// Convenience parameter structure
+typedef struct {
+    double scale;
+    double short_side;
+    double b2a_ratio;
+    double c2a_ratio;
+    double thickness;
+    double sldPipe;
+    double sldSolv;
+    double background;
+} RectangularHollowPrismParameters;
+
+
+RectangularHollowPrismModel :: RectangularHollowPrismModel() {
+    scale      = Parameter(1.0);
+    short_side = Parameter(35.0, true);
+    short_side.set_min(1.0);
+    b2a_ratio   = Parameter(1.0, true);
+    b2a_ratio.set_min(1.0);
+    c2a_ratio   = Parameter(1.0, true);
+    c2a_ratio.set_min(1.0);
+    thickness = Parameter(1.0, true);
+    thickness.set_min(0.0);
+    sldPipe   = Parameter(6.3e-6);
+    sldSolv   = Parameter(1.0e-6);
+    background = Parameter(0.0);
+}
+
+/**
+ * Function to evaluate 1D scattering function
+ * @param q: q-value
+ * @return: function value
+ */
+double RectangularHollowPrismModel :: operator()(double q) {
+
+    double dp[8];
+
+    // Fill parameter array for IGOR library
+    // Add the background after averaging
+    dp[0] = scale();
+    dp[1] = short_side();
+    dp[2] = b2a_ratio();
+    dp[3] = c2a_ratio();
+    dp[4] = thickness();
+    dp[5] = sldPipe();
+    dp[6] = sldSolv();
+    dp[7] = 0.0;
+
+    // Get the dispersion points for a
+    vector<WeightPoint> weights_short_side;
+    short_side.get_weights(weights_short_side);
+
+    // Get the dispersion points for b/a ratio
+    vector<WeightPoint> weights_b2a_ratio;
+    b2a_ratio.get_weights(weights_b2a_ratio);
+
+    // Get the dispersion points for c/a ratio
+    vector<WeightPoint> weights_c2a_ratio;
+    c2a_ratio.get_weights(weights_c2a_ratio);
+
+    // Get the dispersion points for the thickness
+    vector<WeightPoint> weights_thickness;
+    thickness.get_weights(weights_thickness);
+
+    // Perform the computation, with all weight points
+    double sum = 0.0;
+    double norm = 0.0;
+    double vol = 0.0;
+
+    // Loop over short_side weight points
+    for (int i=0; i < (int)weights_short_side.size(); i++) {
+
+        dp[1] = weights_short_side[i].value;
+
+        // Loop over b/a ratios
+        for (int j=0; j < (int)weights_b2a_ratio.size(); j++) {
+
+            dp[2] = weights_b2a_ratio[j].value;
+
+            // Loop over c/a ratios
+            for (int k=0; k < (int)weights_c2a_ratio.size(); k++) {
+
+                dp[3] = weights_c2a_ratio[k].value;
+
+                // Loop over thicknesses
+                for (int l=0; l < (int)weights_thickness.size(); l++) {
+
+                    dp[4] = weights_thickness[l].value;
+
+        	        // Un-normalize  by volume = abc - (a-2*t)*(b-2*t)*(c-2*t)
+
+		            double a = dp[1];
+		            double b = dp[1] * dp[2];
+		            double c = dp[1] * dp[3];
+		            double t = dp[4];
+		            double vol_i = a*b*c - (a-2*t)*(b-2*t)*(c-2*t);
+
+        	        sum += weights_short_side[i].weight *
+        	               weights_b2a_ratio[j].weight *
+        	               weights_c2a_ratio[k].weight *
+			               weights_thickness[l].weight *
+        	               RectangularHollowPrism(dp, q) *
+        	               vol_i;
+
+        	        //Find average volume (ABC)
+
+        	        vol += weights_short_side[i].weight *
+        	               weights_b2a_ratio[j].weight *
+        	               weights_c2a_ratio[k].weight *
+			               weights_thickness[l].weight *
+        	               vol_i;
+
+        	        norm += weights_short_side[i].weight *
+                	        weights_b2a_ratio[j].weight *
+                	        weights_c2a_ratio[k].weight *
+			                weights_thickness[l].weight;
+                }
+            }
+	    }
+    }
+
+    if (vol != 0.0 && norm != 0.0) {
+        //Re-normalize by avg volume
+        sum = sum/(vol/norm);}
+
+    return sum/norm + background();
+    
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param q_x: value of Q along x
+ * @param q_y: value of Q along y
+ * @return: function value
+ */
+double RectangularHollowPrismModel :: operator()(double qx, double qy) {
+    return 1.0;
+}
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double RectangularHollowPrismModel :: evaluate_rphi(double q, double phi) {
+    double qx = q*cos(phi);
+    double qy = q*sin(phi);
+    return (*this).operator()(qx, qy);
+}
+
+/**
+ * Function to calculate effective radius
+ * @return: effective radius value
+ */
+double RectangularHollowPrismModel :: calculate_ER() {
+    return 1.0;
+
+}
+double RectangularHollowPrismModel :: calculate_VR() {
+  return 1.0;
+}
Index: src/sans/models/c_extension/c_models/RectangularHollowPrismInfThinWalls.cpp
===================================================================
--- src/sans/models/c_extension/c_models/RectangularHollowPrismInfThinWalls.cpp	(revision 201af9f481a69dc31baac2847a064fb4036d1223)
+++ src/sans/models/c_extension/c_models/RectangularHollowPrismInfThinWalls.cpp	(revision 201af9f481a69dc31baac2847a064fb4036d1223)
@@ -0,0 +1,181 @@
+/**
+    This software was developed by the University of Tennessee as part of the
+    Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+    project funded by the US National Science Foundation.
+
+    If you use DANSE applications to do scientific research that leads to
+    publication, we ask that you acknowledge the use of the software with the
+    following sentence:
+
+    "This work benefited from DANSE software developed under NSF award DMR-0520547."
+
+    copyright 2008, University of Tennessee
+ */
+
+/**
+ * Scattering model classes
+ * The classes use the IGOR library found in
+ *   sansmodels/src/libigor
+ *
+ *    TODO: add 2D function
+ */
+
+#include <math.h>
+#include "parameters.hh"
+#include <stdio.h>
+#include <iostream>
+using namespace std;
+
+extern "C" {
+    #include "libCylinder.h"
+    #include "libStructureFactor.h"
+    #include "libmultifunc/libfunc.h"
+}
+#include "RectangularHollowPrismInfThinWalls.h"
+
+// Convenience parameter structure
+typedef struct {
+    double scale;
+    double short_side;
+    double b2a_ratio;
+    double c2a_ratio;
+    double sldPipe;
+    double sldSolv;
+    double background;
+} RectangularHollowPrismInfThinWallsParameters;
+
+
+RectangularHollowPrismInfThinWallsModel :: RectangularHollowPrismInfThinWallsModel() {
+    scale      = Parameter(1.0);
+    short_side = Parameter(35.0, true);
+    short_side.set_min(1.0);
+    b2a_ratio   = Parameter(1.0, true);
+    b2a_ratio.set_min(1.0);
+    c2a_ratio   = Parameter(1.0, true);
+    c2a_ratio.set_min(1.0);
+    sldPipe   = Parameter(6.3e-6);
+    sldSolv   = Parameter(1.0e-6);
+    background = Parameter(0.0);
+}
+
+/**
+ * Function to evaluate 1D scattering function
+ * @param q: q-value
+ * @return: function value
+ */
+double RectangularHollowPrismInfThinWallsModel :: operator()(double q) {
+
+    double dp[7];
+
+    // Fill parameter array for IGOR library
+    // Add the background after averaging
+    dp[0] = scale();
+    dp[1] = short_side();
+    dp[2] = b2a_ratio();
+    dp[3] = c2a_ratio();
+    dp[4] = sldPipe();
+    dp[5] = sldSolv();
+    dp[6] = 0.0;
+
+    // Get the dispersion points for a
+    vector<WeightPoint> weights_short_side;
+    short_side.get_weights(weights_short_side);
+
+    // Get the dispersion points for b/a ratio
+    vector<WeightPoint> weights_b2a_ratio;
+    b2a_ratio.get_weights(weights_b2a_ratio);
+
+    // Get the dispersion points for c/a ratio
+    vector<WeightPoint> weights_c2a_ratio;
+    c2a_ratio.get_weights(weights_c2a_ratio);
+
+    // Perform the computation, with all weight points
+    double sum = 0.0;
+    double norm = 0.0;
+    double vol = 0.0;
+
+    // Loop over short_side weight points
+    for (int i=0; i < (int)weights_short_side.size(); i++) {
+
+        dp[1] = weights_short_side[i].value;
+
+        // Loop over b/a ratios
+        for (int j=0; j < (int)weights_b2a_ratio.size(); j++) {
+
+            dp[2] = weights_b2a_ratio[j].value;
+
+            // Loop over c/a ratios
+            for (int k=0; k < (int)weights_c2a_ratio.size(); k++) {
+
+                dp[3] = weights_c2a_ratio[k].value;
+
+        	    // Un-normalize  by volume = 2*a*b + 2*a*c + 2*b*c
+
+	            double a = dp[1];
+	            double b = dp[1] * dp[2];
+	            double c = dp[1] * dp[3];
+                double vol_i = (2.0*a*b) + (2.0*a*c) + (2.0*b*c);
+
+        	    sum += weights_short_side[i].weight *
+        	           weights_b2a_ratio[j].weight *
+        	           weights_c2a_ratio[k].weight *
+        	           RectangularHollowPrismInfThinWalls(dp, q) *
+                       vol_i;
+
+        	    //Find average volume (ABC)
+
+        	    vol += weights_short_side[i].weight *
+        	           weights_b2a_ratio[j].weight *
+        	           weights_c2a_ratio[k].weight *
+                       vol_i;
+
+        	    norm += weights_short_side[i].weight *
+                	    weights_b2a_ratio[j].weight *
+                	    weights_c2a_ratio[k].weight;
+            }
+        }                    
+    }
+
+    if (vol != 0.0 && norm != 0.0) {
+        //Re-normalize by avg volume
+        sum = sum/(vol/norm);}
+
+    return sum/norm + background();
+    
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param q_x: value of Q along x
+ * @param q_y: value of Q along y
+ * @return: function value
+ */
+double RectangularHollowPrismInfThinWallsModel :: operator()(double qx, double qy) {
+    return 1.0;
+}
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double RectangularHollowPrismInfThinWallsModel :: evaluate_rphi(double q, double phi) {
+    double qx = q*cos(phi);
+    double qy = q*sin(phi);
+    return (*this).operator()(qx, qy);
+}
+
+/**
+ * Function to calculate effective radius
+ * @return: effective radius value
+ */
+double RectangularHollowPrismInfThinWallsModel :: calculate_ER() {
+    return 1.0;
+
+}
+double RectangularHollowPrismInfThinWallsModel :: calculate_VR() {
+  return 1.0;
+}
Index: src/sans/models/c_extension/c_models/RectangularPrism.cpp
===================================================================
--- src/sans/models/c_extension/c_models/RectangularPrism.cpp	(revision 201af9f481a69dc31baac2847a064fb4036d1223)
+++ src/sans/models/c_extension/c_models/RectangularPrism.cpp	(revision 201af9f481a69dc31baac2847a064fb4036d1223)
@@ -0,0 +1,177 @@
+/**
+    This software was developed by the University of Tennessee as part of the
+    Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
+    project funded by the US National Science Foundation.
+
+    If you use DANSE applications to do scientific research that leads to
+    publication, we ask that you acknowledge the use of the software with the
+    following sentence:
+
+    "This work benefited from DANSE software developed under NSF award DMR-0520547."
+
+    copyright 2008, University of Tennessee
+ */
+
+/**
+ * Scattering model classes
+ * The classes use the IGOR library found in
+ *   sansmodels/src/libigor
+ *
+ *    TODO: add 2D function
+ */
+
+#include <math.h>
+#include "parameters.hh"
+#include <stdio.h>
+#include <iostream>
+using namespace std;
+
+extern "C" {
+    #include "libCylinder.h"
+    #include "libStructureFactor.h"
+    #include "libmultifunc/libfunc.h"
+}
+#include "RectangularPrism.h"
+
+// Convenience parameter structure
+typedef struct {
+    double scale;
+    double short_side;
+    double b2a_ratio;
+    double c2a_ratio;
+    double sldPipe;
+    double sldSolv;
+    double background;
+} RectangularPrismParameters;
+
+
+RectangularPrismModel :: RectangularPrismModel() {
+    scale      = Parameter(1.0);
+    short_side = Parameter(35.0, true);
+    short_side.set_min(1.0);
+    b2a_ratio   = Parameter(1.0, true);
+    b2a_ratio.set_min(1.0);
+    c2a_ratio   = Parameter(1.0, true);
+    c2a_ratio.set_min(1.0);
+    sldPipe   = Parameter(6.3e-6);
+    sldSolv   = Parameter(1.0e-6);
+    background = Parameter(0.0);
+}
+
+/**
+ * Function to evaluate 1D scattering function
+ * @param q: q-value
+ * @return: function value
+ */
+double RectangularPrismModel :: operator()(double q) {
+
+    double dp[7];
+
+    // Fill parameter array for IGOR library
+    // Add the background after averaging
+    dp[0] = scale();
+    dp[1] = short_side();
+    dp[2] = b2a_ratio();
+    dp[3] = c2a_ratio();
+    dp[4] = sldPipe();
+    dp[5] = sldSolv();
+    dp[6] = 0.0;
+
+    // Get the dispersion points for a
+    vector<WeightPoint> weights_short_side;
+    short_side.get_weights(weights_short_side);
+
+    // Get the dispersion points for b/a ratio
+    vector<WeightPoint> weights_b2a_ratio;
+    b2a_ratio.get_weights(weights_b2a_ratio);
+
+    // Get the dispersion points for c/a ratio
+    vector<WeightPoint> weights_c2a_ratio;
+    c2a_ratio.get_weights(weights_c2a_ratio);
+
+    // Perform the computation, with all weight points
+    double sum = 0.0;
+    double norm = 0.0;
+    double vol = 0.0;
+
+    // Loop over short_side weight points
+    for (int i=0; i < (int)weights_short_side.size(); i++) {
+
+        dp[1] = weights_short_side[i].value;
+
+        // Loop over b/a ratios
+        for (int j=0; j < (int)weights_b2a_ratio.size(); j++) {
+
+            dp[2] = weights_b2a_ratio[j].value;
+
+            // Loop over c/a ratios
+            for (int k=0; k < (int)weights_c2a_ratio.size(); k++) {
+
+                dp[3] = weights_c2a_ratio[k].value;
+
+        	    // Un-normalize  by volume = a * (a * b/a) * (a * c/a)
+                double vol_i = dp[1] * dp[1] * dp[2] * dp[1] * dp[3];
+
+        	    sum += weights_short_side[i].weight *
+        	           weights_b2a_ratio[j].weight *
+        	           weights_c2a_ratio[k].weight *
+        	           RectangularPrism(dp, q) *
+        	           vol_i;
+
+        	    //Find average volume (ABC)
+
+        	    vol += weights_short_side[i].weight *
+        	           weights_b2a_ratio[j].weight *
+        	           weights_c2a_ratio[k].weight *
+        	           vol_i;
+
+        	    norm += weights_short_side[i].weight *
+                	    weights_b2a_ratio[j].weight *
+                	    weights_c2a_ratio[k].weight;
+            }
+        }                    
+    }
+
+    if (vol != 0.0 && norm != 0.0) {
+        //Re-normalize by avg volume
+        sum = sum/(vol/norm);}
+
+    return sum/norm + background();
+    
+}
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param q_x: value of Q along x
+ * @param q_y: value of Q along y
+ * @return: function value
+ */
+double RectangularPrismModel :: operator()(double qx, double qy) {
+    return 1.0;
+}
+
+
+/**
+ * Function to evaluate 2D scattering function
+ * @param pars: parameters of the cylinder
+ * @param q: q-value
+ * @param phi: angle phi
+ * @return: function value
+ */
+double RectangularPrismModel :: evaluate_rphi(double q, double phi) {
+    double qx = q*cos(phi);
+    double qy = q*sin(phi);
+    return (*this).operator()(qx, qy);
+}
+
+/**
+ * Function to calculate effective radius
+ * @return: effective radius value
+ */
+double RectangularPrismModel :: calculate_ER() {
+    return 1.0;
+
+}
+double RectangularPrismModel :: calculate_VR() {
+  return 1.0;
+}
Index: src/sans/models/c_extension/libigor/libCylinder.c
===================================================================
--- src/sans/models/c_extension/libigor/libCylinder.c	(revision 230f479e04002ef42bafef9a227d586459b717c1)
+++ src/sans/models/c_extension/libigor/libCylinder.c	(revision 201af9f481a69dc31baac2847a064fb4036d1223)
@@ -3310,2 +3310,307 @@
 }
 
+
+/*
+ * Massive rectangular prism (a <= b <=c) 
+ * using eqns. (12)+(16) from Nayuk & Huber, Z. Phys. Chem. 226, 837 (2012).
+ * It is totally equivalent to the Parallelepiped model in this library,
+ * but instead of using as input the a, b, and c sides of the prism,
+ * it uses a and the ratios (b/a) and (c/a).
+ * This allows to keep the shape when using the polydispersity.
+ */
+
+double
+RectangularPrism(double dp[], double q)
+{
+    int i, j;
+    double scale, aa, bb, cc, delrho, bkg;		
+    int nordi=76;			        //order of integration
+    int nordj=76;
+    double Pi;
+    double sumi, sumj;
+    double v1a, v1b, v2a, v2b;		// upper and lower integration limits
+    double answer;
+    double theta, phi, arg, termA, termB, termC, AP, ahalf, bhalf, chalf;
+    double vol, sldp, sld;
+
+    Pi = 4.0*atan(1.0);
+
+    scale = dp[0];	
+    aa = dp[1];          // parameter short_side
+    bb = aa * dp[2];     // parameter b/a ratio
+    cc = aa * dp[3];     // parameter c/a ratio
+    sldp = dp[4];        // scattering length density of the object
+    sld = dp[5];         // scattering length density of the solvent
+    delrho = sldp - sld;
+    bkg = dp[6];
+
+    vol = aa*bb*cc;
+    ahalf = 0.5 * aa;
+    bhalf = 0.5 * bb;
+    chalf = 0.5 * cc;
+
+   //Integration limits to use in Gaussian quadrature
+    v1a = 0.;
+    v1b = Pi/2.;  //theta integration limits
+    v2a = 0.;
+    v2b = Pi/2.;  //phi integration limits
+
+    sumi = 0.0;
+    for(i=0;i<nordi;i++) {
+
+	    theta = ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b )/2.0;	//z values for gaussian quadrature
+
+	    arg = q*chalf*cos(theta);
+	    if (fabs(arg) > 1.e-16) {termC = sin(arg) / arg;} else {termC = 1.0;}  
+
+	    sumj = 0.0;
+	    for(j=0;j<nordj;j++) {
+
+                phi = ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b )/2.0; //z values for gaussian quadrature
+
+	        // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0
+
+	        arg = q*ahalf*sin(theta)*sin(phi); 
+	        if (fabs(arg) > 1.e-16) {termA = sin(arg) / arg;} else {termA = 1.0;}
+	       
+	        arg = q*bhalf*sin(theta)*cos(phi); 
+	        if (fabs(arg) > 1.e-16) {termB = sin(arg) / arg;} else {termB = 1.0;}	  
+               
+	        AP = termA * termB * termC;  
+
+	        sumj += Gauss76Wt[j] * (AP*AP);
+
+	    }
+
+	    sumj = (v2b-v2a) * sumj / 2.0;
+	    sumi += Gauss76Wt[i] * sumj * sin(theta);
+
+    }
+
+    answer = (v1b-v1a)/2.0*sumi;
+
+    // Normalize by Pi (Eqn. 16). 
+    // The term (ABC)^2 does not appear because it was introduced before from the defs of termA, termB, termC
+    // The factor 2 appears because the theta integral has been defined between 0 and pi/2, instead of 0 to pi
+    answer = 2.0 * answer / Pi; //Form factor P(q)
+
+    // Multiply by contrast^2
+    answer *= delrho*delrho;
+
+    //normalize by volume
+    answer *= vol;
+	
+    //convert to [cm-1]
+    answer *= 1.0e8;
+
+    //Scale
+    answer *= scale;
+
+    // add in the background
+    answer += bkg;
+
+    return answer;
+}
+
+
+/*
+ * Hollow rectangular prism (a <= b <=c) with infinitely thin walls
+ * using eqn. (7-11) from Nayuk & Huber, Z. Phys. Chem. 226, 837 (2012)
+ */
+
+double
+RectangularHollowPrismInfThinWalls(double dp[], double q)
+{
+    int i,j;
+    double scale,aa,bb,cc,delrho,bkg;
+    int nordi=76;
+    int nordj=76;
+    double Pi;
+    double sumi, sumj;
+    double v1a, v1b, v2a, v2b;		//upper and lower integration limits
+    double answer;
+    double theta, phi, termAL_theta, termAT_theta, AL, AT, ahalf, bhalf, chalf, vol;
+    double sldp, sld ;
+
+    Pi = 4.0*atan(1.0);
+
+    scale = dp[0];		  //make local copies in case memory moves
+    aa = dp[1];
+    bb = aa * dp[2];
+    cc = aa * dp[3];
+    sldp = dp[4];        // scattering length density of the object
+    sld = dp[5];         // scattering length density of the solvent
+    delrho = sldp - sld;
+    bkg = dp[6];
+
+    ahalf = 0.5 * aa;
+    bhalf = 0.5 * bb;
+    chalf = 0.5 * cc;
+
+    //Integration limits to use in Gaussian quadrature
+    v1a = 0.;
+    v1b = Pi/2;   //theta integration limits
+    v2a = 0.;
+    v2b = Pi/2.;  //phi integration limits
+
+    sumi = 0.0;
+    for(i=0;i<nordi;i++) {
+
+        theta = ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b )/2.0;	//z values for gaussian quadrature
+
+        // To check potential problems if denominator goes to zero here !!!
+        termAL_theta = 8.0*cos(q*chalf*cos(theta)) / (q*q*sin(theta)*sin(theta));
+        termAT_theta = 8.0*sin(q*chalf*cos(theta)) / (q*q*sin(theta)*cos(theta));
+
+        sumj = 0.0;
+        for(j=0;j<nordj;j++) {
+
+            phi = ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b )/2.0;	//z values for gaussian quadrature
+
+            // Amplitude AL from eqn. (7c)
+            AL = termAL_theta * sin(q*ahalf*sin(theta)*sin(phi)) * sin(q*bhalf*sin(theta)*cos(phi)) / (sin(phi)*cos(phi));
+
+            // Amplitude AT from eqn. (9)
+            AT = termAT_theta * (  (cos(q*ahalf*sin(theta)*sin(phi))*sin(q*bhalf*sin(theta)*cos(phi))/cos(phi)) + (cos(q*bhalf*sin(theta)*cos(phi))*sin(q*ahalf*sin(theta)*sin(phi))/sin(phi)) );
+
+            sumj += Gauss76Wt[j] * (AL+AT)*(AL+AT);
+
+        }
+
+        sumj = (v2b-v2a) * sumj / 2.0;
+        sumi += Gauss76Wt[i] * sumj * sin(theta);
+
+    }
+
+    answer = (v1b-v1a)/2.0*sumi;
+
+    // Normalize as in Eqn. (11). The factor 2 is due to the different theta integration limit (pi/2 instead of pi)
+    vol = 2.0*aa*bb+2.0*aa*cc+2.0*bb*cc;
+    answer = 2.0 * answer / Pi / vol / vol;
+
+    // Multiply by contrast^2
+    answer *= delrho*delrho;
+
+    //normalize by volume
+    answer *= vol;
+
+    //convert to [cm-1]
+    answer *= 1.0e8;
+
+    //Scale
+    answer *= scale;
+
+    // add in the background
+    answer += bkg;
+
+    return answer;
+}
+
+
+/*
+ * Hollow rectangular prism (a <= b <=c)
+ * using eqn. (13-15) from Nayuk & Huber, Z. Phys. Chem. 226, 837 (2012)
+ */
+ 
+double
+RectangularHollowPrism(double dp[], double q)
+{
+    int i, j;
+    double scale, aa, bb, cc, thickness, delrho, bkg;
+    int nordi=76;
+    int nordj=76;
+    double Pi;
+    double sumi, sumj;
+    double v1a, v1b, v2a, v2b;		//upper and lower integration limits
+    double answer;
+    double theta, phi, arg, termA1, termB1, termC1, termA2, termB2, termC2, AP1, AP2, AP, ahalf, bhalf, chalf, vol;
+    double sldp, sld ;
+
+    Pi = 4.0*atan(1.0);
+
+    scale = dp[0];
+    aa = dp[1];
+    bb = aa * dp[2];
+    cc = aa * dp[3];
+    thickness = dp[4];
+    sldp = dp[5];        // scattering length density of the object
+    sld = dp[6];         // scattering length density of the solvent
+    delrho = sldp - sld;
+    bkg = dp[7];
+
+    ahalf = 0.5 * aa;
+    bhalf = 0.5 * bb;
+    chalf = 0.5 * cc;
+
+    //Integration limits to use in Gaussian quadrature
+    v1a = 0.;
+    v1b = Pi/2;    //theta integration limits
+    v2a = 0.;
+    v2b = Pi/2.;   //phi integration limits
+
+    sumi = 0.0;
+    for(i=0;i<nordi;i++) {
+
+        theta = ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b )/2.0;	//z values for gaussian quadrature
+
+	    arg = q*chalf*cos(theta);
+	    if (fabs(arg) > 1.e-16) {termC1 = sin(arg) / arg;} else {termC1 = 1.0;}
+	    arg = q*(chalf-thickness)*cos(theta);
+	    if (fabs(arg) > 1.e-16) {termC2 = sin(arg) / arg;} else {termC2 = 1.0;}
+
+	    sumj = 0.0;
+	    for(j=0;j<nordj;j++) {
+
+            phi = ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b )/2.0;	//z values for gaussian quadrature
+
+            // Amplitude AP from eqn. (13), rewritten to avoid round-off effects when arg=0
+
+	        arg = q*ahalf*sin(theta)*sin(phi);
+	        if (fabs(arg) > 1.e-16) {termA1 = sin(arg) / arg;} else {termA1 = 1.0;}
+	        arg = q*(ahalf-thickness)*sin(theta)*sin(phi);
+	        if (fabs(arg) > 1.e-16) {termA2 = sin(arg) / arg;} else {termA2 = 1.0;}
+
+	        arg = q*bhalf*sin(theta)*cos(phi);
+	        if (fabs(arg) > 1.e-16) {termB1 = sin(arg) / arg;} else {termB1 = 1.0;}
+	        arg = q*(bhalf-thickness)*sin(theta)*cos(phi);
+	        if (fabs(arg) > 1.e-16) {termB2 = sin(arg) / arg;} else {termB2 = 1.0;}
+
+            AP1 = (aa*bb*cc) * termA1 * termB1 * termC1;
+            AP2 = 8 * (ahalf-thickness) * (bhalf-thickness) * (chalf-thickness) * termA2 * termB2 * termC2;
+            AP = AP1 - AP2;
+
+            sumj += Gauss76Wt[j] * (AP*AP);
+
+	    }
+
+	    sumj = (v2b-v2a) * sumj / 2.0;
+	    sumi += Gauss76Wt[i] * sumj * sin(theta);
+
+    }
+
+    answer = (v1b-v1a)/2.0*sumi;
+
+    // Normalize as in Eqn. (15).
+    // The factor 2 is due to the different theta integration limit (pi/2 instead of pi)
+    vol = (aa*bb*cc-((aa-2*thickness)*(bb-2*thickness)*(cc-2*thickness)));
+    answer = 2.0 * answer / Pi / vol / vol;
+
+    // Multiply by contrast^2
+    answer *= delrho*delrho;
+
+    //normalize by volume
+    answer *= vol;
+
+    //convert to [cm-1]
+    answer *= 1.0e8;
+
+    //Scale
+    answer *= scale;
+
+    // add in the background
+    answer += bkg;
+
+    return answer;
+}
+
+
Index: src/sans/models/c_extension/libigor/libCylinder.h
===================================================================
--- src/sans/models/c_extension/libigor/libCylinder.h	(revision 230f479e04002ef42bafef9a227d586459b717c1)
+++ src/sans/models/c_extension/libigor/libCylinder.h	(revision 201af9f481a69dc31baac2847a064fb4036d1223)
@@ -38,4 +38,7 @@
 double PolyCoreBicelle(double dp[], double q);
 double CSParallelepiped(double dp[], double q);
+double RectangularPrism(double dp[], double q);
+double RectangularHollowPrismInfThinWalls(double dp[], double q);
+double RectangularHollowPrism(double dp[], double q);
 
 /* internal functions */
Index: src/sans/models/include/RectangularHollowPrism.h
===================================================================
--- src/sans/models/include/RectangularHollowPrism.h	(revision 201af9f481a69dc31baac2847a064fb4036d1223)
+++ src/sans/models/include/RectangularHollowPrism.h	(revision 201af9f481a69dc31baac2847a064fb4036d1223)
@@ -0,0 +1,68 @@
+#if !defined(RectangularHollowPrism_h)
+#define RectangularHollowPrism_h
+#include "parameters.hh"
+
+/** Structure definition for RectangularHollowPrism parameters
+ * [PYTHONCLASS] = RectangularHollowPrismModel
+ * [DISP_PARAMS] = short_side, b2a_ratio, c2a_ratio, thickness
+   [DESCRIPTION] = <text> Form factor for a hollow rectangular prism with uniform scattering length density.
+		scale:Scale factor
+		short_side: shortest side of the rectangular prism  [A]
+		b2a_ratio: ratio b/a [adim]
+		c2a_ratio: ratio c/a [adim]
+		thickness: thickness of the walls [A]
+		sldPipe: Pipe_sld
+		sldSolv: solvent_sld
+		background:Incoherent Background [1/cm]
+		</text>
+	[FIXED]= <text> short_side.width; b2a_ratio.width; c2a_ratio.width; thicness.width; </text>
+
+ **/
+
+class RectangularHollowPrismModel{
+public:
+  // Model parameters
+  
+  /// Scale factor
+  //  [DEFAULT]=scale=1.0
+  Parameter scale;
+  
+  ///  Shortest side of rectangular prism [A]
+  //  [DEFAULT]=short_side=35 [A]
+  Parameter short_side;
+  
+  ///  Ratio b/a [adim]
+  //  [DEFAULT]=b2a_ratio=1 [adim]
+  Parameter b2a_ratio;
+  
+  ///  Ratio c/a [adim]
+  //  [DEFAULT]=c2a_ratio=1 [adim]
+  Parameter c2a_ratio;
+  
+  ///  Thickness of walls [A]
+  //  [DEFAULT]=thickness=1 [A]
+  Parameter thickness;
+  
+  /// SLD_Pipe [1/A^(2)]
+  //  [DEFAULT]=sldPipe=6.3e-6 [1/A^(2)]
+  Parameter sldPipe;
+  
+  /// sldSolv [1/A^(2)]
+  //  [DEFAULT]=sldSolv=1.0e-6 [1/A^(2)]
+  Parameter sldSolv;
+
+  /// Incoherent Background [1/cm]
+  //  [DEFAULT]=background=0.0 [1/cm]
+  Parameter background;
+  
+  // Constructor
+  RectangularHollowPrismModel();
+
+  // Operators to get I(Q)
+  double operator()(double q);
+  double operator()(double qx, double qy);
+  double calculate_ER();
+  double calculate_VR();
+  double evaluate_rphi(double q, double phi);
+};
+#endif
Index: src/sans/models/include/RectangularHollowPrismInfThinWalls.h
===================================================================
--- src/sans/models/include/RectangularHollowPrismInfThinWalls.h	(revision 201af9f481a69dc31baac2847a064fb4036d1223)
+++ src/sans/models/include/RectangularHollowPrismInfThinWalls.h	(revision 201af9f481a69dc31baac2847a064fb4036d1223)
@@ -0,0 +1,63 @@
+#if !defined(RectangularHollowPrismInfThinWalls_h)
+#define RectangularHollowPrismInfThinWalls_h
+#include "parameters.hh"
+
+/** Structure definition for RectangularHollowPrismInfThinWalls parameters
+ * [PYTHONCLASS] = RectangularHollowPrismInfThinWallsModel
+ * [DISP_PARAMS] = short_side, b2a_ratio, c2a_ratio
+   [DESCRIPTION] = <text> Form factor for a hollow rectangular prism with uniform scattering length density and infinitely thin walls.
+		scale:Scale factor
+		short_side: shortest side of the rectangular prism  [A]
+		b2a_ratio: ratio b/a [adim]
+		c2a_ratio: ratio c/a [adim]
+		sldPipe: Pipe_sld
+		sldSolv: solvent_sld
+		background:Incoherent Background [1/cm]
+		</text>
+	[FIXED]= <text> short_side.width; b2a_ratio.width; c2a_ratio.width; </text>
+
+ **/
+
+class RectangularHollowPrismInfThinWallsModel{
+public:
+  // Model parameters
+  
+  /// Scale factor
+  //  [DEFAULT]=scale=1.0
+  Parameter scale;
+  
+  ///  Shortest side of rectangular prism [A]
+  //  [DEFAULT]=short_side=35 [A]
+  Parameter short_side;
+  
+  ///  Ratio b/a [adim]
+  //  [DEFAULT]=b2a_ratio=1 [adim]
+  Parameter b2a_ratio;
+  
+  ///  Ratio c/a [adim]
+  //  [DEFAULT]=c2a_ratio=1 [adim]
+  Parameter c2a_ratio;
+  
+  /// SLD_Pipe [1/A^(2)]
+  //  [DEFAULT]=sldPipe=6.3e-6 [1/A^(2)]
+  Parameter sldPipe;
+  
+  /// sldSolv [1/A^(2)]
+  //  [DEFAULT]=sldSolv=1.0e-6 [1/A^(2)]
+  Parameter sldSolv;
+
+  /// Incoherent Background [1/cm]
+  //  [DEFAULT]=background=0.0 [1/cm]
+  Parameter background;
+  
+  // Constructor
+  RectangularHollowPrismInfThinWallsModel();
+
+  // Operators to get I(Q)
+  double operator()(double q);
+  double operator()(double qx, double qy);
+  double calculate_ER();
+  double calculate_VR();
+  double evaluate_rphi(double q, double phi);
+};
+#endif
Index: src/sans/models/include/RectangularPrism.h
===================================================================
--- src/sans/models/include/RectangularPrism.h	(revision 201af9f481a69dc31baac2847a064fb4036d1223)
+++ src/sans/models/include/RectangularPrism.h	(revision 201af9f481a69dc31baac2847a064fb4036d1223)
@@ -0,0 +1,63 @@
+#if !defined(RectangularPrism_h)
+#define RectangularPrism_h
+#include "parameters.hh"
+
+/** Structure definition for RectangularPrism parameters
+ * [PYTHONCLASS] = RectangularPrismModel
+ * [DISP_PARAMS] = short_side, b2a_ratio, c2a_ratio
+   [DESCRIPTION] = <text> Form factor for a massive rectangular prism with uniform scattering length density.
+		scale:Scale factor
+		short_side: shortest side of the rectangular prism  [A]
+		b2a_ratio: ratio b/a [adim]
+		c2a_ratio: ratio c/a [adim]
+		sldPipe: Pipe_sld
+		sldSolv: solvent_sld
+		background:Incoherent Background [1/cm]
+		</text>
+	[FIXED]= <text> short_side.width; b2a_ratio.width; c2a_ratio.width; </text>
+
+ **/
+
+class RectangularPrismModel{
+public:
+  // Model parameters
+  
+  /// Scale factor
+  //  [DEFAULT]=scale=1.0
+  Parameter scale;
+  
+  ///  Shortest side of rectangular prism [A]
+  //  [DEFAULT]=short_side=35 [A]
+  Parameter short_side;
+  
+  ///  Ratio b/a [adim]
+  //  [DEFAULT]=b2a_ratio=1 [adim]
+  Parameter b2a_ratio;
+  
+  ///  Ratio c/a [adim]
+  //  [DEFAULT]=c2a_ratio=1 [adim]
+  Parameter c2a_ratio;
+  
+  /// SLD_Pipe [1/A^(2)]
+  //  [DEFAULT]=sldPipe=6.3e-6 [1/A^(2)]
+  Parameter sldPipe;
+  
+  /// sldSolv [1/A^(2)]
+  //  [DEFAULT]=sldSolv=1.0e-6 [1/A^(2)]
+  Parameter sldSolv;
+
+  /// Incoherent Background [1/cm]
+  //  [DEFAULT]=background=0.0 [1/cm]
+  Parameter background;
+  
+  // Constructor
+  RectangularPrismModel();
+
+  // Operators to get I(Q)
+  double operator()(double q);
+  double operator()(double qx, double qy);
+  double calculate_ER();
+  double calculate_VR();
+  double evaluate_rphi(double q, double phi);
+};
+#endif
Index: src/sans/perspectives/fitting/models.py
===================================================================
--- src/sans/perspectives/fitting/models.py	(revision 27b7accfe21e1b3456e646ee01adff844993eae3)
+++ src/sans/perspectives/fitting/models.py	(revision 201af9f481a69dc31baac2847a064fb4036d1223)
@@ -677,4 +677,22 @@
         self.model_name_list.append(PringlesModel.__name__)
 
+        from sans.models.RectangularPrismModel import RectangularPrismModel
+        self.model_dictionary[RectangularPrismModel.__name__] = RectangularPrismModel
+        self.shape_list.append(RectangularPrismModel)
+        self.multiplication_factor.append(RectangularPrismModel)
+        self.model_name_list.append(RectangularPrismModel.__name__)
+
+        from sans.models.RectangularHollowPrismInfThinWallsModel import RectangularHollowPrismInfThinWallsModel
+        self.model_dictionary[RectangularHollowPrismInfThinWallsModel.__name__] = RectangularHollowPrismInfThinWallsModel
+        self.shape_list.append(RectangularHollowPrismInfThinWallsModel)
+        self.multiplication_factor.append(RectangularHollowPrismInfThinWallsModel)
+        self.model_name_list.append(RectangularHollowPrismInfThinWallsModel.__name__)
+
+        from sans.models.RectangularHollowPrismModel import RectangularHollowPrismModel
+        self.model_dictionary[RectangularHollowPrismModel.__name__] = RectangularHollowPrismModel
+        self.shape_list.append(RectangularHollowPrismModel)
+        self.multiplication_factor.append(RectangularHollowPrismModel)
+        self.model_name_list.append(RectangularHollowPrismModel.__name__)
+
         #from sans.models.FractalO_Z import FractalO_Z
         #self.model_dictionary[FractalO_Z.__name__] = FractalO_Z
