source: sasmodels/sasmodels/models/pringle.c @ 1596de3

double form_volume(double radius, double thickness, double alpha, double beta);

double Iq(double q,
          double radius,
          double thickness,
          double alpha,
          double beta,
          double sld_pringle,
          double sld_solvent);


static
void _integrate_bessel(
    double radius,
    double alpha,
    double beta,
    double q_sin_psi,
    double q_cos_psi,
    double n,
    double *Sn,
    double *Cn)
{
    // translate gauss point z in [-1,1] to a point in [0, radius]
    const double zm = 0.5*radius;
    const double zb = 0.5*radius;

    // evaluate at Gauss points
    double sumS = 0.0;          // initialize integral
    double sumC = 0.0;          // initialize integral
    double r;
    for (int i=0; i < 76; i++) {
        r = Gauss76Z[i]*zm + zb;

        const double qrs = r*q_sin_psi;
        const double qrrc = r*r*q_cos_psi;

        double y = Gauss76Wt[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs);
        double S, C;
        SINCOS(alpha*qrrc, S, C);
        sumS += y*S;
        sumC += y*C;
    }

    #if 1
    // TODO: should the normalization be to R^2 or the last value, r^2
    // the gaussian window does not go all the way from 0 to 1.
    //radius = Gauss76Z[75] * zm + zb;
    *Sn = zm*sumS / (r*r);
    *Cn = zm*sumC / (r*r);
    #else
    *Sn = zm*sumS / (radius*radius);
    *Cn = zm*sumC / (radius*radius);
    #endif
}

static
double _sum_bessel_orders(
    double radius,
    double alpha,
    double beta,
    double q_sin_psi,
    double q_cos_psi)
{
    //calculate sum term from n = -3 to 3
    //Note 1:
    //    S_n(-x) = (-1)^S_n(x)
    //    => S_n^2(-x) = S_n^2(x),
    //    => sum_-k^k Sk = S_0^2 + 2*sum_1^kSk^2
    //Note 2:
    //    better precision to sum terms from smaller to larger
    //    though it doesn't seem to make a difference in this case.
    double Sn, Cn, sum;
    sum = 0.0;
    for (int n=3; n>0; n--) {
      _integrate_bessel(radius, alpha, beta, q_sin_psi, q_cos_psi, n, &Sn, &Cn);
      sum += 2.0*(Sn*Sn + Cn*Cn);
    }
    _integrate_bessel(radius, alpha, beta, q_sin_psi, q_cos_psi, 0, &Sn, &Cn);
    sum += Sn*Sn+ Cn*Cn;
    return sum;
}

static
double _integrate_psi(
    double q,
    double radius,
    double thickness,
    double alpha,
    double beta)
{
    // translate gauss point z in [-1,1] to a point in [0, pi/2]
    const double zm = M_PI_4;
    const double zb = M_PI_4;

    double sum = 0.0;
    for (int i = 0; i < 76; i++) {
        double psi = Gauss76Z[i]*zm + zb;
        double sin_psi, cos_psi;
        SINCOS(psi, sin_psi, cos_psi);
        double bessel_term = _sum_bessel_orders(radius, alpha, beta, q*sin_psi, q*cos_psi);
        double sinc_term = square(sinc(q * thickness * cos_psi / 2.0));
        double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term;
        sum += Gauss76Wt[i] * pringle_kernel;
    }

    return zm * sum;
}

double form_volume(double radius, double thickness, double alpha, double beta)
{
    // TODO: Normalize by form volume
    //return M_PI*radius*radius*thickness;
    return 1.0;
}

double Iq(
    double q,
    double radius,
    double thickness,
    double alpha,
    double beta,
    double sld_pringle,
    double sld_solvent)
{
    double form = _integrate_psi(q, radius, thickness, alpha, beta);
    double contrast = sld_pringle - sld_solvent;
    double volume = M_PI*radius*radius*thickness;
    // TODO: If normalize by form volume, need an extra volume here
    //return 1.0e-4*form * square(contrast * volume);
    return 1.0e-4*form * square(contrast) * volume;
}