[fdb1487] | 1 | |
---|
| 2 | static double |
---|
| 3 | f_exp(double q, double r, double sld_in, double sld_out, |
---|
| 4 | double thickness, double A) |
---|
| 5 | { |
---|
[ce896fd] | 6 | const double vol = M_4PI_3 * cube(r); |
---|
[fdb1487] | 7 | const double qr = q * r; |
---|
| 8 | const double alpha = A * r/thickness; |
---|
| 9 | const double bes = sph_j1c(qr); |
---|
| 10 | const double B = (sld_out - sld_in)/expm1(A); |
---|
| 11 | const double C = sld_in - B; |
---|
| 12 | double fun; |
---|
| 13 | if (qr == 0.0) { |
---|
| 14 | fun = 1.0; |
---|
| 15 | } else if (fabs(A) > 0.0) { |
---|
| 16 | const double qrsq = qr*qr; |
---|
| 17 | const double alphasq = alpha*alpha; |
---|
| 18 | const double sumsq = alphasq + qrsq; |
---|
| 19 | double sinqr, cosqr; |
---|
| 20 | SINCOS(qr, sinqr, cosqr); |
---|
[ce896fd] | 21 | fun = -3.0*( |
---|
[fdb1487] | 22 | ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq |
---|
| 23 | - (alpha*sinqr/qr - cosqr) |
---|
| 24 | ) / sumsq; |
---|
| 25 | } else { |
---|
| 26 | fun = bes; |
---|
| 27 | } |
---|
| 28 | return vol * (B*fun + C*bes); |
---|
| 29 | } |
---|
| 30 | |
---|
| 31 | static double |
---|
| 32 | f_linear(double q, double r, double sld, double slope) |
---|
| 33 | { |
---|
[ce896fd] | 34 | const double vol = M_4PI_3 * cube(r); |
---|
[fdb1487] | 35 | const double qr = q * r; |
---|
| 36 | const double bes = sph_j1c(qr); |
---|
| 37 | double fun = 0.0; |
---|
| 38 | if (qr > 0.0) { |
---|
| 39 | const double qrsq = qr*qr; |
---|
| 40 | double sinqr, cosqr; |
---|
| 41 | SINCOS(qr, sinqr, cosqr); |
---|
| 42 | // Jae-He's code seems to simplify to this |
---|
| 43 | // fun = 3.0 * slope * r * (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq); |
---|
| 44 | // Rederiving the math, we get the following instead: |
---|
| 45 | fun = 3.0 * slope * r * (2.0*cosqr + qr*sinqr)/(qrsq*qrsq); |
---|
| 46 | } |
---|
| 47 | return vol * (sld*bes + fun); |
---|
| 48 | } |
---|
| 49 | |
---|
| 50 | static double |
---|
| 51 | f_constant(double q, double r, double sld) |
---|
| 52 | { |
---|
| 53 | const double bes = sph_j1c(q * r); |
---|
[ce896fd] | 54 | const double vol = M_4PI_3 * cube(r); |
---|
[fdb1487] | 55 | return sld * vol * bes; |
---|
| 56 | } |
---|
| 57 | |
---|
| 58 | static double |
---|
| 59 | form_volume(double core_radius, double n, double thickness[]) |
---|
| 60 | { |
---|
| 61 | int i; |
---|
| 62 | double r = core_radius; |
---|
| 63 | for (i=0; i < n; i++) { |
---|
| 64 | r += thickness[i]; |
---|
| 65 | } |
---|
[ce896fd] | 66 | return M_4PI_3*cube(r); |
---|
[fdb1487] | 67 | } |
---|
| 68 | |
---|
| 69 | static double |
---|
[ce896fd] | 70 | Iq(double q, double sld_core, double core_radius, double sld_solvent, |
---|
| 71 | double n, double sld_in[], double sld_out[], double thickness[], |
---|
[fdb1487] | 72 | double A[]) |
---|
| 73 | { |
---|
| 74 | int i; |
---|
[ce896fd] | 75 | double r = core_radius; |
---|
| 76 | double f = f_constant(q, r, sld_core); |
---|
[fdb1487] | 77 | for (i=0; i<n; i++){ |
---|
| 78 | const double r0 = r; |
---|
| 79 | r += thickness[i]; |
---|
| 80 | if (r == r0) { |
---|
| 81 | // no thickness, so nothing to add |
---|
[abdd01c] | 82 | } else if (fabs(A[i]) < 1.0e-16 || sld_out[i] == sld_in[i]) { |
---|
[fdb1487] | 83 | f -= f_constant(q, r0, sld_in[i]); |
---|
| 84 | f += f_constant(q, r, sld_in[i]); |
---|
[abdd01c] | 85 | } else if (fabs(A[i]) < 1.0e-4) { |
---|
[fdb1487] | 86 | const double slope = (sld_out[i] - sld_in[i])/thickness[i]; |
---|
| 87 | f -= f_linear(q, r0, sld_in[i], slope); |
---|
| 88 | f += f_linear(q, r, sld_out[i], slope); |
---|
| 89 | } else { |
---|
| 90 | f -= f_exp(q, r0, sld_in[i], sld_out[i], thickness[i], A[i]); |
---|
| 91 | f += f_exp(q, r, sld_in[i], sld_out[i], thickness[i], A[i]); |
---|
| 92 | } |
---|
| 93 | } |
---|
[ce896fd] | 94 | f -= f_constant(q, r, sld_solvent); |
---|
| 95 | const double f2 = f * f * 1.0e-4; |
---|
[fdb1487] | 96 | |
---|
| 97 | return f2; |
---|
| 98 | } |
---|