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 | } |
---|