source: sasview/src/sas/models/c_extension/libigor/libCylinder.c @ 90662e6

/*      CylinderFit.c

A simplified project designed to act as a template for your curve fitting function.
The fitting function is a Cylinder form factor. No resolution effects are included (yet)
*/

#include "StandardHeaders.h"                    // Include ANSI headers, Mac headers
#include "GaussWeights.h"
#include "libCylinder.h"

/////////functions for WRC implementation of flexible cylinders
static double
gammaln(double xx)
{
  double x,y,tmp,ser;
  static double cof[6]={76.18009172947146,-86.50532032941677,
    24.01409824083091,-1.231739572450155,
    0.1208650973866179e-2,-0.5395239384953e-5};
  int j;

  y=x=xx;
  tmp=x+5.5;
  tmp -= (x+0.5)*log(tmp);
  ser=1.000000000190015;
  for (j=0;j<=5;j++) ser += cof[j]/++y;
  return -tmp+log(2.5066282746310005*ser/x);
}
/////////functions for WRC implementation of flexible cylinders

//
static double
AlphaSquare(double x)
{
    double yy;
    yy = pow( (1.0 + (x/3.12)*(x/3.12) + (x/8.67)*(x/8.67)*(x/8.67)),(0.176/3.0) );

    return (yy);
}

//
static double
Rgsquarezero(double q, double L, double b)
{
    double yy;
    yy = (L*b/6.0) * (1.0 - 1.5*(b/L) + 1.5*pow((b/L),2) - 0.75*pow((b/L),3)*(1.0 - exp(-2.0*(L/b))));

    return (yy);
}

//
static double
Rgsquareshort(double q, double L, double b)
{
    double yy;
    yy = AlphaSquare(L/b) * Rgsquarezero(q,L,b);

    return (yy);
}

//
static double
Rgsquare(double q, double L, double b)
{
    double yy;
    yy = AlphaSquare(L/b)*L*b/6.0;

    return (yy);
}

// ?? funciton is not used - but should the log actually be log10???
/*
static double
miu(double x)
{
    double yy;
    yy = (1.0/8.0)*(9.0*x - 2.0 + 2.0*log(1.0 + x)/x)*exp(1.0/2.565*(1.0/x + (1.0 - 1.0/(x*x))*log(1.0 + x)));

    return (yy);
}
*/

//WR named this w (too generic)
static double
w_WR(double x)
{
    double yy;
    yy = 0.5*(1 + tanh((x - 1.523)/0.1477));

    return (yy);
}

//
static double
u1(double q, double L, double b)
{
    double yy;

    yy = Rgsquareshort(q,L,b)*q*q;

    return (yy);
}

// was named u
static double
u_WR(double q, double L, double b)
{
    double yy;
    yy = Rgsquare(q,L,b)*q*q;
    return (yy);
}


//
static double
Sdebye1(double q, double L, double b)
{
    double yy;
    yy = 2.0*(exp(-u1(q,L,b)) + u1(q,L,b) -1.0)/( pow((u1(q,L,b)),2.0) );

    return (yy);
}


//
static double
Sdebye(double q, double L, double b)
{
    double yy;
    yy = 2.0*(exp(-u_WR(q,L,b)) + u_WR(q,L,b) -1.0)/(pow((u_WR(q,L,b)),2));

    return (yy);
}

//
static double
Sexv(double q, double L, double b)
{
    double yy,C1,C2,C3,miu,Rg2;
    C1=1.22;
    C2=0.4288;
    C3=-1.651;
    miu = 0.585;

    Rg2 = Rgsquare(q,L,b);

    yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) +  C2*pow((q*sqrt(Rg2)),(-2.0/miu)) +    C3*pow((q*sqrt(Rg2)),(-3.0/miu)));

    return (yy);
}

// this must be WR modified version
static double
Sexvnew(double q, double L, double b)
{
    double yy,C1,C2,C3,miu;
    double del=1.05,C_star2,Rg2;

    C1=1.22;
    C2=0.4288;
    C3=-1.651;
    miu = 0.585;

    //calculating the derivative to decide on the corection (cutoff) term?
    // I have modified this from WRs original code

    if( (Sexv(q*del,L,b)-Sexv(q,L,b))/(q*del - q) >= 0.0 ) {
        C_star2 = 0.0;
    } else {
        C_star2 = 1.0;
    }

    Rg2 = Rgsquare(q,L,b);

  yy = (1.0 - w_WR(q*sqrt(Rg2)))*Sdebye(q,L,b) + C_star2*w_WR(q*sqrt(Rg2))*(C1*pow((q*sqrt(Rg2)),(-1.0/miu)) + C2*pow((q*sqrt(Rg2)),(-2.0/miu)) + C3*pow((q*sqrt(Rg2)),(-3.0/miu)));

    return (yy);
}

// these are the messy ones
static double
a2short(double q, double L, double b, double p1short, double p2short, double q0)
{
    double yy,Rg2_sh;
    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p;
    double pi;

    E = 2.718281828459045091;
  pi = 4.0*atan(1.0);
    Rg2_sh = Rgsquareshort(q,L,b);
    Rg2_sh2 = Rg2_sh*Rg2_sh;
    t1 = ((q0*q0*Rg2_sh)/(b*b));
    Et1 = pow(E,t1);
    Emt1 =pow(E,-t1);
    q02 = q0*q0;
    q0p = pow(q0,(-4.0 + p2short) );

    //E is the number e
    yy = ((-(1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b*b*b*L - 8.0*b*b*b*Et1*L - 2.0*b*b*b*L*p1short + 2.0*b*b*b*Et1*L*p1short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p1short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p1short*pi*q02*q0*Rg2_sh2)))))));

    return (yy);
}

//
static double
a1short(double q, double L, double b, double p1short, double p2short, double q0)
{
    double yy,Rg2_sh;
    double t1,E,Rg2_sh2,Et1,Emt1,q02,q0p,b3;
  double pi;

    E = 2.718281828459045091;
  pi = 4.0*atan(1.0);
    Rg2_sh = Rgsquareshort(q,L,b);
    Rg2_sh2 = Rg2_sh*Rg2_sh;
    b3 = b*b*b;
    t1 = ((q0*q0*Rg2_sh)/(b*b));
    Et1 = pow(E,t1);
    Emt1 =pow(E,-t1);
    q02 = q0*q0;
    q0p = pow(q0,(-4.0 + p1short) );

    yy = ((1.0/(L*((p1short - p2short))*Rg2_sh2)*((b*Emt1*q0p*((8.0*b3*L - 8.0*b3*Et1*L - 2.0*b3*L*p2short + 2.0*b3*Et1*L*p2short + 4.0*b*L*q02*Rg2_sh + 4.0*b*Et1*L*q02*Rg2_sh - 2.0*b*Et1*L*p2short*q02*Rg2_sh - Et1*pi*q02*q0*Rg2_sh2 + Et1*p2short*pi*q02*q0*Rg2_sh2))))));

    return(yy);
}


//need to define this on my own
static double
sech_WR(double x)
{
  return(1/cosh(x));
}

// this one will be lots of trouble
static double
a2long(double q, double L, double b, double p1, double p2, double q0)
{
    double yy,C1,C2,C3,C4,C5,miu,C,Rg2;
    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,pi;
    double E,b2,b3,b4,q02,q03,q04,q05,Rg22;

    pi = 4.0*atan(1.0);
    E = 2.718281828459045091;
    if( L/b > 10.0) {
        C = 3.06/pow((L/b),0.44);
    } else {
        C = 1.0;
    }

    C1 = 1.22;
    C2 = 0.4288;
    C3 = -1.651;
    C4 = 1.523;
    C5 = 0.1477;
    miu = 0.585;

    Rg2 = Rgsquare(q,L,b);
    Rg22 = Rg2*Rg2;
    b2 = b*b;
    b3 = b*b*b;
    b4 = b3*b;
    q02 = q0*q0;
    q03 = q0*q0*q0;
    q04 = q03*q0;
    q05 = q04*q0;

    t1 = (1.0/(b* p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2)) ));

    t2 = (b*C*(((-1.0*((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7*b2)/(15.0*q02*Rg2)))*Rg2)/b)))/L;

    t3 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2.0))/(2.0*C5);

    t4 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2))/(C5*q04*Rg22);

    t5 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);

    t6 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q05*Rg22);

    t7 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu));

    t8 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));

    t9 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))))/L;

    t10 = (2.0*b4*(((-1) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))/(q04*Rg22);


    yy = ((-1.0*(t1* ((-pow(q0,-p1)*(((b2*pi)/(L*q02) + t2 + t3 - t4 + t5 - t6 + 1.0/2.0*t7*t8)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) + t9 + t10 + 1.0/2.0*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))))))));

    return (yy);
}
//
static double
a1long(double q, double L, double b, double p1, double p2, double q0)
{
    double yy,C,C1,C2,C3,C4,C5,miu,Rg2;
    double t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15;
    double E,pi;
    double b2,b3,b4,q02,q03,q04,q05,Rg22;

    pi = 4.0*atan(1.0);
    E = 2.718281828459045091;

    if( L/b > 10.0) {
        C = 3.06/pow((L/b),0.44);
    } else {
        C = 1.0;
    }

    C1 = 1.22;
    C2 = 0.4288;
    C3 = -1.651;
    C4 = 1.523;
    C5 = 0.1477;
    miu = 0.585;

    Rg2 = Rgsquare(q,L,b);
    Rg22 = Rg2*Rg2;
    b2 = b*b;
    b3 = b*b*b;
    b4 = b3*b;
    q02 = q0*q0;
    q03 = q0*q0*q0;
    q04 = q03*q0;
    q05 = q04*q0;

    t1 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2))) + (7.0*b2)/(15.0*q02*Rg2))));

    t2 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));

    t3 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))));

    t4 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));

    t5 = (1.0/(b*p1*pow(q0,((-1.0) - p1 - p2)) - b*p2*pow(q0,((-1.0) - p1 - p2))));

    t6 = (b*C*(((-((14.0*b3)/(15.0*q03*Rg2))) + (14.0*b3*pow(E,(-((q02*Rg2)/b2))))/(15.0*q03*Rg2) + (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*((11.0/15.0 + (7.0*b2)/(15.0*q02*Rg2)))*Rg2)/b)));

    t7 = (sqrt(Rg2)*((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));

    t8 = (b4*sqrt(Rg2)*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*pow(sech_WR(((-C4) + (sqrt(Rg2)*q0)/b)/C5),2));

    t9 = (2.0*b4*(((2.0*q0*Rg2)/b - (2.0*pow(E,(-((q02*Rg2)/b2)))*q0*Rg2)/b))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));

    t10 = (8.0*b4*b*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));

    t11 = (((-((3.0*C3*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 3.0/miu)))/miu)) - (2.0*C2*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 2.0/miu)))/miu - (C1*sqrt(Rg2)*pow((((sqrt(Rg2)*q0)/b)),((-1.0) - 1.0/miu)))/miu));

    t12 = ((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)));

    t13 = (b*C*((4.0/15.0 - pow(E,(-((q02*Rg2)/b2)))*((11.0/15.0 + (7.0*b2)/(15.0*q02* Rg2))) + (7.0*b2)/(15.0*q02*Rg2))));

    t14 = (2.0*b4*(((-1.0) + pow(E,(-((q02*Rg2)/b2))) + (q02*Rg2)/b2))*((1.0 + 1.0/2.0*(((-1.0) - tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5))))));

    t15 = ((C3*pow((((sqrt(Rg2)*q0)/b)),((-3.0)/miu)) + C2*pow((((sqrt(Rg2)*q0)/b)),((-2.0)/miu)) + C1*pow((((sqrt(Rg2)*q0)/b)),((-1.0)/miu))));


    yy = (pow(q0,p1)*(((-((b*pi)/(L*q0))) +t1/L +t2/(q04*Rg22) + 1.0/2.0*t3*t4)) + (t5*((pow(q0,(p1 - p2))*(((-pow(q0,(-p1)))*(((b2*pi)/(L*q02) +t6/L +t7/(2.0*C5) -t8/(C5*q04*Rg22) +t9/(q04*Rg22) -t10/(q05*Rg22) + 1.0/2.0*t11*t12)) - b*p1*pow(q0,((-1.0) - p1))*(((-((b*pi)/(L*q0))) +t13/L +t14/(q04*Rg22) + 1.0/2.0*t15*((1.0 + tanh(((-C4) + (sqrt(Rg2)*q0)/b)/C5)))))))))));

    return (yy);
}


static double
Sk_WR(double q, double L, double b)
{
  //
  double p1,p2,p1short,p2short,q0;
  double C,ans,q0short,Sexvmodify,pi;

  pi = 4.0*atan(1.0);

  p1 = 4.12;
  p2 = 4.42;
  p1short = 5.36;
  p2short = 5.62;
  q0 = 3.1;
  //
  q0short = fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0);

  //
  if(L/b > 10.0) {
    C = 3.06/pow((L/b),0.44);
  } else {
    C = 1.0;
  }
  //

  if( L > 4*b ) { // Longer Chains
    if (q*b <= 3.1) {   //Modified by Yun on Oct. 15,
      Sexvmodify = Sexvnew(q, L, b);
      ans = Sexvmodify + C * (4.0/15.0 + 7.0/(15.0*u_WR(q,L,b)) - (11.0/15.0 + 7.0/(15.0*u_WR(q,L,b)))*exp(-u_WR(q,L,b)))*(b/L);
    } else { //q(i)*b > 3.1
      ans = a1long(q, L, b, p1, p2, q0)/(pow((q*b),p1)) + a2long(q, L, b, p1, p2, q0)/(pow((q*b),p2)) + pi/(q*L);
    }
  } else { //L <= 4*b Shorter Chains
    if (q*b <= fmax(1.9/sqrt(Rgsquareshort(q,L,b)),3.0) ) {
      if (q*b<=0.01) {
        ans = 1.0 - Rgsquareshort(q,L,b)*(q*q)/3.0;
      } else {
        ans = Sdebye1(q,L,b);
      }
    } else {  //q*b > max(1.9/sqrt(Rgsquareshort(q(i),L,b)),3)
      ans = a1short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p1short)) + a2short(q,L,b,p1short,p2short,q0short)/(pow((q*b),p2short)) + pi/(q*L);
    }
  }

  return(ans);
  //return(a2long(q, L, b, p1, p2, q0));
}

/*      CylinderForm  :  calculates the form factor of a cylinder at the give x-value p->x

Warning:
The call to WaveData() below returns a pointer to the middle
of an unlocked Macintosh handle. In the unlikely event that your
calculations could cause memory to move, you should copy the coefficient
values to local variables or an array before such operations.
*/
double
CylinderForm(double dp[], double q)
{

        int i;
        double Pi;
        double scale,radius,length,delrho,bkg,halfheight,sldCyl,sldSolv;                //local variables of coefficient wave
        int nord=76;                    //order of integration
        double uplim,lolim;             //upper and lower integration limits
        double summ,zi,yyy,answer,vcyl;                 //running tally of integration
        
        Pi = 4.0*atan(1.0);
        lolim = 0.0;
        uplim = Pi/2.0;
        
        summ = 0.0;                     //initialize intergral
        
        scale = dp[0];                  //make local copies in case memory moves
        radius = dp[1];
        length = dp[2];
        sldCyl = dp[3];
        sldSolv = dp[4];
        bkg = dp[5];
        
        delrho = sldCyl-sldSolv;
        halfheight = length/2.0;
        for(i=0;i<nord;i++) {
                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                yyy = Gauss76Wt[i] * CylKernel(q, radius, halfheight, zi);
                summ += yyy;
        }
        
        answer = (uplim-lolim)/2.0*summ;
        // Multiply by contrast^2
        answer *= delrho*delrho;
        //normalize by cylinder volume
        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
        vcyl=Pi*radius*radius*length;
        answer *= vcyl;
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
        
        return answer;
}

/*      EllipCyl76X  :  calculates the form factor of a elliptical cylinder at the given x-value p->x

Uses 76 pt Gaussian quadrature for both integrals

Warning:
The call to WaveData() below returns a pointer to the middle
of an unlocked Macintosh handle. In the unlikely event that your
calculations could cause memory to move, you should copy the coefficient
values to local variables or an array before such operations.
*/
double
EllipCyl76(double dp[], double q)
{
        int i,j;
        double Pi,slde,sld;
        double scale,ra,nu,length,delrho,bkg;           //local variables of coefficient wave
        int nord=76;                    //order of integration
        double va,vb;           //upper and lower integration limits
        double summ,zi,yyy,answer,vell;                 //running tally of integration
        double summj,vaj,vbj,zij,arg, si;                       //for the inner integration

        Pi = 4.0*atan(1.0);
        va = 0.0;
        vb = 1.0;               //orintational average, outer integral
        vaj=0.0;
        vbj=Pi;         //endpoints of inner integral
        
        summ = 0.0;                     //initialize intergral
        
        scale = dp[0];                  //make local copies in case memory moves
        ra = dp[1];
        nu = dp[2];
        length = dp[3];
        slde = dp[4];
        sld = dp[5];
        delrho = slde - sld;
        bkg = dp[6];
        
        for(i=0;i<nord;i++) {
                //setup inner integral over the ellipsoidal cross-section
                summj=0;
                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy
                arg = ra*sqrt(1.0-zi*zi);
                for(j=0;j<nord;j++) {
                        //76 gauss points for the inner integral as well
                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "y" dummy
                        yyy = Gauss76Wt[j] * EllipCylKernel(q,arg,nu,zij);
                        summj += yyy;
                }
                //now calculate the value of the inner integral
                answer = (vbj-vaj)/2.0*summj;
                //divide integral by Pi
                answer /=Pi;
                
                //now calculate outer integral
                arg = q*length*zi/2.0;
                if (arg == 0.0){
                        si = 1.0;
                }else{
                        si = sin(arg) * sin(arg) / arg / arg;
                }
                yyy = Gauss76Wt[i] * answer * si;
                summ += yyy;
        }
        answer = (vb-va)/2.0*summ;
        // Multiply by contrast^2
        answer *= delrho*delrho;
        //normalize by cylinder volume
        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
        vell = Pi*ra*(nu*ra)*length;
        answer *= vell;
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
        
        return answer;
}

/*      EllipCyl20X  :  calculates the form factor of a elliptical cylinder at the given x-value p->x

Uses 76 pt Gaussian quadrature for orientational integral
Uses 20 pt quadrature for the inner integral over the elliptical cross-section

Warning:
The call to WaveData() below returns a pointer to the middle
of an unlocked Macintosh handle. In the unlikely event that your
calculations could cause memory to move, you should copy the coefficient
values to local variables or an array before such operations.
*/
double
EllipCyl20(double dp[], double q)
{
        int i,j;
        double Pi,slde,sld;
        double scale,ra,nu,length,delrho,bkg;           //local variables of coefficient wave
        int nordi=76;                   //order of integration
        int nordj=20;
        double va,vb;           //upper and lower integration limits
        double summ,zi,yyy,answer,vell;                 //running tally of integration
        double summj,vaj,vbj,zij,arg,si;                        //for the inner integration

        Pi = 4.0*atan(1.0);
        va = 0.0;
        vb = 1.0;               //orintational average, outer integral
        vaj=0.0;
        vbj=Pi;         //endpoints of inner integral
        
        summ = 0.0;                     //initialize intergral
        
        scale = dp[0];                  //make local copies in case memory moves
        ra = dp[1];
        nu = dp[2];
        length = dp[3];
        slde = dp[4];
        sld = dp[5];
        delrho = slde - sld;
        bkg = dp[6];
        
        for(i=0;i<nordi;i++) {
                //setup inner integral over the ellipsoidal cross-section
                summj=0;
                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy
                arg = ra*sqrt(1.0-zi*zi);
                for(j=0;j<nordj;j++) {
                        //20 gauss points for the inner integral
                        zij = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "y" dummy
                        yyy = Gauss20Wt[j] * EllipCylKernel(q,arg,nu,zij);
                        summj += yyy;
                }
                //now calculate the value of the inner integral
                answer = (vbj-vaj)/2.0*summj;
                //divide integral by Pi
                answer /=Pi;
                
                //now calculate outer integral
                arg = q*length*zi/2.0;
                if (arg == 0.0){
                        si = 1.0;
                }else{
                        si = sin(arg) * sin(arg) / arg / arg;
                }
                yyy = Gauss76Wt[i] * answer * si;
                summ += yyy;
        }
        
        answer = (vb-va)/2.0*summ;
        // Multiply by contrast^2
        answer *= delrho*delrho;
        //normalize by cylinder volume
        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
        vell = Pi*ra*(nu*ra)*length;
        answer *= vell;
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;

        return answer;
}

/*      TriaxialEllipsoidX  :  calculates the form factor of a Triaxial Ellipsoid at the given x-value p->x

Uses 76 pt Gaussian quadrature for both integrals

Warning:
The call to WaveData() below returns a pointer to the middle
of an unlocked Macintosh handle. In the unlikely event that your
calculations could cause memory to move, you should copy the coefficient
values to local variables or an array before such operations.
*/
double
TriaxialEllipsoid(double dp[], double q)
{
        int i,j;
        double Pi;
        double scale,aa,bb,cc,delrho,bkg;               //local variables of coefficient wave
        int nordi=76;                   //order of integration
        int nordj=76;
        double va,vb;           //upper and lower integration limits
        double summ,zi,yyy,answer;                      //running tally of integration
        double summj,vaj,vbj,zij,slde,sld;                      //for the inner integration
        
        Pi = 4.0*atan(1.0);
        va = 0.0;
        vb = 1.0;               //orintational average, outer integral
        vaj = 0.0;
        vbj = 1.0;              //endpoints of inner integral
        
        summ = 0.0;                     //initialize intergral
        
        scale = dp[0];                  //make local copies in case memory moves
        aa = dp[1];
        bb = dp[2];
        cc = dp[3];
        slde = dp[4];
        sld = dp[5];
        delrho = slde - sld;
        bkg = dp[6];
        for(i=0;i<nordi;i++) {
                //setup inner integral over the ellipsoidal cross-section
                summj=0.0;
                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "x" dummy
                for(j=0;j<nordj;j++) {
                        //20 gauss points for the inner integral
                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "y" dummy
                        yyy = Gauss76Wt[j] * TriaxialKernel(q,aa,bb,cc,zi,zij);
                        summj += yyy;
                }
                //now calculate the value of the inner integral
                answer = (vbj-vaj)/2.0*summj;
                
                //now calculate outer integral
                yyy = Gauss76Wt[i] * answer;
                summ += yyy;
        }               //final scaling is done at the end of the function, after the NT_FP64 case
        
        answer = (vb-va)/2.0*summ;
        // Multiply by contrast^2
        answer *= delrho*delrho;
        //normalize by ellipsoid volume
        answer *= 4.0*Pi/3.0*aa*bb*cc;
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
        
        return answer;
}

/*      ParallelepipedX  :  calculates the form factor of a Parallelepiped (a rectangular solid)
at the given x-value p->x

Uses 76 pt Gaussian quadrature for both integrals

Warning:
The call to WaveData() below returns a pointer to the middle
of an unlocked Macintosh handle. In the unlikely event that your
calculations could cause memory to move, you should copy the coefficient
values to local variables or an array before such operations.
*/
double
Parallelepiped(double dp[], double q)
{
        int i,j;
        double scale,aa,bb,cc,delrho,bkg;               //local variables of coefficient wave
        int nordi=76;                   //order of integration
        int nordj=76;
        double va,vb;           //upper and lower integration limits
        double summ,yyy,answer;                 //running tally of integration
        double summj,vaj,vbj;                   //for the inner integration
        double mu,mudum,arg,sigma,uu,vol,sldp,sld;
        
        
        //      Pi = 4.0*atan(1.0);
        va = 0.0;
        vb = 1.0;               //orintational average, outer integral
        vaj = 0.0;
        vbj = 1.0;              //endpoints of inner integral

        summ = 0.0;                     //initialize intergral
        
        scale = dp[0];                  //make local copies in case memory moves
        aa = dp[1];
        bb = dp[2];
        cc = dp[3];
        sldp = dp[4];
        sld = dp[5];
        delrho = sldp - sld;
        bkg = dp[6];
        
        mu = q*bb;
        vol = aa*bb*cc;
        // normalize all WRT bb
        aa = aa/bb;
        cc = cc/bb;
        
        for(i=0;i<nordi;i++) {
                //setup inner integral over the ellipsoidal cross-section
                summj=0.0;
                sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;          //the outer dummy
                
                for(j=0;j<nordj;j++) {
                        //76 gauss points for the inner integral
                        uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;         //the inner dummy
                        mudum = mu*sqrt(1.0-sigma*sigma);
                        yyy = Gauss76Wt[j] * PPKernel(aa,mudum,uu);
                        summj += yyy;
                }
                //now calculate the value of the inner integral
                answer = (vbj-vaj)/2.0*summj;

                arg = mu*cc*sigma/2.0;
                if ( arg == 0.0 ) {
                        answer *= 1.0;
                } else {
                        answer *= sin(arg)*sin(arg)/arg/arg;
                }
                
                //now sum up the outer integral
                yyy = Gauss76Wt[i] * answer;
                summ += yyy;
        }               //final scaling is done at the end of the function, after the NT_FP64 case
        
        answer = (vb-va)/2.0*summ;
        // 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;
}

/*      HollowCylinderX  :  calculates the form factor of a Hollow Cylinder
at the given x-value p->x

Uses 76 pt Gaussian quadrature for the single integral

Warning:
The call to WaveData() below returns a pointer to the middle
of an unlocked Macintosh handle. In the unlikely event that your
calculations could cause memory to move, you should copy the coefficient
values to local variables or an array before such operations.
*/
double
HollowCylinder(double dp[], double q)
{
        int i;
        double scale,rcore,rshell,length,delrho,bkg;            //local variables of coefficient wave
        int nord=76;                    //order of integration
        double va,vb,zi;                //upper and lower integration limits
        double summ,answer,pi,sldc,sld;                 //running tally of integration
        
        pi = 4.0*atan(1.0);
        va = 0.0;
        vb = 1.0;               //limits of numerical integral

        summ = 0.0;                     //initialize intergral
        
        scale = dp[0];          //make local copies in case memory moves
        rcore = dp[1];
        rshell = dp[2];
        length = dp[3];
        sldc = dp[4];
        sld = dp[5];
        delrho = sldc - sld;
        bkg = dp[6];
        
        for(i=0;i<nord;i++) {
                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
                summ += Gauss76Wt[i] * HolCylKernel(q, rcore, rshell, length, zi);
        }
        
        answer = (vb-va)/2.0*summ;
        // Multiply by contrast^2
        answer *= delrho*delrho;
        //normalize by volume
        answer *= pi*(rshell*rshell-rcore*rcore)*length;
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
        
        return answer;
}

/*      EllipsoidFormX  :  calculates the form factor of an ellipsoid of revolution with semiaxes a:a:nua
at the given x-value p->x

Uses 76 pt Gaussian quadrature for the single integral

Warning:
The call to WaveData() below returns a pointer to the middle
of an unlocked Macintosh handle. In the unlikely event that your
calculations could cause memory to move, you should copy the coefficient
values to local variables or an array before such operations.
*/
double
EllipsoidForm(double dp[], double q)
{
        int i;
        double scale,a,nua,delrho,bkg;          //local variables of coefficient wave
        int nord=76;                    //order of integration
        double va,vb,zi;                //upper and lower integration limits
        double summ,answer,pi,slde,sld;                 //running tally of integration
        
        pi = 4.0*atan(1.0);
        va = 0.0;
        vb = 1.0;               //limits of numerical integral

        summ = 0.0;                     //initialize intergral
        
        scale = dp[0];                  //make local copies in case memory moves
        nua = dp[1];
        a = dp[2];
        slde = dp[3];
        sld = dp[4];
        delrho = slde - sld;
        bkg = dp[5];
        
        for(i=0;i<nord;i++) {
                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;
                summ += Gauss76Wt[i] * EllipsoidKernel(q, a, nua, zi);
        }
        
        answer = (vb-va)/2.0*summ;
        // Multiply by contrast^2
        answer *= delrho*delrho;
        //normalize by volume
        answer *= 4.0*pi/3.0*a*a*nua;
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
        
        return answer;
}


/*      Cyl_PolyRadiusX  :  calculates the form factor of a cylinder at the given x-value p->x
the cylinder has a polydisperse cross section

*/
double
Cyl_PolyRadius(double dp[], double q)
{
        int i;
        double scale,radius,length,pd,delrho,bkg;               //local variables of coefficient wave
        int nord=20;                    //order of integration
        double uplim,lolim;             //upper and lower integration limits
        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
        double range,zz,Pi,sldc,sld;
        
        Pi = 4.0*atan(1.0);
        range = 3.4;
        
        summ = 0.0;                     //initialize intergral
        
        scale = dp[0];                  //make local copies in case memory moves
        radius = dp[1];
        length = dp[2];
        pd = dp[3];
        sldc = dp[4];
        sld = dp[5];
        delrho = sldc - sld;
        bkg = dp[6];
        
        zz = (1.0/pd)*(1.0/pd) - 1.0;
        
        lolim = radius*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
        if(lolim<0.0) {
                lolim = 0.0;
        }
        if(pd>0.3) {
                range = 3.4 + (pd-0.3)*18.0;
        }
        uplim = radius*(1.0+range*pd);
        
        for(i=0;i<nord;i++) {
                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                yyy = Gauss20Wt[i] * Cyl_PolyRadKernel(q, radius, length, zz, delrho, zi);
                summ += yyy;
        }
        
        answer = (uplim-lolim)/2.0*summ;
        //normalize by average cylinder volume
        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
        Vpoly=Pi*radius*radius*length*(zz+2.0)/(zz+1.0);
        answer /= Vpoly;
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
        
        return answer;
}

/*      Cyl_PolyLengthX  :  calculates the form factor of a cylinder at the given x-value p->x
the cylinder has a polydisperse Length

*/
double
Cyl_PolyLength(double dp[], double q)
{
        int i;
        double scale,radius,length,pd,delrho,bkg;               //local variables of coefficient wave
        int nord=20;                    //order of integration
        double uplim,lolim;             //upper and lower integration limits
        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
        double range,zz,Pi,sldc,sld;
        
        
        Pi = 4.0*atan(1.0);
        range = 3.4;
        
        summ = 0.0;                     //initialize intergral
        
        scale = dp[0];                  //make local copies in case memory moves
        radius = dp[1];
        length = dp[2];
        pd = dp[3];
        sldc = dp[4];
        sld = dp[5];
        delrho = sldc - sld;
        bkg = dp[6];
        
        zz = (1.0/pd)*(1.0/pd) - 1.0;
        
        lolim = length*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
        if(lolim<0.0) {
                lolim = 0.0;
        }
        if(pd>0.3) {
                range = 3.4 + (pd-0.3)*18.0;
        }
        uplim = length*(1.0+range*pd);
        
        for(i=0;i<nord;i++) {
                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                yyy = Gauss20Wt[i] * Cyl_PolyLenKernel(q, radius, length, zz, delrho, zi);
                summ += yyy;
        }
        
        answer = (uplim-lolim)/2.0*summ;
        //normalize by average cylinder volume (first moment)
        //NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl
        Vpoly=Pi*radius*radius*length;
        answer /= Vpoly;
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
        
        return answer;
}

/*      CoreShellCylinderX  :  calculates the form factor of a cylinder at the given x-value p->x
the cylinder has a core-shell structure

*/
double
CoreShellCylinder(double dp[], double q)
{
        int i;
        double scale,rcore,length,bkg;          //local variables of coefficient wave
        double thick,rhoc,rhos,rhosolv;
        int nord=76;                    //order of integration
        double uplim,lolim,halfheight;          //upper and lower integration limits
        double summ,zi,yyy,answer,Vcyl;                 //running tally of integration
        double Pi;
        
        Pi = 4.0*atan(1.0);
        
        lolim = 0.0;
        uplim = Pi/2.0;
        
        summ = 0.0;                     //initialize intergral
        
        scale = dp[0];          //make local copies in case memory moves
        rcore = dp[1];
        thick = dp[2];
        length = dp[3];
        rhoc = dp[4];
        rhos = dp[5];
        rhosolv = dp[6];
        bkg = dp[7];
        
        halfheight = length/2.0;
        
        for(i=0;i<nord;i++) {
                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                yyy = Gauss76Wt[i] * CoreShellCylKernel(q, rcore, thick, rhoc,rhos,rhosolv, halfheight, zi);
                summ += yyy;
        }
        
        answer = (uplim-lolim)/2.0*summ;
        // length is the total core length 
        Vcyl=Pi*(rcore+thick)*(rcore+thick)*(length+2.0*thick);
        answer /= Vcyl;
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
        
        return answer;
}


/*      PolyCoShCylinderX  :  calculates the form factor of a core-shell cylinder at the given x-value p->x
the cylinder has a polydisperse CORE radius

*/
double
PolyCoShCylinder(double dp[], double q)
{
        int i;
        double scale,radius,length,sigma,bkg;           //local variables of coefficient wave
        double rad,radthick,facthick,rhoc,rhos,rhosolv;
        int nord=20;                    //order of integration
        double uplim,lolim;             //upper and lower integration limits
        double summ,yyy,answer,Vpoly;                   //running tally of integration
        double Pi,AR,Rsqrsumm,Rsqryyy,Rsqr;
                
        Pi = 4.0*atan(1.0);
        
        summ = 0.0;                     //initialize intergral
        Rsqrsumm = 0.0;
        
        scale = dp[0];
        radius = dp[1];
        sigma = dp[2];                          //sigma is the standard mean deviation
        length = dp[3];
        radthick = dp[4];
        facthick= dp[5];
        rhoc = dp[6];
        rhos = dp[7];
        rhosolv = dp[8];
        bkg = dp[9];
        
        lolim = exp(log(radius)-(4.*sigma));
        if (lolim<0.0) {
                lolim=0.0;              //to avoid numerical error when  va<0 (-ve r value)
        }
        uplim = exp(log(radius)+(4.*sigma));
        
        for(i=0;i<nord;i++) {
                rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                AR=(1.0/(rad*sigma*sqrt(2.0*Pi)))*exp(-(0.5*((log(radius/rad))/sigma)*((log(radius/rad))/sigma)));
                yyy = AR* Gauss20Wt[i] * CSCylIntegration(q,rad,radthick,facthick,rhoc,rhos,rhosolv,length);
                Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick);             //SRK normalize to total dimensions
                summ += yyy;
                Rsqrsumm += Rsqryyy;
        }
        
        answer = (uplim-lolim)/2.0*summ;
        Rsqr = (uplim-lolim)/2.0*Rsqrsumm;
        //normalize by average cylinder volume
        Vpoly = Pi*Rsqr*(length+2*facthick);
        answer /= Vpoly;
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
        
        return answer;
}

/*      OblateFormX  :  calculates the form factor of a core-shell Oblate ellipsoid at the given x-value p->x
the ellipsoid has a core-shell structure

*/
double
OblateForm(double dp[], double q)
{
        int i;
        double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg;
        int nord=76;                    //order of integration
        double uplim,lolim;             //upper and lower integration limits
        double summ,zi,yyy,answer,oblatevol;                    //running tally of integration
        double Pi,sldc,slds,sld;
        
        Pi = 4.0*atan(1.0);
        
        lolim = 0.0;
        uplim = 1.0;
        
        summ = 0.0;                     //initialize intergral
        
        
        scale = dp[0];          //make local copies in case memory moves
        crmaj = dp[1];
        crmin = dp[2];
        trmaj = dp[3];
        trmin = dp[4];
        sldc = dp[5];
        slds = dp[6];
        sld = dp[7];
        delpc = sldc - slds;    //core - shell
        delps = slds - sld;             //shell - solvent
        bkg = dp[8];

        for(i=0;i<nord;i++) {
                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                yyy = Gauss76Wt[i] * gfn4(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q);
                summ += yyy;
        }
        
        answer = (uplim-lolim)/2.0*summ;
        // normalize by particle volume
        oblatevol = 4.0*Pi/3.0*trmaj*trmaj*trmin;
        answer /= oblatevol;
        
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
        
        return answer;
}

/*      ProlateFormX  :  calculates the form factor of a core-shell Prolate ellipsoid at the given x-value p->x
the ellipsoid has a core-shell structure

*/
double
ProlateForm(double dp[], double q)
{
        int i;
        double scale,crmaj,crmin,trmaj,trmin,delpc,delps,bkg;
        int nord=76;                    //order of integration
        double uplim,lolim;             //upper and lower integration limits
        double summ,zi,yyy,answer,prolatevol;                   //running tally of integration
        double Pi,sldc,slds,sld;
        
        Pi = 4.0*atan(1.0);
        
        lolim = 0.0;
        uplim = 1.0;
        
        summ = 0.0;                     //initialize intergral
        
        scale = dp[0];          //make local copies in case memory moves
        crmaj = dp[1];
        crmin = dp[2];
        trmaj = dp[3];
        trmin = dp[4];
        sldc = dp[5];
        slds = dp[6];
        sld = dp[7];
        delpc = sldc - slds;            //core - shell
        delps = slds - sld;             //shell  - sovent
        bkg = dp[8];

        for(i=0;i<nord;i++) {
                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                yyy = Gauss76Wt[i] * gfn2(zi,crmaj,crmin,trmaj,trmin,delpc,delps,q);
                summ += yyy;
        }
        
        answer = (uplim-lolim)/2.0*summ;
        // normalize by particle volume
        prolatevol = 4.0*Pi/3.0*trmaj*trmin*trmin;
        answer /= prolatevol;
        
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
        
        return answer;
}


/*      StackedDiscsX  :  calculates the form factor of a stacked "tactoid" of core shell disks
like clay platelets that are not exfoliated

*/
double
StackedDiscs(double dp[], double q)
{
        int i;
        double scale,length,bkg,rcore,thick,rhoc,rhol,rhosolv,N,gsd;            //local variables of coefficient wave
        double va,vb,vcyl,summ,yyy,zi,halfheight,d,answer;
        int nord=76;                    //order of integration
        double Pi;
        
        
        Pi = 4.0*atan(1.0);
        
        va = 0.0;
        vb = Pi/2.0;
        
        summ = 0.0;                     //initialize intergral
        
        scale = dp[0];
        rcore = dp[1];
        length = dp[2];
        thick = dp[3];
        rhoc = dp[4];
        rhol = dp[5];
        rhosolv = dp[6];
        N = dp[7];
        gsd = dp[8];
        bkg = dp[9];
        
        d=2.0*thick+length;
        halfheight = length/2.0;
        
        for(i=0;i<nord;i++) {
                zi = ( Gauss76Z[i]*(vb-va) + vb + va )/2.0;
                yyy = Gauss76Wt[i] * Stackdisc_kern(q, rcore, rhoc,rhol,rhosolv, halfheight,thick,zi,gsd,d,N);
                summ += yyy;
        }
        
        answer = (vb-va)/2.0*summ;
        // length is the total core length 
        vcyl=Pi*rcore*rcore*(2.0*thick+length)*N;
        answer /= vcyl;
        //Convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
        
        return answer;
}


/*      LamellarFFX  :  calculates the form factor of a lamellar structure - no S(q) effects included

*/
double
LamellarFF(double dp[], double q)
{
        double scale,del,sig,contr,bkg;         //local variables of coefficient wave
        double inten, qval,Pq;
        double Pi,sldb,sld;
        
        
        Pi = 4.0*atan(1.0);
        scale = dp[0];
        del = dp[1];
        sig = dp[2]*del;
        sldb = dp[3];
        sld = dp[4];
        contr = sldb - sld;
        bkg = dp[5];
        qval=q;
        
        Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig));
        
        inten = 2.0*Pi*scale*Pq/(qval*qval);            //this is now dimensionless...
        
        inten /= del;                   //normalize by the thickness (in A)
        
        inten *= 1.0e8;         // 1/A to 1/cm
        
        return(inten+bkg);
}

/*      LamellarPSX  :  calculates the form factor of a lamellar structure - with S(q) effects included
--- now the proper resolution effects are used - the "default" resolution is turned off (= 0) and the
model is smeared just like any other function
*/
double
LamellarPS(double dp[], double q)
{
        double scale,dd,del,sig,contr,NN,Cp,bkg;                //local variables of coefficient wave
        double inten,qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ;
        double Pi,Euler,dQDefault,fii,sldb,sld;
        int ii,NNint;
//      char buf[256];

        
        Euler = 0.5772156649;           // Euler's constant
//      dQDefault = 0.0025;             //[=] 1/A, q-resolution, default value
        dQDefault = 0.0;
        dQ = dQDefault;
        
        Pi = 4.0*atan(1.0);
        qval = q;
        
        scale = dp[0];
        dd = dp[1];
        del = dp[2];
        sig = dp[3]*del;
        sldb = dp[4];
        sld = dp[5];
        contr = sldb - sld;
        NN = trunc(dp[6]);              //be sure that NN is an integer
        Cp = dp[7];
        bkg = dp[8];
        
        Pq = 2.0*contr*contr/qval/qval*(1.0-cos(qval*del)*exp(-0.5*qval*qval*sig*sig));
        
        NNint = (int)NN;                //cast to an integer for the loop
        
//      sprintf(buf, "qval = %g\r", qval);
//      XOPNotice(buf);
        
        ii=0;
        Sq = 0.0;
        for(ii=1;ii<=(NNint-1);ii+=1) {
        
                fii = (double)ii;               //do I really need to do this?
                
                temp = 0.0;
                alpha = Cp/4.0/Pi/Pi*(log(Pi*fii) + Euler);
                t1 = 2.0*dQ*dQ*dd*dd*alpha;
                t2 = 2.0*qval*qval*dd*dd*alpha;
                t3 = dQ*dQ*dd*dd*fii*fii;
                
                temp = 1.0-fii/NN;
                temp *= cos(dd*qval*fii/(1.0+t1));
                temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
                temp /= sqrt(1.0+t1);
                
                Sq += temp;
        }
        
        Sq *= 2.0;
        Sq += 1.0;
        
        inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval);
        
        inten *= 1.0e8;         // 1/A to 1/cm
        
    return(inten+bkg);
}


/*      LamellarPS_HGX  :  calculates the form factor of a lamellar structure - with S(q) effects included
--- now the proper resolution effects are used - the "default" resolution is turned off (= 0) and the
model is smeared just like any other function
*/
double
LamellarPS_HG(double dp[], double q)
{
        double scale,dd,delT,delH,SLD_T,SLD_H,SLD_S,NN,Cp,bkg;          //local variables of coefficient wave
        double inten,qval,Pq,Sq,alpha,temp,t1,t2,t3,dQ,drh,drt;
        double Pi,Euler,dQDefault,fii;
        int ii,NNint;
        
        
        Euler = 0.5772156649;           // Euler's constant
//      dQDefault = 0.0025;             //[=] 1/A, q-resolution, default value
        dQDefault = 0.0;
        dQ = dQDefault;
        
        Pi = 4.0*atan(1.0);
        qval= q;
        
        scale = dp[0];
        dd = dp[1];
        delT = dp[2];
        delH = dp[3];
        SLD_T = dp[4];
        SLD_H = dp[5];
        SLD_S = dp[6];
        NN = trunc(dp[7]);              //be sure that NN is an integer
        Cp = dp[8];
        bkg = dp[9];
        
        
        drh = SLD_H - SLD_S;
        drt = SLD_T - SLD_S;    //correction 13FEB06 by L.Porcar
        
        Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT);
        Pq *= Pq;
        Pq *= 4.0/(qval*qval);
        
        NNint = (int)NN;                //cast to an integer for the loop
        ii=0;
        Sq = 0.0;
        for(ii=1;ii<=(NNint-1);ii+=1) {
        
                fii = (double)ii;               //do I really need to do this?
                
                temp = 0.0;
                alpha = Cp/4.0/Pi/Pi*(log(Pi*ii) + Euler);
                t1 = 2.0*dQ*dQ*dd*dd*alpha;
                t2 = 2.0*qval*qval*dd*dd*alpha;
                t3 = dQ*dQ*dd*dd*ii*ii;
                
                temp = 1.0-ii/NN;
                temp *= cos(dd*qval*ii/(1.0+t1));
                temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
                temp /= sqrt(1.0+t1);
                
                Sq += temp;
        }
        
        Sq *= 2.0;
        Sq += 1.0;
        
        inten = 2.0*Pi*scale*Pq*Sq/(dd*qval*qval);
        
        inten *= 1.0e8;         // 1/A to 1/cm
        
        return(inten+bkg);
}

/*      LamellarFF_HGX  :  calculates the form factor of a lamellar structure - no S(q) effects included
but extra SLD for head groups is included

*/
double
LamellarFF_HG(double dp[], double q)
{
        double scale,delT,delH,slds,sldh,sldt,bkg;              //local variables of coefficient wave
        double inten, qval,Pq,drh,drt;
        double Pi;
        
        
        Pi = 4.0*atan(1.0);
        qval= q;
        scale = dp[0];
        delT = dp[1];
        delH = dp[2];
        sldt = dp[3];
        sldh = dp[4];
        slds = dp[5];
        bkg = dp[6];
        
        
        drh = sldh - slds;
        drt = sldt - slds;              //correction 13FEB06 by L.Porcar
        
        Pq = drh*(sin(qval*(delH+delT))-sin(qval*delT)) + drt*sin(qval*delT);
        Pq *= Pq;
        Pq *= 4.0/(qval*qval);
        
        inten = 2.0*Pi*scale*Pq/(qval*qval);            //dimensionless...
        
        inten /= 2.0*(delT+delH);                       //normalize by the bilayer thickness
        
        inten *= 1.0e8;         // 1/A to 1/cm
        
        return(inten+bkg);
}

/*      FlexExclVolCylX  :  calculates the form factor of a flexible cylinder with a circular cross section
-- incorporates Wei-Ren Chen's fixes - 2006

        */
double
FlexExclVolCyl(double dp[], double q)
{
        double scale,L,B,bkg,rad,qr,cont,sldc,slds;
        double Pi,flex,crossSect,answer;
        
        
        Pi = 4.0*atan(1.0);
        
        scale = dp[0];          //make local copies in case memory moves
        L = dp[1];
        B = dp[2];
        rad = dp[3];
        sldc = dp[4];
        slds = dp[5];
        cont = sldc-slds;
        bkg = dp[6];
        
    
        qr = q*rad;
        
        flex = Sk_WR(q,L,B);
        
        crossSect = (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
        flex *= crossSect;
        flex *= Pi*rad*rad*L;
        flex *= cont*cont;
        flex *= 1.0e8;
        answer = scale*flex + bkg;
        
        return answer;
}

/*      FlexCyl_EllipX  :  calculates the form factor of a flexible cylinder with an elliptical cross section
-- incorporates Wei-Ren Chen's fixes - 2006

        */
double
FlexCyl_Ellip(double dp[], double q)
{
        double scale,L,B,bkg,rad,qr,cont,ellRatio,slds,sldc;
        double Pi,flex,crossSect,answer;
        
        
        Pi = 4.0*atan(1.0);
        scale = dp[0];          //make local copies in case memory moves
        L = dp[1];
        B = dp[2];
        rad = dp[3];
        ellRatio = dp[4];
        sldc = dp[5];
        slds = dp[6];
        bkg = dp[7];
        
        cont = sldc - slds;
        qr = q*rad;
        
        flex = Sk_WR(q,L,B);
        
        crossSect = EllipticalCross_fn(q,rad,(rad*ellRatio));
        flex *= crossSect;
        flex *= Pi*rad*rad*ellRatio*L;
        flex *= cont*cont;
        flex *= 1.0e8;
        answer = scale*flex + bkg;
        
        return answer;
}

double
EllipticalCross_fn(double qq, double a, double b)
{
    double uplim,lolim,Pi,summ,arg,zi,yyy,answer;
    int i,nord=76;
    
    Pi = 4.0*atan(1.0);
    lolim=0.0;
    uplim=Pi/2.0;
    summ=0.0;
    
    for(i=0;i<nord;i++) {
                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                arg = qq*sqrt(a*a*sin(zi)*sin(zi)+b*b*cos(zi)*cos(zi));
                yyy = pow((2.0 * NR_BessJ1(arg) / arg),2);
                yyy *= Gauss76Wt[i];
                summ += yyy;
    }
    answer = (uplim-lolim)/2.0*summ;
    answer *= 2.0/Pi;
    return(answer);
        
}
/*      FlexCyl_PolyLenX  :  calculates the form factor of a flecible cylinder at the given x-value p->x
the cylinder has a polydisperse Length

*/
double
FlexCyl_PolyLen(double dp[], double q)
{
        int i;
        double scale,radius,length,pd,bkg,lb,delrho,sldc,slds;          //local variables of coefficient wave
        int nord=20;                    //order of integration
        double uplim,lolim;             //upper and lower integration limits
        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
        double range,zz,Pi;
        
        Pi = 4.0*atan(1.0);
        range = 3.4;
        
        summ = 0.0;                     //initialize intergral
        scale = dp[0];                  //make local copies in case memory moves
        length = dp[1];                 //radius
        pd = dp[2];                     // average length
        lb = dp[3];
        radius = dp[4];
        sldc = dp[5];
        slds = dp[6];
        bkg = dp[7];
        
        delrho = sldc - slds;
        zz = (1.0/pd)*(1.0/pd) - 1.0;
        
        lolim = length*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
        if(lolim<0.0) {
                lolim = 0.0;
        }
        if(pd>0.3) {
                range = 3.4 + (pd-0.3)*18.0;
        }
        uplim = length*(1.0+range*pd);
        
        for(i=0;i<nord;i++) {
                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                yyy = Gauss20Wt[i] * FlePolyLen_kernel(q,radius,length,lb,zz,delrho,zi);
                summ += yyy;
        }
        
        answer = (uplim-lolim)/2.0*summ;
        //normalize by average cylinder volume (first moment), using the average length
        Vpoly=Pi*radius*radius*length;
        answer /= Vpoly;
        
        answer *=delrho*delrho;
        
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
        
        return answer;
}

/*      FlexCyl_PolyLenX  :  calculates the form factor of a flexible cylinder at the given x-value p->x
the cylinder has a polydisperse cross sectional radius

*/
double
FlexCyl_PolyRad(double dp[], double q)
{
        int i;
        double scale,radius,length,pd,delrho,bkg,lb,sldc,slds;          //local variables of coefficient wave
        int nord=76;                    //order of integration
        double uplim,lolim;             //upper and lower integration limits
        double summ,zi,yyy,answer,Vpoly;                        //running tally of integration
        double range,zz,Pi;
        
        
        Pi = 4.0*atan(1.0);
        range = 3.4;
        
        summ = 0.0;                     //initialize intergral
        
        scale = dp[0];                  //make local copies in case memory moves
        length = dp[1];                 //radius
        lb = dp[2];                     // average length
        radius = dp[3];
        pd = dp[4];
        sldc = dp[5];
        slds = dp[6];
        bkg = dp[7];
        
        delrho = sldc-slds;
        zz = (1.0/pd)*(1.0/pd) - 1.0;
        
        lolim = radius*(1.0-range*pd);          //set the upper/lower limits to cover the distribution
        if(lolim<0.0) {
                lolim = 0.0;
        }
        if(pd>0.3) {
                range = 3.4 + (pd-0.3)*18.0;
        }
        uplim = radius*(1.0+range*pd);
        
        for(i=0;i<nord;i++) {
                //zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                //yyy = Gauss20Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi);
                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                yyy = Gauss76Wt[i] * FlePolyRad_kernel(q,radius,length,lb,zz,delrho,zi);
                summ += yyy;
        }
        
        answer = (uplim-lolim)/2.0*summ;
        //normalize by average cylinder volume (second moment), using the average radius
        Vpoly = Pi*radius*radius*length*(zz+2.0)/(zz+1.0);
        answer /= Vpoly;
        
        answer *=delrho*delrho;
        
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
        
        return answer;
}



///////////////

//
//     FUNCTION gfn2:    CONTAINS F(Q,A,B,mu)**2  AS GIVEN
//                       BY (53) AND (56,57) IN CHEN AND 
//                       KOTLARCHYK REFERENCE
//
//     <PROLATE ELLIPSOIDS>
//
double
gfn2(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq)
{
        // local variables
        double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,gfn2,pi43,gfn,Pi;

        Pi = 4.0*atan(1.0);
        
        pi43=4.0/3.0*Pi;
        aa = crmaj;
        bb = crmin;
        u2 = (aa*aa*xx*xx + bb*bb*(1.0-xx*xx));
        ut2 = (trmaj*trmaj*xx*xx + trmin*trmin*(1.0-xx*xx));
        uq = sqrt(u2)*qq;
        ut= sqrt(ut2)*qq;
        vc = pi43*aa*bb*bb;
        vt = pi43*trmaj*trmin*trmin;
        if (uq == 0.0){
                siq = 1.0/3.0;
        }else{
                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
        }
        if (ut == 0.0){
                sit = 1.0/3.0;
        }else{
                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut;
        }
        gfnc = 3.0*siq*vc*delpc;
        gfnt = 3.0*sit*vt*delps;
        gfn = gfnc+gfnt;
        gfn2 = gfn*gfn;
        
        return (gfn2);
}

//
//     FUNCTION gfn4:    CONTAINS F(Q,A,B,MU)**2  AS GIVEN
//                       BY (53) & (58-59) IN CHEN AND
//                       KOTLARCHYK REFERENCE
//
//       <OBLATE ELLIPSOID>
// function gfn4 for oblate ellipsoids 
double
gfn4(double xx, double crmaj, double crmin, double trmaj, double trmin, double delpc, double delps, double qq)
{
        // local variables
        double aa,bb,u2,ut2,uq,ut,vc,vt,siq,sit,gfnc,gfnt,tgfn,gfn4,pi43,Pi;

        Pi = 4.0*atan(1.0);
        pi43=4.0/3.0*Pi;
        aa = crmaj;
        bb = crmin;
        u2 = (bb*bb*xx*xx + aa*aa*(1.0-xx*xx));
        ut2 = (trmin*trmin*xx*xx + trmaj*trmaj*(1.0-xx*xx));
        uq = sqrt(u2)*qq;
        ut= sqrt(ut2)*qq;
        vc = pi43*aa*aa*bb;
        vt = pi43*trmaj*trmaj*trmin;
        if (uq == 0.0){
                siq = 1.0/3.0;
        }else{
                siq = (sin(uq)/uq/uq - cos(uq)/uq)/uq;
        }
        if (ut == 0.0){
                sit = 1.0/3.0;
        }else{
                sit = (sin(ut)/ut/ut - cos(ut)/ut)/ut;
        }
        gfnc = 3.0*siq*vc*delpc;
        gfnt = 3.0*sit*vt*delps;
        tgfn = gfnc+gfnt;
        gfn4 = tgfn*tgfn;
        
        return (gfn4);
}

double
FlePolyLen_kernel(double q, double radius, double length, double lb, double zz, double delrho, double zi)
{
    double Pq,vcyl,dl;
    double Pi,qr;
    
    Pi = 4.0*atan(1.0);
    qr = q*radius;
    
    Pq = Sk_WR(q,zi,lb);                //does not have cross section term
    if (qr !=0){
        Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
    } 
    vcyl=Pi*radius*radius*zi;
    Pq *= vcyl*vcyl;
    
    dl = SchulzPoint_cpr(zi,length,zz);
    return (Pq*dl);     
        
}

double
FlePolyRad_kernel(double q, double ravg, double Lc, double Lb, double zz, double delrho, double zi)
{
    double Pq,vcyl,dr;
    double Pi,qr;
    
    Pi = 4.0*atan(1.0);
    qr = q*zi;
    
    Pq = Sk_WR(q,Lc,Lb);                //does not have cross section term
    if (qr !=0){
        Pq *= (2.0*NR_BessJ1(qr)/qr)*(2.0*NR_BessJ1(qr)/qr);
    }

    vcyl=Pi*zi*zi*Lc;
    Pq *= vcyl*vcyl;
    
    dr = SchulzPoint_cpr(zi,ravg,zz);
    return (Pq*dr);     
        
}

double
CSCylIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length)
{
        double answer,halfheight,Pi;
        double lolim,uplim,summ,yyy,zi;
        int nord,i;
        
        // set up the integration end points 
        Pi = 4.0*atan(1.0);
        nord = 76;
        lolim = 0.0;
        uplim = Pi/2.0;
        halfheight = length/2.0;
        
        summ = 0.0;                             // initialize integral
        i=0;
        for(i=0;i<nord;i++) {
                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                yyy = Gauss76Wt[i] * CScyl(qq, rad, radthick, facthick, rhoc,rhos,rhosolv, halfheight, zi);
                summ += yyy;
        }
        
        // calculate value of integral to return
        answer = (uplim-lolim)/2.0*summ;
        return (answer);
}

double
CScyl(double qq, double rad, double radthick, double facthick, double rhoc, double rhos, double rhosolv, double length, double dum)
{       
        // qq is the q-value for the calculation (1/A)
        // radius is the core radius of the cylinder (A)
        //  radthick and facthick are the radial and face layer thicknesses
        // rho(n) are the respective SLD's
        // length is the *Half* CORE-LENGTH of the cylinder 
        // dum is the dummy variable for the integration (theta)

        double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval;
        double Pi;
        
        Pi = 4.0*atan(1.0); 
        
        dr1 = rhoc-rhos;
        dr2 = rhos-rhosolv;
        vol1 = Pi*rad*rad*(2.0*length);
        vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick);
        
        besarg1 = qq*rad*sin(dum);
        besarg2 = qq*(rad+radthick)*sin(dum);
        sinarg1 = qq*length*cos(dum);
        sinarg2 = qq*(length+facthick)*cos(dum);
        if (besarg1 == 0.0){
                be1 = 0.5;
        }else{
                be1 = NR_BessJ1(besarg1)/besarg1;
        }
        if (besarg2 == 0.0){
                be2 = 0.5;
        }else{
                be2 = NR_BessJ1(besarg2)/besarg2;
        }
        if (sinarg1 == 0.0){
                si1 = 1.0;
        }else{
                si1 = sin(sinarg1)/sinarg1;
        }
        if (sinarg2 == 0.0){
                si2 = 1.0;
        }else{
                si2 = sin(sinarg2)/sinarg2;
        }

        t1 = 2.0*vol1*dr1*si1*be1;
        t2 = 2.0*vol2*dr2*si2*be2;

        retval = ((t1+t2)*(t1+t2))*sin(dum);
        return (retval);
    
}


double
CoreShellCylKernel(double qq, double rcore, double thick, double rhoc, double rhos, double rhosolv, double length, double dum)
{

    double dr1,dr2,besarg1,besarg2,vol1,vol2,sinarg1,sinarg2,si1,si2,be1,be2,t1,t2,retval;
    double Pi;
    
    Pi = 4.0*atan(1.0);
    
    dr1 = rhoc-rhos;
    dr2 = rhos-rhosolv;
    vol1 = Pi*rcore*rcore*(2.0*length);
    vol2 = Pi*(rcore+thick)*(rcore+thick)*(2.0*length+2.0*thick);
    
    besarg1 = qq*rcore*sin(dum);
    besarg2 = qq*(rcore+thick)*sin(dum);
    sinarg1 = qq*length*cos(dum);
    sinarg2 = qq*(length+thick)*cos(dum);

        if (besarg1 == 0.0){
                be1 = 0.5;
        }else{
                be1 = NR_BessJ1(besarg1)/besarg1;
        }
        if (besarg2 == 0.0){
                be2 = 0.5;
        }else{
                be2 = NR_BessJ1(besarg2)/besarg2;
        }
        if (sinarg1 == 0.0){
                si1 = 1.0;
        }else{
                si1 = sin(sinarg1)/sinarg1;
        }
        if (sinarg2 == 0.0){
                si2 = 1.0;
        }else{
                si2 = sin(sinarg2)/sinarg2;
        }

    t1 = 2.0*vol1*dr1*si1*be1;
    t2 = 2.0*vol2*dr2*si2*be2;

    retval = ((t1+t2)*(t1+t2))*sin(dum);
        
    return (retval);
}

double
Cyl_PolyLenKernel(double q, double radius, double len_avg, double zz, double delrho, double dumLen)
{
        
    double halfheight,uplim,lolim,zi,summ,yyy,Pi;
    double answer,dr,Vcyl;
    int i,nord;
    
    Pi = 4.0*atan(1.0);
    lolim = 0.0;
    uplim = Pi/2.0;
    halfheight = dumLen/2.0;
    nord=20;
    summ = 0.0;
    
    //do the cylinder orientational average
    for(i=0;i<nord;i++) {
                zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                yyy = Gauss20Wt[i] * CylKernel(q, radius, halfheight, zi);
                summ += yyy;
    }
    answer = (uplim-lolim)/2.0*summ;
    // Multiply by contrast^2
    answer *= delrho*delrho;
    // don't do the normal scaling to volume here
    // instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder
    Vcyl = Pi*radius*radius*dumLen;
    answer *= Vcyl*Vcyl;
    
    dr = SchulzPoint_cpr(dumLen,len_avg,zz);
    return(dr*answer);
}


double
Stackdisc_kern(double qq, double rcore, double rhoc, double rhol, double rhosolv, double length, double thick, double dum, double gsd, double d, double N)
{               
        // qq is the q-value for the calculation (1/A)
        // rcore is the core radius of the cylinder (A)
        // rho(n) are the respective SLD's
        // length is the *Half* CORE-LENGTH of the cylinder = L (A)
        // dum is the dummy variable for the integration (x in Feigin's notation)

        //Local variables
        double totald,dr1,dr2,besarg1,besarg2,be1,be2,si1,si2,area,sinarg1,sinarg2,t1,t2,retval,sqq,dexpt;
        double Pi;
        int kk;
        
        Pi = 4.0*atan(1.0);
        
        dr1 = rhoc-rhosolv;
        dr2 = rhol-rhosolv;
        area = Pi*rcore*rcore;
        totald=2.0*(thick+length);
        
        besarg1 = qq*rcore*sin(dum);
        besarg2 = qq*rcore*sin(dum);
        
        sinarg1 = qq*length*cos(dum);
        sinarg2 = qq*(length+thick)*cos(dum);

        if (besarg1 == 0.0){
                be1 = 0.5;
        }else{
                be1 = NR_BessJ1(besarg1)/besarg1;
        }
        if (besarg2 == 0.0){
                be2 = 0.5;
        }else{
                be2 = NR_BessJ1(besarg2)/besarg2;
        }
        if (sinarg1 == 0.0){
                si1 = 1.0;
        }else{
                si1 = sin(sinarg1)/sinarg1;
        }
        if (sinarg2 == 0.0){
                si2 = 1.0;
        }else{
                si2 = sin(sinarg2)/sinarg2;
        }

        t1 = 2.0*area*(2.0*length)*dr1*(si1)*(be1);
        t2 = 2.0*area*dr2*(totald*si2-2.0*length*si1)*(be2);

        retval =((t1+t2)*(t1+t2))*sin(dum);
        
        // loop for the structure facture S(q)
        sqq=0.0;
        for(kk=1;kk<N;kk+=1) {
                dexpt=qq*cos(dum)*qq*cos(dum)*d*d*gsd*gsd*kk/2.0;
                sqq=sqq+(N-kk)*cos(qq*cos(dum)*d*kk)*exp(-1.*dexpt);
        }                       
        
        // end of loop for S(q)
        sqq=1.0+2.0*sqq/N;
        retval *= sqq;
    
        return(retval);
}


double
Cyl_PolyRadKernel(double q, double radius, double length, double zz, double delrho, double dumRad)
{
        
    double halfheight,uplim,lolim,zi,summ,yyy,Pi;
    double answer,dr,Vcyl;
    int i,nord;
    
    Pi = 4.0*atan(1.0);
    lolim = 0.0;
    uplim = Pi/2.0;
    halfheight = length/2.0;
        //    nord=20;
    nord=76;
    summ = 0.0;
    
    //do the cylinder orientational average
        //    for(i=0;i<nord;i++) {
        //            zi = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
        //            yyy = Gauss20Wt[i] * CylKernel(q, dumRad, halfheight, zi);
        //            summ += yyy;
        //    }
    for(i=0;i<nord;i++) {
                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                yyy = Gauss76Wt[i] * CylKernel(q, dumRad, halfheight, zi);
                summ += yyy;
    }
    answer = (uplim-lolim)/2.0*summ;
    // Multiply by contrast^2
    answer *= delrho*delrho;
    // don't do the normal scaling to volume here
    // instead, multiply by VCyl^2 to get rid of the normalization for this radius of cylinder
    Vcyl = Pi*dumRad*dumRad*length;
    answer *= Vcyl*Vcyl;
    
    dr = SchulzPoint_cpr(dumRad,radius,zz);
    return(dr*answer);
}

double
SchulzPoint_cpr(double dumRad, double radius, double zz)
{
    double dr;
    
    dr = zz*log(dumRad) - gammaln(zz+1.0) + (zz+1.0)*log((zz+1.0)/radius)-(dumRad/radius*(zz+1.0));
    return(exp(dr));
}


double
EllipsoidKernel(double qq, double a, double nua, double dum)
{
    double arg,nu,retval;               //local variables
        
    nu = nua/a;
    arg = qq*a*sqrt(1.0+dum*dum*(nu*nu-1.0));
    if (arg == 0.0){
        retval =1.0/3.0;
    }else{
        retval = (sin(arg)-arg*cos(arg))/(arg*arg*arg);
    }
    retval *= retval;
    retval *= 9.0;

    return(retval);
}//Function EllipsoidKernel()

double
HolCylKernel(double qq, double rcore, double rshell, double length, double dum)
{
    double gamma,arg1,arg2,lam1,lam2,psi,sinarg,t2,retval;              //local variables
    
    gamma = rcore/rshell;
    arg1 = qq*rshell*sqrt(1.0-dum*dum);         //1=shell (outer radius)
    arg2 = qq*rcore*sqrt(1.0-dum*dum);                  //2=core (inner radius)
    if (arg1 == 0.0){
        lam1 = 1.0;
    }else{
        lam1 = 2.0*NR_BessJ1(arg1)/arg1;
    }
    if (arg2 == 0.0){
        lam2 = 1.0;
    }else{
        lam2 = 2.0*NR_BessJ1(arg2)/arg2;
    }
    //Todo: Need to check psi behavior as gamma goes to 1.
    psi = 1.0/(1.0-gamma*gamma)*(lam1 -  gamma*gamma*lam2);             //SRK 10/19/00
    sinarg = qq*length*dum/2.0;
    if (sinarg == 0.0){
        t2 = 1.0;
    }else{
        t2 = sin(sinarg)/sinarg;
    }

    retval = psi*psi*t2*t2;
    
    return(retval);
}//Function HolCylKernel()

double
PPKernel(double aa, double mu, double uu)
{
    // mu passed in is really mu*sqrt(1-sig^2)
    double arg1,arg2,Pi,tmp1,tmp2;                      //local variables
        
    Pi = 4.0*atan(1.0);
    
    //handle arg=0 separately, as sin(t)/t -> 1 as t->0
    arg1 = (mu/2.0)*cos(Pi*uu/2.0);
    arg2 = (mu*aa/2.0)*sin(Pi*uu/2.0);
    if(arg1==0.0) {
                tmp1 = 1.0;
        } else {
                tmp1 = sin(arg1)*sin(arg1)/arg1/arg1;
    }

    if (arg2==0.0) {
                tmp2 = 1.0;
        } else {
                tmp2 = sin(arg2)*sin(arg2)/arg2/arg2;
    }
        
    return (tmp1*tmp2);
        
}//Function PPKernel()


double
TriaxialKernel(double q, double aa, double bb, double cc, double dx, double dy)
{
        
    double arg,val,pi;                  //local variables
        
    pi = 4.0*atan(1.0);
    
    arg = aa*aa*cos(pi*dx/2)*cos(pi*dx/2);
    arg += bb*bb*sin(pi*dx/2)*sin(pi*dx/2)*(1-dy*dy);
    arg += cc*cc*dy*dy;
    arg = q*sqrt(arg);
    if (arg == 0.0){
        val = 1.0; // as arg --> 0, val should go to 1.0
    }else{
        val = 9.0 * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) ) * ( (sin(arg) - arg*cos(arg))/(arg*arg*arg) );
    }
    return (val);
        
}//Function TriaxialKernel()


double
CylKernel(double qq, double rr,double h, double theta)
{
        
        // qq is the q-value for the calculation (1/A)
        // rr is the radius of the cylinder (A)
        // h is the HALF-LENGTH of the cylinder = L/2 (A)

    double besarg,bj,retval,d1,t1,b1,t2,b2,siarg,be,si;          //Local variables


    besarg = qq*rr*sin(theta);
    siarg = qq * h * cos(theta);
    bj =NR_BessJ1(besarg);
        
        //* Computing 2nd power */
    d1 = sin(siarg);
    t1 = d1 * d1;
        //* Computing 2nd power */
    d1 = bj;
    t2 = d1 * d1 * 4.0 * sin(theta);
        //* Computing 2nd power */
    d1 = siarg;
    b1 = d1 * d1;
        //* Computing 2nd power */
    d1 = qq * rr * sin(theta);
    b2 = d1 * d1;
    if (besarg == 0.0){
        be = sin(theta);
    }else{
        be = t2 / b2;
    }
    if (siarg == 0.0){
        si = 1.0;
    }else{
        si = t1 / b1;
    }
    retval = be * si;

    return (retval);
        
}//Function CylKernel()

double
EllipCylKernel(double qq, double ra,double nu, double theta)
{
        //this is the function LAMBDA1^2 in Feigin's notation   
        // qq is the q-value for the calculation (1/A)
        // ra is the transformed radius"a" in Feigin's notation
        // nu is the ratio (major radius)/(minor radius) of the Ellipsoid [=] ---
        // theta is the dummy variable of the integration
        
        double retval,arg;               //Local variables 
        
        arg = qq*ra*sqrt((1.0+nu*nu)/2+(1.0-nu*nu)*cos(theta)/2);
        if (arg == 0.0){
                retval = 1.0;
        }else{
                retval = 2.0*NR_BessJ1(arg)/arg;
        }

        //square it
        retval *= retval;
        
    return(retval);
        
}//Function EllipCylKernel()

double NR_BessJ1(double x)
{
        double ax,z;
        double xx,y,ans,ans1,ans2;
        
        if ((ax=fabs(x)) < 8.0) {
                y=x*x;
                ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
                                                                                                  +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
                ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
                                                                                           +y*(99447.43394+y*(376.9991397+y*1.0))));
                ans=ans1/ans2;
        } else {
                z=8.0/ax;
                y=z*z;
                xx=ax-2.356194491;
                ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
                                                                   +y*(0.2457520174e-5+y*(-0.240337019e-6))));
                ans2=0.04687499995+y*(-0.2002690873e-3
                                                          +y*(0.8449199096e-5+y*(-0.88228987e-6
                                                                                                         +y*0.105787412e-6)));
                ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
                if (x < 0.0) ans = -ans;
        }
        
        return(ans);
}

/*      Lamellar_ParaCrystal - Pedersen's model

*/
double
Lamellar_ParaCrystal(double w[], double q)
{
//       Input (fitting) variables are:
        //[0] scale factor
        //[1]   thickness
        //[2]   number of layers
        //[3]   spacing between layers
        //[4]   polydispersity of spacing
        //[5] SLD lamellar
        //[6] SLD solvent
        //[7] incoherent background
//      give them nice names
        double inten,qval,scale,th,nl,davg,pd,contr,bkg,xn;
        double xi,ww,Pbil,Znq,Snq,an,sldLayer,sldSolvent,pi;
        long n1,n2;
        
        pi = 4.0*atan(1.0);
        scale = w[0];
        th = w[1];
        nl = w[2];
        davg = w[3];
        pd = w[4];
        sldLayer = w[5];
        sldSolvent = w[6];
        bkg = w[7];
        
        contr = w[5] - w[6];
        qval = q;
                
        //get the fractional part of nl, to determine the "mixing" of N's
        
        n1 = trunc(nl);         //rounds towards zero
        n2 = n1 + 1;
        xn = (double)n2 - nl;                   //fractional contribution of n1
        
        ww = exp(-qval*qval*pd*pd*davg*davg/2.0);
        
        //calculate the n1 contribution
        an = paraCryst_an(ww,qval,davg,n1);
        Snq = paraCryst_sn(ww,qval,davg,n1,an);
        
        Znq = xn*Snq;
        
        //calculate the n2 contribution
        an = paraCryst_an(ww,qval,davg,n2);
        Snq = paraCryst_sn(ww,qval,davg,n2,an);
        
        Znq += (1.0-xn)*Snq;
        
        //and the independent contribution
        Znq += (1.0-ww*ww)/(1.0+ww*ww-2.0*ww*cos(qval*davg));
        
        //the limit when NL approaches infinity
//      Zq = (1-ww^2)/(1+ww^2-2*ww*cos(qval*davg))
        
        xi = th/2.0;            //use 1/2 the bilayer thickness
        Pbil = (sin(qval*xi)/(qval*xi))*(sin(qval*xi)/(qval*xi));
        
        inten = 2.0*pi*contr*contr*Pbil*Znq/(qval*qval);
        inten *= 1.0e8;
        
        return(scale*inten+bkg);
}

// functions for the lamellar paracrystal model
double
paraCryst_sn(double ww, double qval, double davg, long nl, double an) {
        
        double Snq;
        
        Snq = an/( (double)nl*pow((1.0+ww*ww-2.0*ww*cos(qval*davg)),2) );
        
        return(Snq);
}


double
paraCryst_an(double ww, double qval, double davg, long nl) {
        
        double an;
        
        an = 4.0*ww*ww - 2.0*(ww*ww*ww+ww)*cos(qval*davg);
        an -= 4.0*pow(ww,(nl+2))*cos((double)nl*qval*davg);
        an += 2.0*pow(ww,(nl+3))*cos((double)(nl-1)*qval*davg);
        an += 2.0*pow(ww,(nl+1))*cos((double)(nl+1)*qval*davg);
        
        return(an);
}


/*      Spherocylinder  :

Uses 76 pt Gaussian quadrature for both integrals
*/
double
Spherocylinder(double w[], double x)
{
        int i,j;
        double Pi;
        double scale,contr,bkg,sldc,slds;
        double len,rad,hDist,endRad;
        int nordi=76;                   //order of integration
        int nordj=76;
        double va,vb;           //upper and lower integration limits
        double summ,zi,yyy,answer;                      //running tally of integration
        double summj,vaj,vbj,zij;                       //for the inner integration
        double SphCyl_tmp[7],arg1,arg2,inner,be;
        
        
        scale = w[0];
        rad = w[1];
        len = w[2];
        sldc = w[3];
        slds = w[4];
        bkg = w[5];
        
        SphCyl_tmp[0] = w[0];
        SphCyl_tmp[1] = w[1];
        SphCyl_tmp[2] = w[2];
        SphCyl_tmp[3] = w[1];           //end radius is same as cylinder radius
        SphCyl_tmp[4] = w[3];
        SphCyl_tmp[5] = w[4];
        SphCyl_tmp[6] = w[5];
        
        hDist = 0;              //by definition for this model
        endRad = rad;
        
        contr = sldc-slds;
        
        Pi = 4.0*atan(1.0);
        va = 0.0;
        vb = Pi/2.0;            //orintational average, outer integral
        vaj = -1.0*hDist/endRad;
        vbj = 1.0;              //endpoints of inner integral

        summ = 0.0;                     //initialize intergral

        for(i=0;i<nordi;i++) {
                //setup inner integral over the ellipsoidal cross-section
                summj=0.0;
                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
                
                for(j=0;j<nordj;j++) {
                        //20 gauss points for the inner integral
                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
                        yyy = Gauss76Wt[j] * SphCyl_kernel(SphCyl_tmp,x,zij,zi);
                        summj += yyy;
                }
                //now calculate the value of the inner integral
                inner = (vbj-vaj)/2.0*summj;
                inner *= 4.0*Pi*endRad*endRad*endRad;
                
                //now calculate outer integrand
                arg1 = x*len/2.0*cos(zi);
                arg2 = x*rad*sin(zi);
                yyy = inner;

                if(arg2 == 0) {
                        be = 0.5;
                } else {
                        be = NR_BessJ1(arg2)/arg2;
                }
                
                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
                        yyy += Pi*rad*rad*len*2.0*be;
                } else {
                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
                }
                yyy *= yyy;
                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
                yyy *= Gauss76Wt[i];
                summ += yyy;
        }               //final scaling is done at the end of the function, after the NT_FP64 case
        
        answer = (vb-va)/2.0*summ;

        answer /= Pi*rad*rad*len + Pi*4.0*endRad*endRad*endRad/3.0;             //divide by volume
        answer *= 1.0e8;                //convert to cm^-1
        answer *= contr*contr;
        answer *= scale;
        answer += bkg;
                
        return answer;
}


// inner integral of the sphereocylinder model, special case of hDist = 0
//
double
SphCyl_kernel(double w[], double x, double tt, double theta) {

        double val,arg1,arg2;
        double scale,bkg,sldc,slds;
        double len,rad,hDist,endRad,be;
        scale = w[0];
        rad = w[1];
        len = w[2];
        endRad = w[3];
        sldc = w[4];
        slds = w[5];
        bkg = w[6];
        
        hDist = 0.0;
                
        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0);
        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt);
        
        if(arg2 == 0) {
                be = 0.5;
        } else {
                be = NR_BessJ1(arg2)/arg2;
        }
        val = cos(arg1)*(1.0-tt*tt)*be;
        
        return(val);
}


/*      Convex Lens  : special case where L ~ 0 and hDist < 0

Uses 76 pt Gaussian quadrature for both integrals
*/
double
ConvexLens(double w[], double x)
{
        int i,j;
        double Pi;
        double scale,contr,bkg,sldc,slds;
        double len,rad,hDist,endRad;
        int nordi=76;                   //order of integration
        int nordj=76;
        double va,vb;           //upper and lower integration limits
        double summ,zi,yyy,answer;                      //running tally of integration
        double summj,vaj,vbj,zij;                       //for the inner integration
        double CLens_tmp[7],arg1,arg2,inner,hh,be;
        
        
        scale = w[0];
        rad = w[1];
//      len = w[2]
        endRad = w[2];
        sldc = w[3];
        slds = w[4];
        bkg = w[5];
        
        len = 0.01;
        
        CLens_tmp[0] = w[0];
        CLens_tmp[1] = w[1];
        CLens_tmp[2] = len;                     //length is some small number, essentially zero
        CLens_tmp[3] = w[2];
        CLens_tmp[4] = w[3];
        CLens_tmp[5] = w[4];
        CLens_tmp[6] = w[5];
                
        hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad));         //by definition for this model
        
        contr = sldc-slds;
        
        Pi = 4.0*atan(1.0);
        va = 0.0;
        vb = Pi/2.0;            //orintational average, outer integral
        vaj = -1.0*hDist/endRad;
        vbj = 1.0;              //endpoints of inner integral

        summ = 0.0;                     //initialize intergral

        for(i=0;i<nordi;i++) {
                //setup inner integral over the ellipsoidal cross-section
                summj=0.0;
                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
                
                for(j=0;j<nordj;j++) {
                        //20 gauss points for the inner integral
                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
                        yyy = Gauss76Wt[j] * ConvLens_kernel(CLens_tmp,x,zij,zi);
                        summj += yyy;
                }
                //now calculate the value of the inner integral
                inner = (vbj-vaj)/2.0*summj;
                inner *= 4.0*Pi*endRad*endRad*endRad;
                
                //now calculate outer integrand
                arg1 = x*len/2.0*cos(zi);
                arg2 = x*rad*sin(zi);
                yyy = inner;
                
                if(arg2 == 0) {
                        be = 0.5;
                } else {
                        be = NR_BessJ1(arg2)/arg2;
                }
                
                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
                        yyy += Pi*rad*rad*len*2.0*be;
                } else {
                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
                }
                yyy *= yyy;
                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
                yyy *= Gauss76Wt[i];
                summ += yyy;
        }               //final scaling is done at the end of the function, after the NT_FP64 case
        
        answer = (vb-va)/2.0*summ;

        hh = fabs(hDist);               //need positive value for spherical cap volume
        answer /= 2.0*(1.0/3.0*Pi*(endRad-hh)*(endRad-hh)*(2.0*endRad+hh));             //divide by volume
        answer *= 1.0e8;                //convert to cm^-1
        answer *= contr*contr;
        answer *= scale;
        answer += bkg;
                
        return answer;
}

/*      Capped Cylinder  : special case where L is nonzero and hDist < 0

 -- uses the same Kernel as the Convex Lens
 
Uses 76 pt Gaussian quadrature for both integrals
*/
double
CappedCylinder(double w[], double x)
{
        int i,j;
        double Pi;
        double scale,contr,bkg,sldc,slds;
        double len,rad,hDist,endRad;
        int nordi=76;                   //order of integration
        int nordj=76;
        double va,vb;           //upper and lower integration limits
        double summ,zi,yyy,answer;                      //running tally of integration
        double summj,vaj,vbj,zij;                       //for the inner integration
        double arg1,arg2,inner,hh,be;
        
        
        scale = w[0];
        rad = w[1];
        len = w[2];
        endRad = w[3];
        sldc = w[4];
        slds = w[5];
        bkg = w[6];
                
        hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad));         //by definition for this model
        
        contr = sldc-slds;
        
        Pi = 4.0*atan(1.0);
        va = 0.0;
        vb = Pi/2.0;            //orintational average, outer integral
        vaj = -1.0*hDist/endRad;
        vbj = 1.0;              //endpoints of inner integral

        summ = 0.0;                     //initialize intergral

        for(i=0;i<nordi;i++) {
                //setup inner integral over the ellipsoidal cross-section
                summj=0.0;
                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
                
                for(j=0;j<nordj;j++) {
                        //20 gauss points for the inner integral
                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
                        yyy = Gauss76Wt[j] * ConvLens_kernel(w,x,zij,zi);               //uses the same kernel as ConvexLens, except here L != 0
                        summj += yyy;
                }
                //now calculate the value of the inner integral
                inner = (vbj-vaj)/2.0*summj;
                inner *= 4.0*Pi*endRad*endRad*endRad;
                
                //now calculate outer integrand
                arg1 = x*len/2.0*cos(zi);
                arg2 = x*rad*sin(zi);
                yyy = inner;
                
                if(arg2 == 0) {
                        be = 0.5;
                } else {
                        be = NR_BessJ1(arg2)/arg2;
                }
                
                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
                        yyy += Pi*rad*rad*len*2.0*be;
                } else {
                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
                }
                
                
                
                yyy *= yyy;
                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
                yyy *= Gauss76Wt[i];
                summ += yyy;
        }               //final scaling is done at the end of the function, after the NT_FP64 case
        
        answer = (vb-va)/2.0*summ;

        hh = fabs(hDist);               //need positive value for spherical cap volume
        answer /= Pi*rad*rad*len + 2.0*(1.0/3.0*Pi*(endRad-hh)*(endRad-hh)*(2.0*endRad+hh));            //divide by volume
        answer *= 1.0e8;                //convert to cm^-1
        answer *= contr*contr;
        answer *= scale;
        answer += bkg;
                
        return answer;
}



// inner integral of the ConvexLens model, special case where L ~ 0 and hDist < 0
//
double
ConvLens_kernel(double w[], double x, double tt, double theta) {

        double val,arg1,arg2;
        double scale,bkg,sldc,slds;
        double len,rad,hDist,endRad,be;
        scale = w[0];
        rad = w[1];
        len = w[2];
        endRad = w[3];
        sldc = w[4];
        slds = w[5];
        bkg = w[6];
        
        hDist = -1.0*sqrt(fabs(endRad*endRad-rad*rad));
                
        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0);
        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt);
        
        if(arg2 == 0) {
                be = 0.5;
        } else {
                be = NR_BessJ1(arg2)/arg2;
        }
        
        val = cos(arg1)*(1.0-tt*tt)*be;
        
        return(val);
}


/*      Dumbbell  : special case where L ~ 0 and hDist > 0

Uses 76 pt Gaussian quadrature for both integrals
*/
double
Dumbbell(double w[], double x)
{
        int i,j;
        double Pi;
        double scale,contr,bkg,sldc,slds;
        double len,rad,hDist,endRad;
        int nordi=76;                   //order of integration
        int nordj=76;
        double va,vb;           //upper and lower integration limits
        double summ,zi,yyy,answer;                      //running tally of integration
        double summj,vaj,vbj,zij;                       //for the inner integration
        double Dumb_tmp[7],arg1,arg2,inner,be;
        
        
        scale = w[0];
        rad = w[1];
//      len = w[2]
        endRad = w[2];
        sldc = w[3];
        slds = w[4];
        bkg = w[5];
        
        len = 0.01;
        
        Dumb_tmp[0] = w[0];
        Dumb_tmp[1] = w[1];
        Dumb_tmp[2] = len;              //length is some small number, essentially zero
        Dumb_tmp[3] = w[2];
        Dumb_tmp[4] = w[3];
        Dumb_tmp[5] = w[4];
        Dumb_tmp[6] = w[5];
                        
        hDist = sqrt(fabs(endRad*endRad-rad*rad));              //by definition for this model
        
        contr = sldc-slds;
        
        Pi = 4.0*atan(1.0);
        va = 0.0;
        vb = Pi/2.0;            //orintational average, outer integral
        vaj = -1.0*hDist/endRad;
        vbj = 1.0;              //endpoints of inner integral

        summ = 0.0;                     //initialize intergral

        for(i=0;i<nordi;i++) {
                //setup inner integral over the ellipsoidal cross-section
                summj=0.0;
                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
                
                for(j=0;j<nordj;j++) {
                        //20 gauss points for the inner integral
                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
                        yyy = Gauss76Wt[j] * Dumb_kernel(Dumb_tmp,x,zij,zi);
                        summj += yyy;
                }
                //now calculate the value of the inner integral
                inner = (vbj-vaj)/2.0*summj;
                inner *= 4.0*Pi*endRad*endRad*endRad;
                
                //now calculate outer integrand
                arg1 = x*len/2.0*cos(zi);
                arg2 = x*rad*sin(zi);
                yyy = inner;
                
                if(arg2 == 0) {
                        be = 0.5;
                } else {
                        be = NR_BessJ1(arg2)/arg2;
                }
                
                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
                        yyy += Pi*rad*rad*len*2.0*be;
                } else {
                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
                }
                yyy *= yyy;
                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
                yyy *= Gauss76Wt[i];
                summ += yyy;
        }               //final scaling is done at the end of the function, after the NT_FP64 case
        
        answer = (vb-va)/2.0*summ;

        answer /= 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*hDist*hDist/3.0);              //divide by volume
        answer *= 1.0e8;                //convert to cm^-1
        answer *= contr*contr;
        answer *= scale;
        answer += bkg;
                
        return answer;
}


/*      Barbell  : "normal" case where L is nonzero 0 and hDist > 0

-- uses the same kernel as the Dumbbell case

Uses 76 pt Gaussian quadrature for both integrals
*/
double
Barbell(double w[], double x)
{
        int i,j;
        double Pi;
        double scale,contr,bkg,sldc,slds;
        double len,rad,hDist,endRad;
        int nordi=76;                   //order of integration
        int nordj=76;
        double va,vb;           //upper and lower integration limits
        double summ,zi,yyy,answer;                      //running tally of integration
        double summj,vaj,vbj,zij;                       //for the inner integration
        double arg1,arg2,inner,be;
        
        
        scale = w[0];
        rad = w[1];
        len = w[2];
        endRad = w[3];
        sldc = w[4];
        slds = w[5];
        bkg = w[6];
                        
        hDist = sqrt(fabs(endRad*endRad-rad*rad));              //by definition for this model
        
        contr = sldc-slds;
        
        Pi = 4.0*atan(1.0);
        va = 0.0;
        vb = Pi/2.0;            //orintational average, outer integral
        vaj = -1.0*hDist/endRad;
        vbj = 1.0;              //endpoints of inner integral

        summ = 0.0;                     //initialize intergral

        for(i=0;i<nordi;i++) {
                //setup inner integral over the ellipsoidal cross-section
                summj=0.0;
                zi = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;             //the "theta" dummy
                
                for(j=0;j<nordj;j++) {
                        //20 gauss points for the inner integral
                        zij = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;                //the "t" dummy
                        yyy = Gauss76Wt[j] * Dumb_kernel(w,x,zij,zi);           //uses the same Kernel as the Dumbbell, here L>0
                        summj += yyy;
                }
                //now calculate the value of the inner integral
                inner = (vbj-vaj)/2.0*summj;
                inner *= 4.0*Pi*endRad*endRad*endRad;
                
                //now calculate outer integrand
                arg1 = x*len/2.0*cos(zi);
                arg2 = x*rad*sin(zi);
                yyy = inner;
                
                if(arg2 == 0) {
                        be = 0.5;
                } else {
                        be = NR_BessJ1(arg2)/arg2;
                }
                
                if(arg1 == 0.0) {               //limiting value of sinc(0) is 1; sinc is not defined in math.h
                        yyy += Pi*rad*rad*len*2.0*be;
                } else {
                        yyy += Pi*rad*rad*len*sin(arg1)/arg1*2.0*be;
                }
                yyy *= yyy;
                yyy *= sin(zi);         // = |A(q)|^2*sin(theta)
                yyy *= Gauss76Wt[i];
                summ += yyy;
        }               //final scaling is done at the end of the function, after the NT_FP64 case
        
        answer = (vb-va)/2.0*summ;

        answer /= Pi*rad*rad*len + 2.0*Pi*(2.0*endRad*endRad*endRad/3.0+endRad*endRad*hDist-hDist*hDist*hDist/3.0);             //divide by volume
        answer *= 1.0e8;                //convert to cm^-1
        answer *= contr*contr;
        answer *= scale;
        answer += bkg;
                
        return answer;
}



// inner integral of the Dumbbell model, special case where L ~ 0 and hDist > 0
//
// inner integral of the Barbell model if L is nonzero
//
double
Dumb_kernel(double w[], double x, double tt, double theta) {

        double val,arg1,arg2;
        double scale,bkg,sldc,slds;
        double len,rad,hDist,endRad,be;
        scale = w[0];
        rad = w[1];
        len = w[2];
        endRad = w[3];
        sldc = w[4];
        slds = w[5];
        bkg = w[6];
        
        hDist = sqrt(fabs(endRad*endRad-rad*rad));
                
        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2.0);
        arg2 = x*endRad*sin(theta)*sqrt(1.0-tt*tt);
        
        if(arg2 == 0) {
                be = 0.5;
        } else {
                be = NR_BessJ1(arg2)/arg2;
        }
        val = cos(arg1)*(1.0-tt*tt)*be;
        
        return(val);
}

double PolyCoreBicelle(double dp[], double q)
{
        int i;
        int nord = 20;
        double scale, length, sigma, bkg, radius, radthick, facthick;
        double rhoc, rhoh, rhor, rhosolv;
        double answer, Vpoly;
        double Pi,lolim,uplim,summ,yyy,rad,AR,Rsqr,Rsqrsumm,Rsqryyy;
        
        scale = dp[0];
        radius = dp[1];
        sigma = dp[2];                          //sigma is the standard mean deviation
        length = dp[3];
        radthick = dp[4];
        facthick= dp[5];
        rhoc = dp[6];
        rhoh = dp[7];
        rhor=dp[8];
        rhosolv = dp[9];
        bkg = dp[10];
        
        Pi = 4.0*atan(1.0);
        
        lolim = exp(log(radius)-(4.*sigma));
        if (lolim<0.0) {
                lolim=0.0;              //to avoid numerical error when  va<0 (-ve r value)
        }
        uplim = exp(log(radius)+(4.*sigma));
        
        summ = 0.0;
        Rsqrsumm = 0.0;
        
        for(i=0;i<nord;i++) {
                rad = ( Gauss20Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                AR=(1.0/(rad*sigma*sqrt(2.0*Pi)))*exp(-(0.5*((log(radius/rad))/sigma)*((log(radius/rad))/sigma)));
                yyy = AR* Gauss20Wt[i] * BicelleIntegration(q,rad,radthick,facthick,rhoc,rhoh,rhor,rhosolv,length);
                Rsqryyy= Gauss20Wt[i] * AR * (rad+radthick)*(rad+radthick);             //SRK normalize to total dimensions
                summ += yyy;
                Rsqrsumm += Rsqryyy;
        }
        
        answer = (uplim-lolim)/2.0*summ;
        Rsqr = (uplim-lolim)/2.0*Rsqrsumm;
        //normalize by average cylinder volume
        Vpoly = Pi*Rsqr*(length+2*facthick);
        answer /= Vpoly;
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
                
        return answer;
        
}

double
BicelleIntegration(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length){

        double answer,halfheight,Pi;
        double lolim,uplim,summ,yyy,zi;
        int nord,i;
        
        // set up the integration end points 
        Pi = 4.0*atan(1.0);
        nord = 76;
        lolim = 0.0;
        uplim = Pi/2;
        halfheight = length/2.0;
        
        summ = 0.0;                             // initialize integral
        i=0;
        for(i=0;i<nord;i++) {
                zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0;
                yyy = Gauss76Wt[i] * BicelleKernel(qq, rad, radthick, facthick, rhoc, rhoh, rhor,rhosolv, halfheight, zi);
                summ += yyy;
        }
        
        // calculate value of integral to return
        answer = (uplim-lolim)/2.0*summ;
        return(answer); 
}

double
BicelleKernel(double qq, double rad, double radthick, double facthick, double rhoc, double rhoh, double rhor, double rhosolv, double length, double dum)
{
        double dr1,dr2,dr3;
        double besarg1,besarg2;
        double vol1,vol2,vol3;
        double sinarg1,sinarg2;
        double t1,t2,t3;
        double retval,si1,si2,be1,be2;
        
        double Pi = 4.0*atan(1.0);
        
        dr1 = rhoc-rhoh;
        dr2 = rhor-rhosolv;
        dr3=  rhoh-rhor;
        vol1 = Pi*rad*rad*(2.0*length);
        vol2 = Pi*(rad+radthick)*(rad+radthick)*(2.0*length+2.0*facthick);
        vol3= Pi*(rad)*(rad)*(2.0*length+2.0*facthick);
        besarg1 = qq*rad*sin(dum);
        besarg2 = qq*(rad+radthick)*sin(dum);
        sinarg1 = qq*length*cos(dum);
        sinarg2 = qq*(length+facthick)*cos(dum);
        
        if(besarg1 == 0) {
                be1 = 0.5;
        } else {
                be1 = NR_BessJ1(besarg1)/besarg1;
        }
        if(besarg2 == 0) {
                be2 = 0.5;
        } else {
                be2 = NR_BessJ1(besarg2)/besarg2;
        }       
        if(sinarg1 == 0) {
                si1 = 1.0;
        } else {
                si1 = sin(sinarg1)/sinarg1;
        }
        if(sinarg2 == 0) {
                si2 = 1.0;
        } else {
                si2 = sin(sinarg2)/sinarg2;
        }
        t1 = 2.0*vol1*dr1*si1*be1;
        t2 = 2.0*vol2*dr2*si2*be2;
        t3 = 2.0*vol3*dr3*si2*be1;
        
        retval = ((t1+t2+t3)*(t1+t2+t3))*sin(dum);
        return(retval);
        
}


double
CSPPKernel(double dp[], double mu, double uu)
{       
        double aa,bb,cc, ta,tb,tc; 
        double Vin,Vot,V1,V2;
        double rhoA,rhoB,rhoC, rhoP, rhosolv;
        double dr0, drA,drB, drC;
        double arg1,arg2,arg3,arg4,t1,t2, t3, t4;
        double Pi,retVal;

        aa = dp[1];
        bb = dp[2];
        cc = dp[3];
        ta = dp[4];
        tb = dp[5];
        tc = dp[6];
        rhoA=dp[7];
        rhoB=dp[8];
        rhoC=dp[9];
        rhoP=dp[10];
        rhosolv=dp[11];
        dr0=rhoP-rhosolv;
        drA=rhoA-rhosolv;
        drB=rhoB-rhosolv;
        drC=rhoC-rhosolv; 
        Vin=(aa*bb*cc);
        Vot=(aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc);
        V1=(2.0*ta*bb*cc);   //  incorrect V1 (aa*bb*cc+2*ta*bb*cc)
        V2=(2.0*aa*tb*cc);  // incorrect V2(aa*bb*cc+2*aa*tb*cc)
        aa = aa/bb;
        ta=(aa+2.0*ta)/bb;
        tb=(aa+2.0*tb)/bb;
        
        Pi = 4.0*atan(1.0);
        
        arg1 = (mu*aa/2.0)*sin(Pi*uu/2.0);
        arg2 = (mu/2.0)*cos(Pi*uu/2.0);
        arg3=  (mu*ta/2.0)*sin(Pi*uu/2.0);
        arg4=  (mu*tb/2.0)*cos(Pi*uu/2.0);
                         
        if(arg1==0.0){
                t1 = 1.0;
        } else {
                t1 = (sin(arg1)/arg1);                //defn for CSPP model sin(arg1)/arg1    test:  (sin(arg1)/arg1)*(sin(arg1)/arg1)   
        }
        if(arg2==0.0){
                t2 = 1.0;
        } else {
                t2 = (sin(arg2)/arg2);           //defn for CSPP model sin(arg2)/arg2   test: (sin(arg2)/arg2)*(sin(arg2)/arg2)    
        }       
        if(arg3==0.0){
                t3 = 1.0;
        } else {
                t3 = sin(arg3)/arg3;
        }
        if(arg4==0.0){
                t4 = 1.0;
        } else {
                t4 = sin(arg4)/arg4;
        }
        retVal =( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 )*( dr0*t1*t2*Vin + drA*(t3-t1)*t2*V1+ drB*t1*(t4-t2)*V2 );   //  correct FF : square of sum of phase factors
        return(retVal); 

}

/*      CSParallelepiped  :  calculates the form factor of a Parallelepiped with a core-shell structure
 -- different SLDs can be used for the face and rim

Uses 76 pt Gaussian quadrature for both integrals
*/
double
CSParallelepiped(double dp[], double q)
{
        int i,j;
        double scale,aa,bb,cc,ta,tb,tc,rhoA,rhoB,rhoC,rhoP,rhosolv,bkg;         //local variables of coefficient wave
        int nordi=76;                   //order of integration
        int nordj=76;
        double va,vb;           //upper and lower integration limits
        double summ,yyy,answer;                 //running tally of integration
        double summj,vaj,vbj;                   //for the inner integration
        double mu,mudum,arg,sigma,uu,vol;
        
        
        //      Pi = 4.0*atan(1.0);
        va = 0.0;
        vb = 1.0;               //orintational average, outer integral
        vaj = 0.0;
        vbj = 1.0;              //endpoints of inner integral
        
        summ = 0.0;                     //initialize intergral
        
        scale = dp[0];
        aa = dp[1];
        bb = dp[2];
        cc = dp[3];
        ta  = dp[4];
        tb  = dp[5];
        tc  = dp[6];   // is 0 at the moment  
        rhoA=dp[7];   //rim A SLD
        rhoB=dp[8];   //rim B SLD
        rhoC=dp[9];    //rim C SLD
        rhoP = dp[10];   //Parallelpiped core SLD
        rhosolv=dp[11];  // Solvent SLD
        bkg = dp[12];
        
        mu = q*bb;
        vol = aa*bb*cc+2.0*ta*bb*cc+2.0*aa*tb*cc+2.0*aa*bb*tc;          //calculate volume before rescaling
        
        // do the rescaling here, not in the kernel
        // normalize all WRT bb
        aa = aa/bb;
        cc = cc/bb;
        
        for(i=0;i<nordi;i++) {
                //setup inner integral over the ellipsoidal cross-section
                summj=0.0;
                sigma = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0;          //the outer dummy
                
                for(j=0;j<nordj;j++) {
                        //76 gauss points for the inner integral
                        uu = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;         //the inner dummy
                        mudum = mu*sqrt(1.0-sigma*sigma);
                        yyy = Gauss76Wt[j] * CSPPKernel(dp,mudum,uu);
                        summj += yyy;
                }
                //now calculate the value of the inner integral
                answer = (vbj-vaj)/2.0*summj;
                
                //finish the outer integral cc already scaled
                arg = mu*cc*sigma/2.0;
                if ( arg == 0.0 ) {
                        answer *= 1.0;
                } else {
                        answer *= sin(arg)*sin(arg)/arg/arg;
                }
                
                //now sum up the outer integral
                yyy = Gauss76Wt[i] * answer;
                summ += yyy;
        }               //final scaling is done at the end of the function, after the NT_FP64 case
        
        answer = (vb-va)/2.0*summ;

        //normalize by volume
        answer /= vol;
        //convert to [cm-1]
        answer *= 1.0e8;
        //Scale
        answer *= scale;
        // add in the background
        answer += bkg;
        
        return answer;
}


/*
 * 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;
}