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