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

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since ee60aa7 was ee60aa7, checked in by Paul Kienzle <pkienzle@…>, 6 years ago

clean up effective radius functions; improve mono_gauss_coil accuracy; start moving VR into C

  • Property mode set to 100644
File size: 2.1 KB
RevLine 
[fdb1487]1
2static double
3f_exp(double q, double r, double sld_in, double sld_out,
[d119f34]4    double thickness, double A, double side)
[fdb1487]5{
[ce896fd]6  const double vol = M_4PI_3 * cube(r);
[fdb1487]7  const double qr = q * r;
[925ad6e]8  const double bes = sas_3j1x_x(qr);
[d119f34]9  const double alpha = A * r/thickness;
10  double result;
[fdb1487]11  if (qr == 0.0) {
[d119f34]12    result = 1.0;
[fdb1487]13  } else if (fabs(A) > 0.0) {
[d119f34]14    const double qrsq = qr * qr;
15    const double alphasq = alpha * alpha;
[fdb1487]16    const double sumsq = alphasq + qrsq;
17    double sinqr, cosqr;
18    SINCOS(qr, sinqr, cosqr);
[d119f34]19    const double t1 = (alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr;
20    const double t2 = alpha*sinqr/qr - cosqr;
21    const double fun = -3.0*(t1/sumsq - t2)/sumsq;
22    const double slope = (sld_out - sld_in)/expm1(A);
23    const double contrast = slope*exp(A*side);
24    result = contrast*fun + (sld_in-slope)*bes;
[fdb1487]25  } else {
[d119f34]26    result = sld_in*bes;
[fdb1487]27  }
[d119f34]28  return vol * result;
[fdb1487]29}
30
31static double
[ee60aa7]32outer_radius(double radius_core, double n_shells, double thickness[])
[fdb1487]33{
[768c0c4]34  int n = (int)(n_shells+0.5);
[9762341]35  double r = radius_core;
[768c0c4]36  for (int i=0; i < n; i++) {
[fdb1487]37    r += thickness[i];
38  }
[ee60aa7]39  return r;
40}
41
42static double
43form_volume(double radius_core, double n_shells, double thickness[])
44{
45  return M_4PI_3*cube(outer_radius(radius_core, n_shells, thickness));
[fdb1487]46}
47
[d277229]48static double
49effective_radius(int mode, double radius_core, double n_shells, double thickness[])
50{
[ee60aa7]51  // case 1: outer radius
52  return outer_radius(radius_core, n_shells, thickness);
[d277229]53}
54
[71b751d]55static void
56Fq(double q, double *F1, double *F2, double sld_core, double radius_core, double sld_solvent,
[d119f34]57    double n_shells, double sld_in[], double sld_out[], double thickness[],
[fdb1487]58    double A[])
59{
[d119f34]60  int n = (int)(n_shells+0.5);
[9762341]61  double r_out = radius_core;
[d119f34]62  double f = f_exp(q, r_out, sld_core, 0.0, 0.0, 0.0, 0.0);
63  for (int i=0; i < n; i++){
64    const double r_in = r_out;
65    r_out += thickness[i];
66    f -= f_exp(q, r_in, sld_in[i], sld_out[i], thickness[i], A[i], 0.0);
67    f += f_exp(q, r_out, sld_in[i], sld_out[i], thickness[i], A[i], 1.0);
[fdb1487]68  }
[d119f34]69  f -= f_exp(q, r_out, sld_solvent, 0.0, 0.0, 0.0, 0.0);
[fdb1487]70
[71b751d]71  *F1 = 1e-2 * f;
72  *F2 = 1e-4 * f * f;
[fdb1487]73}
Note: See TracBrowser for help on using the repository browser.