source: sasmodels/sasmodels/models/onion.c @ d277229

static double
f_exp(double q, double r, double sld_in, double sld_out,
    double thickness, double A, double side)
{
  const double vol = M_4PI_3 * cube(r);
  const double qr = q * r;
  const double bes = sas_3j1x_x(qr);
  const double alpha = A * r/thickness;
  double result;
  if (qr == 0.0) {
    result = 1.0;
  } else if (fabs(A) > 0.0) {
    const double qrsq = qr * qr;
    const double alphasq = alpha * alpha;
    const double sumsq = alphasq + qrsq;
    double sinqr, cosqr;
    SINCOS(qr, sinqr, cosqr);
    const double t1 = (alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr;
    const double t2 = alpha*sinqr/qr - cosqr;
    const double fun = -3.0*(t1/sumsq - t2)/sumsq;
    const double slope = (sld_out - sld_in)/expm1(A);
    const double contrast = slope*exp(A*side);
    result = contrast*fun + (sld_in-slope)*bes;
  } else {
    result = sld_in*bes;
  }
  return vol * result;
}

static double
form_volume(double radius_core, double n_shells, double thickness[])
{
  int n = (int)(n_shells+0.5);
  double r = radius_core;
  for (int i=0; i < n; i++) {
    r += thickness[i];
  }
  return M_4PI_3*cube(r);
}

static double
effective_radius(int mode, double radius_core, double n_shells, double thickness[])
{
  int n = (int)(n_shells+0.5);
  double r = radius_core;
  for (int i=0; i < n; i++) {
    r += thickness[i];
  }
  return r;
}

static void
Fq(double q, double *F1, double *F2, double sld_core, double radius_core, double sld_solvent,
    double n_shells, double sld_in[], double sld_out[], double thickness[],
    double A[])
{
  int n = (int)(n_shells+0.5);
  double r_out = radius_core;
  double f = f_exp(q, r_out, sld_core, 0.0, 0.0, 0.0, 0.0);
  for (int i=0; i < n; i++){
    const double r_in = r_out;
    r_out += thickness[i];
    f -= f_exp(q, r_in, sld_in[i], sld_out[i], thickness[i], A[i], 0.0);
    f += f_exp(q, r_out, sld_in[i], sld_out[i], thickness[i], A[i], 1.0);
  }
  f -= f_exp(q, r_out, sld_solvent, 0.0, 0.0, 0.0, 0.0);

  *F1 = 1e-2 * f;
  *F2 = 1e-4 * f * f;
}