[5d4777d] | 1 | real form_volume(real radius, real cap_radius, real length); |
---|
| 2 | real Iq(real q, real sld, real solvent_sld, |
---|
| 3 | real radius, real cap_radius, real length); |
---|
| 4 | real Iqxy(real qx, real qy, real sld, real solvent_sld, |
---|
| 5 | real radius, real cap_radius, real length, real theta, real phi); |
---|
| 6 | |
---|
| 7 | // Integral over a convex lens kernel for t in [h/R,1]. See the docs for |
---|
| 8 | // the definition of the function being integrated. |
---|
| 9 | // q is the magnitude of the q vector. |
---|
| 10 | // h is the length of the lens "inside" the cylinder. This negative wrt the |
---|
| 11 | // definition of h in the docs. |
---|
| 12 | // cap_radius is the radius of the lens |
---|
| 13 | // length is the cylinder length, or the separation between the lens halves |
---|
| 14 | // alpha is the angle of the cylinder wrt q. |
---|
| 15 | real _cap_kernel(real q, real h, real cap_radius, real length, |
---|
| 16 | real sin_alpha, real cos_alpha); |
---|
| 17 | real _cap_kernel(real q, real h, real cap_radius, real length, |
---|
| 18 | real sin_alpha, real cos_alpha) |
---|
| 19 | { |
---|
| 20 | // For speed, we are pre-calculating terms which are constant over the |
---|
| 21 | // kernel. |
---|
| 22 | const real upper = REAL(1.0); |
---|
| 23 | const real lower = h/cap_radius; // integral lower bound |
---|
| 24 | const real m = q*cos_alpha*cap_radius; // cos argument slope |
---|
| 25 | const real b = q*cos_alpha*(REAL(0.5)*length-h); // cos argument intercept |
---|
| 26 | const real qrst = q*sin_alpha*cap_radius; // Q*R*sin(theta) |
---|
| 27 | real total = REAL(0.0); |
---|
| 28 | for (int i=0; i<76 ;i++) { |
---|
| 29 | // translate a point in [-1,1] to a point in [lower,upper] |
---|
| 30 | //const real t = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; |
---|
| 31 | const real t = REAL(0.5)*(Gauss76Z[i]*(upper-lower)+upper+lower); |
---|
| 32 | const real radical = REAL(1.0) - t*t; |
---|
| 33 | const real caparg = qrst*sqrt(radical); // cap bessel function arg |
---|
| 34 | const real be = (caparg == REAL(0.0) ? REAL(0.5) : J1(caparg)/caparg); |
---|
| 35 | const real Fq = cos(m*t + b) * radical * be; |
---|
| 36 | total += Gauss76Wt[i] * Fq; |
---|
| 37 | } |
---|
| 38 | // translate dx in [-1,1] to dx in [lower,upper] |
---|
| 39 | //const real form = (upper-lower)/2.0*total; |
---|
| 40 | const real integral = REAL(0.5)*(upper-lower)*total; |
---|
| 41 | return REAL(4.0)*M_PI*cap_radius*cap_radius*cap_radius*integral; |
---|
| 42 | } |
---|
| 43 | |
---|
| 44 | real form_volume(real radius, real cap_radius, real length) |
---|
| 45 | { |
---|
| 46 | // cap radius should never be less than radius when this is called |
---|
| 47 | // Note: cap volume = pi hc/6 * (3 a^2 + hc^2), where a is the cylinder |
---|
| 48 | // radius and hc is the height of the cap. Multiply by two for both ends. |
---|
| 49 | // So: |
---|
| 50 | // cap V = pi hc (r^2 + hc^2/3) |
---|
| 51 | // cylinder V = pi r^2 L |
---|
| 52 | // V = cylinder V + cap V |
---|
| 53 | // = pi r^2 L + pi hc (r^2 + hc^2/3) |
---|
| 54 | // = pi * (r^2 (L+hc) + hc^3/3) |
---|
| 55 | const real hc = cap_radius - sqrt(cap_radius*cap_radius - radius*radius); |
---|
| 56 | return M_PI*(radius*radius*(length+hc) + REAL(0.333333333333333)*hc*hc*hc); |
---|
| 57 | } |
---|
| 58 | |
---|
| 59 | real Iq(real q, |
---|
| 60 | real sld, |
---|
| 61 | real solvent_sld, |
---|
| 62 | real radius, |
---|
| 63 | real cap_radius, |
---|
| 64 | real length) |
---|
| 65 | { |
---|
| 66 | real sn, cn; // slots to hold sincos function output |
---|
| 67 | |
---|
| 68 | // Exclude invalid inputs. |
---|
| 69 | if (cap_radius < radius) return REAL(-1.0); |
---|
| 70 | |
---|
| 71 | const real lower = REAL(0.0); |
---|
| 72 | const real upper = M_PI_2; |
---|
| 73 | const real h = sqrt(cap_radius*cap_radius - radius*radius); // negative h |
---|
| 74 | real total = REAL(0.0); |
---|
| 75 | for (int i=0; i<76 ;i++) { |
---|
| 76 | // translate a point in [-1,1] to a point in [lower,upper] |
---|
| 77 | const real alpha= REAL(0.5)*(Gauss76Z[i]*(upper-lower) + upper + lower); |
---|
| 78 | SINCOS(alpha, sn, cn); |
---|
| 79 | |
---|
| 80 | const real cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn); |
---|
| 81 | |
---|
| 82 | // The following is CylKernel() / sin(alpha), but we are doing it in place |
---|
| 83 | // to avoid sin(alpha)/sin(alpha) for alpha = 0. It is also a teensy bit |
---|
| 84 | // faster since we don't multiply and divide sin(alpha). |
---|
| 85 | const real besarg = q*radius*sn; |
---|
| 86 | const real siarg = REAL(0.5)*q*length*cn; |
---|
| 87 | // lim_{x->0} J1(x)/x = 1/2, lim_{x->0} sin(x)/x = 1 |
---|
| 88 | const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); |
---|
| 89 | const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg); |
---|
| 90 | const real cyl_Fq = M_PI*radius*radius*length*REAL(2.0)*bj*si; |
---|
| 91 | |
---|
| 92 | // Volume weighted average F(q) |
---|
| 93 | const real Aq = cyl_Fq + cap_Fq; |
---|
| 94 | total += Gauss76Wt[i] * Aq * Aq * sn; // sn for spherical coord integration |
---|
| 95 | } |
---|
| 96 | // translate dx in [-1,1] to dx in [lower,upper] |
---|
| 97 | const real form = total * REAL(0.5)*(upper-lower); |
---|
| 98 | |
---|
| 99 | // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 |
---|
| 100 | // NOTE that for this (Fournet) definition of the integral, one must MULTIPLY by Vcyl |
---|
| 101 | // The additional volume factor is for polydisperse volume normalization. |
---|
| 102 | const real s = (sld - solvent_sld); |
---|
| 103 | return REAL(1.0e-4) * form * s * s; // form_volume(radius, cap_radius, length); |
---|
| 104 | } |
---|
| 105 | |
---|
| 106 | |
---|
| 107 | real Iqxy(real qx, real qy, |
---|
| 108 | real sld, |
---|
| 109 | real solvent_sld, |
---|
| 110 | real radius, |
---|
| 111 | real cap_radius, |
---|
| 112 | real length, |
---|
| 113 | real theta, |
---|
| 114 | real phi) |
---|
| 115 | { |
---|
| 116 | real sn, cn; // slots to hold sincos function output |
---|
| 117 | |
---|
| 118 | // Exclude invalid inputs. |
---|
| 119 | if (cap_radius < radius) return REAL(-1.0); |
---|
| 120 | |
---|
| 121 | // Compute angle alpha between q and the cylinder axis |
---|
| 122 | SINCOS(theta*M_PI_180, sn, cn); |
---|
| 123 | // # The following correction factor exists in sasview, but it can't be |
---|
| 124 | // # right, so we are leaving it out for now. |
---|
| 125 | const real q = sqrt(qx*qx+qy*qy); |
---|
| 126 | const real cos_val = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); |
---|
| 127 | const real alpha = acos(cos_val); // rod angle relative to q |
---|
| 128 | SINCOS(alpha, sn, cn); |
---|
| 129 | |
---|
| 130 | const real h = sqrt(cap_radius*cap_radius - radius*radius); // negative h |
---|
| 131 | const real cap_Fq = _cap_kernel(q, h, cap_radius, length, sn, cn); |
---|
| 132 | |
---|
| 133 | // The following is CylKernel() / sin(alpha), but we are doing it in place |
---|
| 134 | // to avoid sin(alpha)/sin(alpha) for alpha = 0. It is also a teensy bit |
---|
| 135 | // faster since we don't multiply and divide sin(alpha). |
---|
| 136 | const real besarg = q*radius*sn; |
---|
| 137 | const real siarg = REAL(0.5)*q*length*cn; |
---|
| 138 | // lim_{x->0} J1(x)/x = 1/2, lim_{x->0} sin(x)/x = 1 |
---|
| 139 | const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); |
---|
| 140 | const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg); |
---|
| 141 | const real cyl_Fq = M_PI*radius*radius*length*REAL(2.0)*bj*si; |
---|
| 142 | |
---|
| 143 | // Volume weighted average F(q) |
---|
| 144 | const real Aq = cyl_Fq + cap_Fq; |
---|
| 145 | |
---|
| 146 | // Multiply by contrast^2, normalize by cylinder volume and convert to cm-1 |
---|
| 147 | const real s = (sld - solvent_sld); |
---|
| 148 | return REAL(1.0e-4) * Aq * Aq * s * s; // form_volume(radius, cap_radius, length); |
---|
| 149 | } |
---|