[fdb1487] | 1 | |
---|
| 2 | static double |
---|
| 3 | f_exp(double q, double r, double sld_in, double sld_out, |
---|
| 4 | double thickness, double A) |
---|
| 5 | { |
---|
| 6 | const double vol = 4.0/3.0 * M_PI * r * r * r; |
---|
| 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); |
---|
| 21 | fun = -3.0( |
---|
| 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 | { |
---|
| 34 | const double vol = 4.0/3.0 * M_PI * r * r * r; |
---|
| 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); |
---|
| 54 | const double vol = 4.0/3.0 * M_PI * r * r * r; |
---|
| 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 | } |
---|
| 66 | return 4.0/3.0 * M_PI * r * r * r; |
---|
| 67 | } |
---|
| 68 | |
---|
| 69 | static double |
---|
| 70 | Iq(double q, double core_sld, double core_radius, double solvent_sld, |
---|
| 71 | double n, double in_sld[], double out_sld[], double thickness[], |
---|
| 72 | double A[]) |
---|
| 73 | { |
---|
| 74 | int i; |
---|
| 75 | r = core_radius; |
---|
| 76 | f = f_constant(q, r, core_sld); |
---|
| 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 |
---|
| 82 | } else if (fabs(A[i]) < 1e-16 || sld_out[i] == sld_in[i]) { |
---|
| 83 | f -= f_constant(q, r0, sld_in[i]); |
---|
| 84 | f += f_constant(q, r, sld_in[i]); |
---|
| 85 | } else if (fabs(A[i]) < 1e-4) { |
---|
| 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 | } |
---|
| 94 | f -= f_constant(q, r, solvent_sld); |
---|
| 95 | f2 = f * f * 1.0e-4; |
---|
| 96 | |
---|
| 97 | return f2; |
---|
| 98 | } |
---|