Changes in / [07c8d46:fe8ff99] in sasmodels
- Location:
- sasmodels
- Files:
-
- 4 added
- 2 deleted
- 74 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/generate.py
r7fcdc9f r1e7b0db0 139 139 square(x) = x*x 140 140 cube(x) = x*x*x 141 s inc(x) = sin(x)/x, with sin(0)/0 -> 1141 sas_sinx_x(x) = sin(x)/x, with sin(0)/0 -> 1 142 142 all double precision constants must include the decimal point 143 143 all double declarations may be converted to half, float, or long double -
sasmodels/kernel_header.c
r218cdbc r1e7b0db0 146 146 inline double square(double x) { return x*x; } 147 147 inline double cube(double x) { return x*x*x; } 148 inline double s inc(double x) { return x==0 ? 1.0 : sin(x)/x; }148 inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; } 149 149 150 #if 0 150 #if 1 151 //think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 151 152 #define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 152 153 SINCOS(phi*M_PI_180, sn, cn); \ 153 154 q = sqrt(qx*qx + qy*qy); \ 154 cn = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * cos(theta*M_PI_180)); \155 cn = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180)); \ 155 156 sn = sqrt(1 - cn*cn); \ 156 157 } while (0) -
sasmodels/kernel_template.c
r0d6e865 r1e7b0db0 133 133 inline double square(double x) { return x*x; } 134 134 inline double cube(double x) { return x*x*x; } 135 inline double s inc(double x) { return x==0 ? 1.0 : sin(x)/x; }135 inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; } 136 136 137 137 -
sasmodels/models/barbell.c
r3a48772 r592343f 33 33 const double t = Gauss76Z[i]*zm + zb; 34 34 const double radical = 1.0 - t*t; 35 const double bj = sas_ J1c(qrst*sqrt(radical));35 const double bj = sas_2J1x_x(qrst*sqrt(radical)); 36 36 const double Fq = cos(m*t + b) * radical * bj; 37 37 total += Gauss76Wt[i] * Fq; … … 49 49 { 50 50 const double bell_fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 51 const double bj = sas_ J1c(q*radius*sin_alpha);52 const double si = s inc(q*half_length*cos_alpha);51 const double bj = sas_2J1x_x(q*radius*sin_alpha); 52 const double si = sas_sinx_x(q*half_length*cos_alpha); 53 53 const double cyl_fq = 2.0*M_PI*radius*radius*half_length*bj*si; 54 54 const double Aq = bell_fq + cyl_fq; -
sasmodels/models/barbell.py
r0d6e865 rfcb33e4 20 20 .. math:: 21 21 22 I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q )\right>22 I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q,\alpha).sin(\alpha)\right> 23 23 24 where the amplitude $A(q )$ is given as24 where the amplitude $A(q,\alpha)$ with the rod axis at angle $\alpha$ to $q$ is given as 25 25 26 26 .. math:: 27 27 28 28 A(q) =&\ \pi r^2L 29 \frac{\sin\left(\tfrac12 qL\cos\ theta\right)}30 {\tfrac12 qL\cos\ theta}31 \frac{2 J_1(qr\sin\ theta)}{qr\sin\theta} \\29 \frac{\sin\left(\tfrac12 qL\cos\alpha\right)} 30 {\tfrac12 qL\cos\alpha} 31 \frac{2 J_1(qr\sin\alpha)}{qr\sin\alpha} \\ 32 32 &\ + 4 \pi R^3 \int_{-h/R}^1 dt 33 \cos\left[ q\cos\ theta33 \cos\left[ q\cos\alpha 34 34 \left(Rt + h + {\tfrac12} L\right)\right] 35 35 \times (1-t^2) 36 \frac{J_1\left[qR\sin\ theta \left(1-t^2\right)^{1/2}\right]}37 {qR\sin\ theta \left(1-t^2\right)^{1/2}}36 \frac{J_1\left[qR\sin\alpha \left(1-t^2\right)^{1/2}\right]} 37 {qR\sin\alpha \left(1-t^2\right)^{1/2}} 38 38 39 39 The $\left<\ldots\right>$ brackets denote an average of the structure over 40 all orientations. $\left<A^2(q )\right>$ is then the form factor, $P(q)$.40 all orientations. $\left<A^2(q,\alpha)\right>$ is then the form factor, $P(q)$. 41 41 The scale factor is equivalent to the volume fraction of cylinders, each of 42 42 volume, $V$. Contrast $\Delta\rho$ is the difference of scattering length … … 85 85 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 86 86 * **Last Modified by:** Paul Butler **Date:** March 20, 2016 87 * **Last Reviewed by:** Paul Butler **Date:** March 20, 201687 * **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 88 88 """ 89 89 from numpy import inf -
sasmodels/models/bcc_paracrystal.py
rb0c4271 r925ad6e 128 128 # pylint: enable=bad-whitespace, line-too-long 129 129 130 source = ["lib/s ph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc_paracrystal.c"]130 source = ["lib/sas_3j1x_x.c", "lib/gauss150.c", "lib/sphere_form.c", "bcc_paracrystal.c"] 131 131 132 132 # parameters for demo -
sasmodels/models/binary_hard_sphere.c
r4f79d94 r925ad6e 58 58 qr2 = r2*q; 59 59 60 sc1 = s ph_j1c(qr1);61 sc2 = s ph_j1c(qr2);60 sc1 = sas_3j1x_x(qr1); 61 sc2 = sas_3j1x_x(qr2); 62 62 b1 = r1*r1*r1*(rho1-rhos)*M_4PI_3*sc1; 63 63 b2 = r2*r2*r2*(rho2-rhos)*M_4PI_3*sc2; -
sasmodels/models/binary_hard_sphere.py
rb0c4271 r925ad6e 108 108 ] 109 109 110 source = ["lib/s ph_j1c.c", "binary_hard_sphere.c"]110 source = ["lib/sas_3j1x_x.c", "binary_hard_sphere.c"] 111 111 112 112 # parameters for demo and documentation -
sasmodels/models/capped_cylinder.c
r3a48772 r592343f 39 39 const double t = Gauss76Z[i]*zm + zb; 40 40 const double radical = 1.0 - t*t; 41 const double bj = sas_ J1c(qrst*sqrt(radical));41 const double bj = sas_2J1x_x(qrst*sqrt(radical)); 42 42 const double Fq = cos(m*t + b) * radical * bj; 43 43 total += Gauss76Wt[i] * Fq; … … 54 54 { 55 55 const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 56 const double bj = sas_ J1c(q*radius*sin_alpha);57 const double si = s inc(q*half_length*cos_alpha);56 const double bj = sas_2J1x_x(q*radius*sin_alpha); 57 const double si = sas_sinx_x(q*half_length*cos_alpha); 58 58 const double cyl_Fq = 2.0*M_PI*radius*radius*half_length*bj*si; 59 59 const double Aq = cap_Fq + cyl_Fq; -
sasmodels/models/capped_cylinder.py
r0d6e865 rfcb33e4 21 21 .. math:: 22 22 23 I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q )\right>23 I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q,\alpha).sin(\alpha)\right> 24 24 25 where the amplitude $A(q )$ is given as25 where the amplitude $A(q,\alpha)$ with the rod axis at angle $\alpha$ to $q$ is given as 26 26 27 27 .. math:: 28 28 29 29 A(q) =&\ \pi r^2L 30 \frac{\sin\left(\tfrac12 qL\cos\ theta\right)}31 {\tfrac12 qL\cos\ theta}32 \frac{2 J_1(qr\sin\ theta)}{qr\sin\theta} \\30 \frac{\sin\left(\tfrac12 qL\cos\alpha\right)} 31 {\tfrac12 qL\cos\alpha} 32 \frac{2 J_1(qr\sin\alpha)}{qr\sin\alpha} \\ 33 33 &\ + 4 \pi R^3 \int_{-h/R}^1 dt 34 \cos\left[ q\cos\ theta34 \cos\left[ q\cos\alpha 35 35 \left(Rt + h + {\tfrac12} L\right)\right] 36 36 \times (1-t^2) 37 \frac{J_1\left[qR\sin\ theta \left(1-t^2\right)^{1/2}\right]}38 {qR\sin\ theta \left(1-t^2\right)^{1/2}}37 \frac{J_1\left[qR\sin\alpha \left(1-t^2\right)^{1/2}\right]} 38 {qR\sin\alpha \left(1-t^2\right)^{1/2}} 39 39 40 40 The $\left<\ldots\right>$ brackets denote an average of the structure over … … 88 88 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 89 89 * **Last Modified by:** Paul Butler **Date:** September 30, 2016 90 * **Last Reviewed by:** Richard Heenan **Date:** March 19, 2016 90 * **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 91 91 92 """ 92 93 from numpy import inf -
sasmodels/models/core_multi_shell.c
rc5ac2b2 r925ad6e 3 3 f_constant(double q, double r, double sld) 4 4 { 5 const double bes = s ph_j1c(q * r);5 const double bes = sas_3j1x_x(q * r); 6 6 const double vol = M_4PI_3 * cube(r); 7 7 return sld * vol * bes; … … 33 33 f = 0.; 34 34 for (int i=0; i<n; i++) { 35 f += M_4PI_3 * cube(r) * (sld[i] - last_sld) * s ph_j1c(q*r);35 f += M_4PI_3 * cube(r) * (sld[i] - last_sld) * sas_3j1x_x(q*r); 36 36 last_sld = sld[i]; 37 37 r += thickness[i]; 38 38 } 39 f += M_4PI_3 * cube(r) * (solvent_sld - last_sld) * s ph_j1c(q*r);39 f += M_4PI_3 * cube(r) * (solvent_sld - last_sld) * sas_3j1x_x(q*r); 40 40 return f * f * 1.0e-4; 41 41 } -
sasmodels/models/core_multi_shell.py
r2d73a53 r925ad6e 101 101 ] 102 102 103 source = ["lib/s ph_j1c.c", "core_multi_shell.c"]103 source = ["lib/sas_3j1x_x.c", "core_multi_shell.c"] 104 104 105 105 def profile(sld_core, radius, sld_solvent, n, sld, thickness): -
sasmodels/models/core_shell_bicelle.c
r5bddd89 r592343f 55 55 double sinarg2 = qq*(length+facthick)*cos_alpha; 56 56 57 be1 = sas_ J1c(besarg1);58 be2 = sas_ J1c(besarg2);59 si1 = s inc(sinarg1);60 si2 = s inc(sinarg2);57 be1 = sas_2J1x_x(besarg1); 58 be2 = sas_2J1x_x(besarg2); 59 si1 = sas_sinx_x(sinarg1); 60 si2 = sas_sinx_x(sinarg2); 61 61 62 62 const double t = vol1*dr1*si1*be1 + -
sasmodels/models/core_shell_bicelle.py
rb829b16 r8afefae 41 41 42 42 I(Q,\alpha) = \frac{\text{scale}}{V_t} \cdot 43 F(Q,\alpha)^2 + \text{background}43 F(Q,\alpha)^2.sin(\alpha) + \text{background} 44 44 45 45 where … … 85 85 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 86 86 * **Last Modified by:** Paul Butler **Date:** September 30, 2016 87 * **Last Reviewed by:** Richard Heenan **Date:** October 5, 201687 * **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 88 88 """ 89 89 … … 141 141 # pylint: enable=bad-whitespace, line-too-long 142 142 143 source = ["lib/ Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c",143 source = ["lib/sas_Si.c", "lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", 144 144 "core_shell_bicelle.c"] 145 145 … … 156 156 phi=0) 157 157 158 qx, qy = 0.4 * cos(90), 0.5 * sin(0)158 #qx, qy = 0.4 * cos(pi/2.0), 0.5 * sin(0) -
sasmodels/models/core_shell_cylinder.c
r9aa4881 r592343f 11 11 double _cyl(double vd, double besarg, double siarg) 12 12 { 13 return vd * s inc(siarg) * sas_J1c(besarg);13 return vd * sas_sinx_x(siarg) * sas_2J1x_x(besarg); 14 14 } 15 15 -
sasmodels/models/core_shell_cylinder.py
r755ecc2 r8e68ea0 9 9 .. math:: 10 10 11 I(q,\alpha) = \frac{\text{scale}}{V_s} F^2(q ) + \text{background}11 I(q,\alpha) = \frac{\text{scale}}{V_s} F^2(q,\alpha).sin(\alpha) + \text{background} 12 12 13 13 where … … 15 15 .. math:: 16 16 17 F(q ) = &\ (\rho_c - \rho_s) V_c17 F(q,\alpha) = &\ (\rho_c - \rho_s) V_c 18 18 \frac{\sin \left( q \tfrac12 L\cos\alpha \right)} 19 19 {q \tfrac12 L\cos\alpha} -
sasmodels/models/core_shell_ellipsoid.py
r73e08ae r8e68ea0 137 137 # pylint: enable=bad-whitespace, line-too-long 138 138 139 source = ["lib/s ph_j1c.c", "lib/gfn.c", "lib/gauss76.c",139 source = ["lib/sas_3j1x_x.c", "lib/gfn.c", "lib/gauss76.c", 140 140 "core_shell_ellipsoid.c"] 141 141 … … 162 162 163 163 q = 0.1 164 phi = pi/6 165 qx = q*cos(p hi)166 qy = q*sin(p hi)167 # After redefinition of angles find new reasonable values for unit test168 #tests = [169 ## Accuracy tests based on content in test/utest_coreshellellipsoidXTmodel.py170 #[{'radius_equat_core': 200.0,171 #'x_core': 0.1,172 #'thick_shell': 50.0,173 #'x_polar_shell': 0.2,174 #'sld_core': 2.0,175 #'sld_shell': 1.0,176 #'sld_solvent': 6.3,177 #'background': 0.001,178 #'scale': 1.0,179 #}, 1.0, 0.00189402],164 # tests had in old coords theta=0, phi=0; new coords theta=90, phi=0 165 qx = q*cos(pi/6.0) 166 qy = q*sin(pi/6.0) 167 # 11Jan2017 RKH sorted tests after redefinition of angles 168 tests = [ 169 # Accuracy tests based on content in test/utest_coreshellellipsoidXTmodel.py 170 [{'radius_equat_core': 200.0, 171 'x_core': 0.1, 172 'thick_shell': 50.0, 173 'x_polar_shell': 0.2, 174 'sld_core': 2.0, 175 'sld_shell': 1.0, 176 'sld_solvent': 6.3, 177 'background': 0.001, 178 'scale': 1.0, 179 }, 1.0, 0.00189402], 180 180 181 181 # Additional tests with larger range of parameters 182 # [{'background': 0.01}, 0.1, 11.6915], 183 184 # [{'radius_equat_core': 20.0, 185 # 'x_core': 200.0, 186 # 'thick_shell': 54.0, 187 # 'x_polar_shell': 3.0, 188 # 'sld_core': 20.0, 189 # 'sld_shell': 10.0, 190 # 'sld_solvent': 6.0, 191 # 'background': 0.0, 192 # 'scale': 1.0, 193 # }, 0.01, 8688.53], 194 195 # [{'background': 0.001}, (0.4, 0.5), 0.00690673], 196 197 # [{'radius_equat_core': 20.0, 198 # 'x_core': 200.0, 199 # 'thick_shell': 54.0, 200 # 'x_polar_shell': 3.0, 201 # 'sld_core': 20.0, 202 # 'sld_shell': 10.0, 203 # 'sld_solvent': 6.0, 204 # 'background': 0.01, 205 # 'scale': 0.01, 206 # }, (qx, qy), 0.0100002], 207 # ] 182 [{'background': 0.01}, 0.1, 11.6915], 183 184 [{'radius_equat_core': 20.0, 185 'x_core': 200.0, 186 'thick_shell': 54.0, 187 'x_polar_shell': 3.0, 188 'sld_core': 20.0, 189 'sld_shell': 10.0, 190 'sld_solvent': 6.0, 191 'background': 0.0, 192 'scale': 1.0, 193 }, 0.01, 8688.53], 194 195 # 2D tests 196 [{'background': 0.001, 197 'theta': 90.0, 198 'phi': 0.0, 199 }, (0.4, 0.5), 0.00690673], 200 201 [{'radius_equat_core': 20.0, 202 'x_core': 200.0, 203 'thick_shell': 54.0, 204 'x_polar_shell': 3.0, 205 'sld_core': 20.0, 206 'sld_shell': 10.0, 207 'sld_solvent': 6.0, 208 'background': 0.01, 209 'scale': 0.01, 210 'theta': 90.0, 211 'phi': 0.0, 212 }, (qx, qy), 0.01000025], 213 ] -
sasmodels/models/core_shell_parallelepiped.c
r14838a3 r1e7b0db0 87 87 double sin_uu, cos_uu; 88 88 SINCOS(M_PI_2*uu, sin_uu, cos_uu); 89 const double si1 = s inc(mu_proj * sin_uu * a_scaled);90 const double si2 = s inc(mu_proj * cos_uu);91 const double si3 = s inc(mu_proj * sin_uu * ta);92 const double si4 = s inc(mu_proj * cos_uu * tb);89 const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 90 const double si2 = sas_sinx_x(mu_proj * cos_uu); 91 const double si3 = sas_sinx_x(mu_proj * sin_uu * ta); 92 const double si4 = sas_sinx_x(mu_proj * cos_uu * tb); 93 93 94 94 // Expression in libCylinder.c (neither drC nor Vot are used) … … 109 109 110 110 // now sum up the outer integral 111 const double si = s inc(mu * c_scaled * sigma);111 const double si = sas_sinx_x(mu * c_scaled * sigma); 112 112 outer_total += Gauss76Wt[i] * inner_total * si * si; 113 113 } … … 160 160 double tc = length_a + 2.0*thick_rim_c; 161 161 //handle arg=0 separately, as sin(t)/t -> 1 as t->0 162 double siA = s inc(0.5*q*length_a*cos_val_a);163 double siB = s inc(0.5*q*length_b*cos_val_b);164 double siC = s inc(0.5*q*length_c*cos_val_c);165 double siAt = s inc(0.5*q*ta*cos_val_a);166 double siBt = s inc(0.5*q*tb*cos_val_b);167 double siCt = s inc(0.5*q*tc*cos_val_c);162 double siA = sas_sinx_x(0.5*q*length_a*cos_val_a); 163 double siB = sas_sinx_x(0.5*q*length_b*cos_val_b); 164 double siC = sas_sinx_x(0.5*q*length_c*cos_val_c); 165 double siAt = sas_sinx_x(0.5*q*ta*cos_val_a); 166 double siBt = sas_sinx_x(0.5*q*tb*cos_val_b); 167 double siCt = sas_sinx_x(0.5*q*tc*cos_val_c); 168 168 169 169 -
sasmodels/models/core_shell_sphere.py
r4962519 r925ad6e 75 75 # pylint: enable=bad-whitespace, line-too-long 76 76 77 source = ["lib/s ph_j1c.c", "lib/core_shell.c", "core_shell_sphere.c"]77 source = ["lib/sas_3j1x_x.c", "lib/core_shell.c", "core_shell_sphere.c"] 78 78 79 79 demo = dict(scale=1, background=0, radius=60, thickness=10, -
sasmodels/models/cylinder.c
rb829b16 r592343f 18 18 const double qr = q*radius; 19 19 const double qh = q*0.5*length; 20 return sas_ J1c(qr*sn) * sinc(qh*cn);20 return sas_2J1x_x(qr*sn) * sas_sinx_x(qh*cn); 21 21 } 22 22 -
sasmodels/models/cylinder.py
r4cdd0cc rb7e8b94 2 2 # Note: model title and parameter table are inserted automatically 3 3 r""" 4 The form factor is normalized by the particle volume V = \piR^2L. 4 5 5 For information about polarised and magnetic scattering, see 6 6 the :ref:`magnetism` documentation. … … 14 14 .. math:: 15 15 16 P(q,\alpha) = \frac{\text{scale}}{V} F^2(q,\alpha) + \text{background}16 P(q,\alpha) = \frac{\text{scale}}{V} F^2(q,\alpha).sin(\alpha) + \text{background} 17 17 18 18 where … … 25 25 \frac{J_1 \left(q R \sin \alpha\right)}{q R \sin \alpha} 26 26 27 and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, $V $27 and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, $V =\pi R^2L$ 28 28 is the volume of the cylinder, $L$ is the length of the cylinder, $R$ is the 29 29 radius of the cylinder, and $\Delta\rho$ (contrast) is the scattering length … … 35 35 .. math:: 36 36 37 F^2(q)=\int_{0}^{\pi/2}{F^2(q,\ theta)\sin(\theta)d\theta}37 F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}=\int_{0}^{1}{F^2(q,u)du} 38 38 39 39 40 To provide easy access to the orientation of the cylinder, we define the 41 axis of the cylinder using two angles $\theta$ and $\phi$. Those angles 40 Numerical integration is simplified by a change of variable to $u = cos(\alpha)$ with 41 $sin(\alpha)=\sqrt{1-u^2}$. 42 43 The output of the 1D scattering intensity function for randomly oriented 44 cylinders is thus given by 45 46 .. math:: 47 48 P(q) = \frac{\text{scale}}{V} 49 \int_0^{\pi/2} F^2(q,\alpha) \sin \alpha\ d\alpha + \text{background} 50 51 52 NB: The 2nd virial coefficient of the cylinder is calculated based on the 53 radius and length values, and used as the effective radius for $S(q)$ 54 when $P(q) \cdot S(q)$ is applied. 55 56 For oriented cylinders, we define the direction of the 57 axis of the cylinder using two angles $\theta$ (note this is not the 58 same as the scattering angle used in q) and $\phi$. Those angles 42 59 are defined in :numref:`cylinder-angle-definition` . 43 60 … … 48 65 Definition of the angles for oriented cylinders. 49 66 50 51 NB: The 2nd virial coefficient of the cylinder is calculated based on the 52 radius and length values, and used as the effective radius for $S(q)$ 53 when $P(q) \cdot S(q)$ is applied. 54 55 The output of the 1D scattering intensity function for randomly oriented 56 cylinders is then given by 57 58 .. math:: 59 60 P(q) = \frac{\text{scale}}{V} 61 \int_0^{\pi/2} F^2(q,\alpha) \sin \alpha\ d\alpha + \text{background} 62 63 The $\theta$ and $\phi$ parameters are not used for the 1D output. 67 The $\theta$ and $\phi$ parameters only appear in the model when fitting 2d data. 64 68 65 69 Validation … … 74 78 75 79 P(q) = \int_0^{\pi/2} d\phi 76 \int_0^\pi p(\ alpha) P_0(q,\alpha) \sin \alpha\ d\alpha80 \int_0^\pi p(\theta) P_0(q,\theta) \sin \theta\ d\theta 77 81 78 82 79 83 where $p(\theta,\phi) = 1$ is the probability distribution for the orientation 80 and $P_0(q,\ alpha)$ is the scattering intensity for the fully oriented84 and $P_0(q,\theta)$ is the scattering intensity for the fully oriented 81 85 system, and then comparing to the 1D result. 82 86 … … 145 149 146 150 qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 147 # After redefinition of angles, find new tests values 148 #tests = [[{}, 0.2, 0.042761386790780453], 149 # [{}, [0.2], [0.042761386790780453]], 150 # [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852], 151 # [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]], 152 # ] 151 # After redefinition of angles, find new tests values. Was 10 10 in old coords 152 tests = [[{}, 0.2, 0.042761386790780453], 153 [{}, [0.2], [0.042761386790780453]], 154 # new coords 155 [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852], 156 [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]], 157 # old coords [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852], 158 # [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]], 159 ] 153 160 del qx, qy # not necessary to delete, but cleaner 154 161 # ADDED by: RKH ON: 18Mar2016 renamed sld's etc -
sasmodels/models/ellipsoid.c
r73e08ae r925ad6e 18 18 const double r = radius_equatorial 19 19 * sqrt(1.0 + sin_alpha*sin_alpha*(ratio*ratio - 1.0)); 20 const double f = s ph_j1c(q*r);20 const double f = sas_3j1x_x(q*r); 21 21 22 22 return f*f; -
sasmodels/models/ellipsoid.py
r0d6e865 r925ad6e 135 135 ] 136 136 137 source = ["lib/s ph_j1c.c", "lib/gauss76.c", "ellipsoid.c"]137 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"] 138 138 139 139 def ER(radius_polar, radius_equatorial): -
sasmodels/models/elliptical_cylinder.c
r251f54b r592343f 39 39 const double theta = ( Gauss20Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 40 40 const double r = sin_val*sqrt(rA - rB*cos(theta)); 41 const double be = sas_ J1c(q*r);41 const double be = sas_2J1x_x(q*r); 42 42 inner_sum += Gauss20Wt[j] * be * be; 43 43 } … … 46 46 47 47 //now calculate outer integral 48 const double si = s inc(q*0.5*length*cos_val);48 const double si = sas_sinx_x(q*0.5*length*cos_val); 49 49 outer_sum += Gauss76Wt[i] * inner_sum * si * si; 50 50 } … … 73 73 // Given: radius_major = r_ratio * radius_minor 74 74 const double r = radius_minor*sqrt(square(r_ratio*cos_nu) + cos_mu*cos_mu); 75 const double be = sas_ J1c(q*r);76 const double si = s inc(q*0.5*length*cos_val);75 const double be = sas_2J1x_x(q*r); 76 const double si = sas_sinx_x(q*0.5*length*cos_val); 77 77 const double Aq = be * si; 78 78 const double delrho = sld - solvent_sld; -
sasmodels/models/elliptical_cylinder.py
ra807206 rfcb33e4 20 20 21 21 I(\vec q)=\frac{1}{V_\text{cyl}}\int{d\psi}\int{d\phi}\int{ 22 p(\theta,\phi,\psi)F^2(\vec q,\alpha,\psi)\sin(\ theta)d\theta}22 p(\theta,\phi,\psi)F^2(\vec q,\alpha,\psi)\sin(\alpha)d\alpha} 23 23 24 24 with the functions … … 26 26 .. math:: 27 27 28 F( \vecq,\alpha,\psi) = 2\frac{J_1(a)\sin(b)}{ab}28 F(q,\alpha,\psi) = 2\frac{J_1(a)\sin(b)}{ab} 29 29 30 30 where … … 32 32 .. math:: 33 33 34 a &= \vec q\sin(\alpha)\left[ 35 r^2_\text{major}\sin^2(\psi)+r^2_\text{minor}\cos(\psi) \right]^{1/2} 34 a = qr'\sin(\alpha) 35 36 b = q\frac{L}{2}\cos(\alpha) 37 38 r'=\frac{r_{minor}}{\sqrt{2}}\sqrt{(1+\nu^{2}) + (1-\nu^{2})cos(\psi)} 36 39 37 b &= \vec q\frac{L}{2}\cos(\alpha)38 40 39 and the angle $\ Psi$ is defined as the orientation of the major axis of the41 and the angle $\psi$ is defined as the orientation of the major axis of the 40 42 ellipse with respect to the vector $\vec q$. The angle $\alpha$ is the angle 41 43 between the axis of the cylinder and $\vec q$. … … 95 97 96 98 L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and 97 Neutron Scattering*, Plenum, New York, (1987) 99 Neutron Scattering*, Plenum, New York, (1987) [see table 3.4] 100 101 Authorship and Verification 102 ---------------------------- 103 104 * **Author:** 105 * **Last Modified by:** 106 * **Last Reviewed by:** Richard Heenan - corrected equation in docs **Date:** December 21, 2016 107 98 108 """ 99 109 -
sasmodels/models/fcc_paracrystal.py
r0bef47b r925ad6e 116 116 # pylint: enable=bad-whitespace, line-too-long 117 117 118 source = ["lib/s ph_j1c.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc_paracrystal.c"]118 source = ["lib/sas_3j1x_x.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc_paracrystal.c"] 119 119 120 120 # parameters for demo -
sasmodels/models/flexible_cylinder.c
re6408d0 r592343f 14 14 { 15 15 const double contrast = sld - solvent_sld; 16 const double cross_section = sas_ J1c(q*radius);16 const double cross_section = sas_2J1x_x(q*radius); 17 17 const double volume = M_PI*radius*radius*length; 18 18 const double flex = Sk_WR(q, length, kuhn_length); -
sasmodels/models/flexible_cylinder_elliptical.c
r92ce163 r592343f 22 22 SINCOS(zi, sn, cn); 23 23 const double arg = q*sqrt(a*a*sn*sn + b*b*cn*cn); 24 const double yyy = sas_ J1c(arg);24 const double yyy = sas_2J1x_x(arg); 25 25 sum += Gauss76Wt[i] * yyy * yyy; 26 26 } -
sasmodels/models/fractal.c
r217590b r925ad6e 14 14 //calculate P(q) for the spherical subunits 15 15 const double V = M_4PI_3*cube(radius); 16 const double pq = V * square((sld_block-sld_solvent)*s ph_j1c(q*radius));16 const double pq = V * square((sld_block-sld_solvent)*sas_3j1x_x(q*radius)); 17 17 18 18 // scale to units cm-1 sr-1 (assuming data on absolute scale) -
sasmodels/models/fractal.py
rd1cfa86 r925ad6e 97 97 # pylint: enable=bad-whitespace, line-too-long 98 98 99 source = ["lib/s ph_j1c.c", "lib/sas_gamma.c", "lib/fractal_sq.c", "fractal.c"]99 source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "lib/fractal_sq.c", "fractal.c"] 100 100 101 101 demo = dict(volfraction=0.05, -
sasmodels/models/fractal_core_shell.py
rd6f60c3 r925ad6e 95 95 # pylint: enable=bad-whitespace, line-too-long 96 96 97 source = ["lib/s ph_j1c.c", "lib/sas_gamma.c", "lib/core_shell.c",97 source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "lib/core_shell.c", 98 98 "lib/fractal_sq.c", "fractal_core_shell.c"] 99 99 -
sasmodels/models/fuzzy_sphere.py
r3a48772 r925ad6e 81 81 # pylint: enable=bad-whitespace,line-too-long 82 82 83 source = ["lib/s ph_j1c.c"]83 source = ["lib/sas_3j1x_x.c"] 84 84 85 85 # No volume normalization despite having a volume parameter … … 91 91 Iq = """ 92 92 const double qr = q*radius; 93 const double bes = s ph_j1c(qr);93 const double bes = sas_3j1x_x(qr); 94 94 const double qf = q*fuzziness; 95 95 const double fq = bes * (sld - sld_solvent) * form_volume(radius) * exp(-0.5*qf*qf); -
sasmodels/models/guinier_porod.py
ra807206 rcdcebf1 4 4 and dimensionality of scattering objects, including asymmetric objects 5 5 such as rods or platelets, and shapes intermediate between spheres 6 and rods or between rods and platelets. 6 and rods or between rods and platelets, and overcomes some of the 7 deficiencies of the (Beaucage) Unified_Power_Rg model (see Hammouda, 2010). 7 8 8 9 Definition … … 59 60 --------- 60 61 61 A Guinier, G Fournet, *Small-Angle Scattering of X-Rays*, 62 John Wiley and Sons, New York, (1955) 62 B Hammouda, *A new Guinier-Porod model, J. Appl. Cryst.*, (2010), 43, 716-719 63 63 64 O Glatter, O Kratky, *Small-Angle X-Ray Scattering*, Academic Press (1982) 65 Check out Chapter 4 on Data Treatment, pages 155-156. 64 B Hammouda, *Analysis of the Beaucage model, J. Appl. Cryst.*, (2010), 43, 1474-1478 65 66 66 """ 67 67 -
sasmodels/models/hollow_cylinder.c
rf8f0991 r592343f 20 20 { 21 21 const double qs = q*sin_val; 22 const double lam1 = sas_ J1c((radius+thickness)*qs);23 const double lam2 = sas_ J1c(radius*qs);22 const double lam1 = sas_2J1x_x((radius+thickness)*qs); 23 const double lam2 = sas_2J1x_x(radius*qs); 24 24 const double gamma_sq = square(radius/(radius+thickness)); 25 //Note: lim_{thickness -> 0} psi = J0(radius*qs)26 //Note: lim_{radius -> 0} psi = sas_ J1c(thickness*qs)25 //Note: lim_{thickness -> 0} psi = sas_J0(radius*qs) 26 //Note: lim_{radius -> 0} psi = sas_2J1x_x(thickness*qs) 27 27 const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00 28 const double t2 = s inc(0.5*q*length*cos_val);28 const double t2 = sas_sinx_x(0.5*q*length*cos_val); 29 29 return psi*t2; 30 30 } -
sasmodels/models/hollow_rectangular_prism.c
r6f676fb r1e7b0db0 45 45 SINCOS(theta, sin_theta, cos_theta); 46 46 47 const double termC1 = s inc(q * c_half * cos(theta));48 const double termC2 = s inc(q * (c_half-thickness)*cos(theta));47 const double termC1 = sas_sinx_x(q * c_half * cos(theta)); 48 const double termC2 = sas_sinx_x(q * (c_half-thickness)*cos(theta)); 49 49 50 50 double inner_sum = 0.0; … … 57 57 // Amplitude AP from eqn. (13), rewritten to avoid round-off effects when arg=0 58 58 59 const double termA1 = s inc(q * a_half * sin_theta * sin_phi);60 const double termA2 = s inc(q * (a_half-thickness) * sin_theta * sin_phi);59 const double termA1 = sas_sinx_x(q * a_half * sin_theta * sin_phi); 60 const double termA2 = sas_sinx_x(q * (a_half-thickness) * sin_theta * sin_phi); 61 61 62 const double termB1 = s inc(q * b_half * sin_theta * cos_phi);63 const double termB2 = s inc(q * (b_half-thickness) * sin_theta * cos_phi);62 const double termB1 = sas_sinx_x(q * b_half * sin_theta * cos_phi); 63 const double termB2 = sas_sinx_x(q * (b_half-thickness) * sin_theta * cos_phi); 64 64 65 65 const double AP1 = vol_total * termA1 * termB1 * termC1; -
sasmodels/models/lib/core_shell.c
r3a48772 r925ad6e 16 16 const double core_qr = q * radius; 17 17 const double core_contrast = core_sld - shell_sld; 18 const double core_bes = s ph_j1c(core_qr);18 const double core_bes = sas_3j1x_x(core_qr); 19 19 const double core_volume = M_4PI_3 * cube(radius); 20 20 double f = core_volume * core_bes * core_contrast; … … 23 23 const double shell_qr = q * (radius + thickness); 24 24 const double shell_contrast = shell_sld - solvent_sld; 25 const double shell_bes = s ph_j1c(shell_qr);25 const double shell_bes = sas_3j1x_x(shell_qr); 26 26 const double shell_volume = M_4PI_3 * cube(radius + thickness); 27 27 f += shell_volume * shell_bes * shell_contrast; -
sasmodels/models/lib/gfn.c
r3a48772 r925ad6e 18 18 // changing to more accurate sph_j1c since the following inexplicably fails on Radeon Nano. 19 19 //const double siq = (uq == 0.0 ? 1.0 : 3.0*(sin(uq)/uq/uq - cos(uq)/uq)/uq); 20 const double siq = s ph_j1c(uq);20 const double siq = sas_3j1x_x(uq); 21 21 const double vc = M_4PI_3*aa*aa*bb; 22 22 const double gfnc = siq*vc*delpc; … … 26 26 const double vt = M_4PI_3*trmaj*trmaj*trmin; 27 27 //const double sit = (ut == 0.0 ? 1.0 : 3.0*(sin(ut)/ut/ut - cos(ut)/ut)/ut); 28 const double sit = s ph_j1c(ut);28 const double sit = sas_3j1x_x(ut); 29 29 const double gfnt = sit*vt*delps; 30 30 -
sasmodels/models/lib/polevl.c
r0278e3f r447e9aa 28 28 * N 0 29 29 * 30 * The function p1evl() assumes that coef[N]= 1.0 and is30 * The function p1evl() assumes that C_N = 1.0 and is 31 31 * omitted from the array. Its calling arguments are 32 32 * otherwise the same as polevl(). -
sasmodels/models/lib/sas_J0.c
r1596de3 rc8902ac 55 55 #if FLOAT_SIZE>4 56 56 //Cephes double precission 57 double j0(double x);57 double cephes_j0(double x); 58 58 59 59 constant double PPJ0[8] = { … … 147 147 }; 148 148 149 double j0(double x)149 double cephes_j0(double x) 150 150 { 151 151 double w, z, p, q, xn; … … 186 186 #else 187 187 //Cephes single precission 188 float j0f(float x);188 float cephes_j0f(float x); 189 189 190 190 constant float MOJ0[8] = { … … 221 221 }; 222 222 223 float j0f(float x)223 float cephes_j0f(float x) 224 224 { 225 225 float xx, w, z, p, q, xn; … … 257 257 258 258 #if FLOAT_SIZE>4 259 #define sas_J0 j0259 #define sas_J0 cephes_j0 260 260 #else 261 #define sas_J0 j0f261 #define sas_J0 cephes_j0f 262 262 #endif -
sasmodels/models/lib/sas_J1.c
r1596de3 r473a9f1 42 42 #if FLOAT_SIZE>4 43 43 //Cephes double pression function 44 double j1(double x);44 double cephes_j1(double x); 45 45 46 46 constant double RPJ1[8] = { … … 106 106 0.0 }; 107 107 108 double j1(double x)108 double cephes_j1(double x) 109 109 { 110 110 … … 144 144 #else 145 145 //Single precission version of cephes 146 float j1f(float x);146 float cephes_j1f(float x); 147 147 148 148 constant float JPJ1[8] = { … … 179 179 }; 180 180 181 float j1f(float x)181 float cephes_j1f(float x) 182 182 { 183 183 … … 211 211 212 212 #if FLOAT_SIZE>4 213 #define sas_J1 j1213 #define sas_J1 cephes_j1 214 214 #else 215 #define sas_J1 j1f215 #define sas_J1 cephes_j1f 216 216 #endif 217 217 218 218 //Finally J1c function that equals 2*J1(x)/x 219 double sas_ J1c(double x);220 double sas_ J1c(double x)219 double sas_2J1x_x(double x); 220 double sas_2J1x_x(double x) 221 221 { 222 222 return (x != 0.0 ) ? 2.0*sas_J1(x)/x : 1.0; -
sasmodels/models/lib/sas_JN.c
r1596de3 rc8902ac 50 50 #if FLOAT_SIZE > 4 51 51 52 double jn( int n, double x );53 double jn( int n, double x ) {52 double cephes_jn( int n, double x ); 53 double cephes_jn( int n, double x ) { 54 54 55 55 // PAK: seems to be machine epsilon/2 … … 75 75 76 76 if( n == 0 ) 77 return( sign * j0(x) );77 return( sign * cephes_j0(x) ); 78 78 if( n == 1 ) 79 return( sign * j1(x) );79 return( sign * cephes_j1(x) ); 80 80 if( n == 2 ) 81 return( sign * (2.0 * j1(x) / x -j0(x)) );81 return( sign * (2.0 * cephes_j1(x) / x - cephes_j0(x)) ); 82 82 83 83 if( x < MACHEP ) … … 112 112 113 113 if( fabs(pk) > fabs(pkm1) ) 114 ans = j1(x)/pk;114 ans = cephes_j1(x)/pk; 115 115 else 116 ans = j0(x)/pkm1;116 ans = cephes_j0(x)/pkm1; 117 117 118 118 return( sign * ans ); … … 121 121 #else 122 122 123 float jnf(int n, float x);124 float jnf(int n, float x)123 float cephes_jnf(int n, float x); 124 float cephes_jnf(int n, float x) 125 125 { 126 126 // PAK: seems to be machine epsilon/2 … … 146 146 147 147 if( n == 0 ) 148 return( sign * j0f(x) );148 return( sign * cephes_j0f(x) ); 149 149 if( n == 1 ) 150 return( sign * j1f(x) );150 return( sign * cephes_j1f(x) ); 151 151 if( n == 2 ) 152 return( sign * (2.0 * j1f(x) / x -j0f(x)) );152 return( sign * (2.0 * cephes_j1f(x) / x - cephes_j0f(x)) ); 153 153 154 154 if( x < MACHEP ) … … 189 189 190 190 if( r > ans ) /* if( fabs(pk) > fabs(pkm1) ) */ 191 ans = sign * j1f(x)/pk;191 ans = sign * cephes_j1f(x)/pk; 192 192 else 193 ans = sign * j0f(x)/pkm1;193 ans = sign * cephes_j0f(x)/pkm1; 194 194 return( ans ); 195 195 } … … 197 197 198 198 #if FLOAT_SIZE>4 199 #define sas_JN jn199 #define sas_JN cephes_jn 200 200 #else 201 #define sas_JN jnf201 #define sas_JN cephes_jnf 202 202 #endif -
sasmodels/models/lib/sphere_form.c
rba32cdd r925ad6e 9 9 double sphere_form(double q, double radius, double sld, double solvent_sld) 10 10 { 11 const double fq = sphere_volume(radius) * s ph_j1c(q*radius);11 const double fq = sphere_volume(radius) * sas_3j1x_x(q*radius); 12 12 const double contrast = (sld - solvent_sld); 13 13 return 1.0e-4*square(contrast * fq); -
sasmodels/models/linear_pearls.c
r4962519 r925ad6e 44 44 45 45 //sine functions of a pearl 46 double psi = s ph_j1c(q * radius);46 double psi = sas_3j1x_x(q * radius); 47 47 48 48 // N pearls contribution … … 50 50 n_contrib = num_pearls; 51 51 for(int num=1; num<=n_max; num++) { 52 n_contrib += (2.0*(num_pearls-num)*s inc(q*separation*num));52 n_contrib += (2.0*(num_pearls-num)*sas_sinx_x(q*separation*num)); 53 53 } 54 54 // form factor for num_pearls -
sasmodels/models/linear_pearls.py
r4962519 r925ad6e 63 63 single = False 64 64 65 source = ["lib/s ph_j1c.c", "linear_pearls.c"]65 source = ["lib/sas_3j1x_x.c", "linear_pearls.c"] 66 66 67 67 demo = dict(scale=1.0, background=0.0, -
sasmodels/models/mass_fractal.c
r6d96b66 r925ad6e 5 5 { 6 6 //calculate P(q) 7 const double pq = square(s ph_j1c(q*radius));7 const double pq = square(sas_3j1x_x(q*radius)); 8 8 9 9 //calculate S(q) -
sasmodels/models/mass_fractal.py
r6d96b66 r925ad6e 86 86 # pylint: enable=bad-whitespace, line-too-long 87 87 88 source = ["lib/s ph_j1c.c", "lib/sas_gamma.c", "mass_fractal.c"]88 source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "mass_fractal.c"] 89 89 90 90 demo = dict(scale=1, background=0, -
sasmodels/models/multilayer_vesicle.c
r3a48772 r925ad6e 20 20 // layer 1 21 21 voli = M_4PI_3*ri*ri*ri; 22 fval += voli*sldi*s ph_j1c(ri*q);22 fval += voli*sldi*sas_3j1x_x(ri*q); 23 23 24 24 ri += thick_shell; … … 26 26 // layer 2 27 27 voli = M_4PI_3*ri*ri*ri; 28 fval -= voli*sldi*s ph_j1c(ri*q);28 fval -= voli*sldi*sas_3j1x_x(ri*q); 29 29 30 30 //do 2 layers at a time -
sasmodels/models/multilayer_vesicle.py
r041bc75 r925ad6e 107 107 # pylint: enable=bad-whitespace, line-too-long 108 108 109 source = ["lib/s ph_j1c.c", "multilayer_vesicle.c"]109 source = ["lib/sas_3j1x_x.c", "multilayer_vesicle.c"] 110 110 111 111 polydispersity = ["radius", "n_pairs"] -
sasmodels/models/onion.c
r9762341 r925ad6e 6 6 const double vol = M_4PI_3 * cube(r); 7 7 const double qr = q * r; 8 const double bes = s ph_j1c(qr);8 const double bes = sas_3j1x_x(qr); 9 9 const double alpha = A * r/thickness; 10 10 double result; -
sasmodels/models/onion.py
r9762341 r925ad6e 314 314 # pylint: enable=bad-whitespace, line-too-long 315 315 316 source = ["lib/s ph_j1c.c", "onion.c"]316 source = ["lib/sas_3j1x_x.c", "onion.c"] 317 317 single = False 318 318 -
sasmodels/models/parallelepiped.c
r14838a3 r1e7b0db0 39 39 double sin_uu, cos_uu; 40 40 SINCOS(M_PI_2*uu, sin_uu, cos_uu); 41 const double si1 = s inc(mu_proj * sin_uu * a_scaled);42 const double si2 = s inc(mu_proj * cos_uu);41 const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 42 const double si2 = sas_sinx_x(mu_proj * cos_uu); 43 43 inner_total += Gauss76Wt[j] * square(si1 * si2); 44 44 } 45 45 inner_total *= 0.5; 46 46 47 const double si = s inc(mu * c_scaled * sigma);47 const double si = sas_sinx_x(mu * c_scaled * sigma); 48 48 outer_total += Gauss76Wt[i] * inner_total * si * si; 49 49 } … … 70 70 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a); 71 71 72 const double siA = s inc(0.5*q*length_a*cos_val_a);73 const double siB = s inc(0.5*q*length_b*cos_val_b);74 const double siC = s inc(0.5*q*length_c*cos_val_c);72 const double siA = sas_sinx_x(0.5*q*length_a*cos_val_a); 73 const double siB = sas_sinx_x(0.5*q*length_b*cos_val_b); 74 const double siC = sas_sinx_x(0.5*q*length_c*cos_val_c); 75 75 const double V = form_volume(length_a, length_b, length_c); 76 76 const double drho = (sld - solvent_sld); -
sasmodels/models/pearl_necklace.c
r2126131 r4b541ac 34 34 // But there is a 1/(1-sinc) term below which blows up so don't bother 35 35 const double q_edge = q * edge_sep; 36 const double beta = ( Si(q*(A_s-radius)) -Si(q*radius)) / q_edge;37 const double gamma = Si(q_edge) / q_edge;38 const double psi = s ph_j1c(q*radius);36 const double beta = (sas_Si(q*(A_s-radius)) - sas_Si(q*radius)) / q_edge; 37 const double gamma = sas_Si(q_edge) / q_edge; 38 const double psi = sas_3j1x_x(q*radius); 39 39 40 40 // Precomputed sinc terms 41 const double si = s inc(q*A_s);41 const double si = sas_sinx_x(q*A_s); 42 42 const double omsi = 1.0 - si; 43 43 const double pow_si = pow(si, num_pearls); … … 54 54 - 2.0 * (1.0 - pow_si/si)*beta*beta / (omsi*omsi) 55 55 + 2.0 * num_strings*beta*beta / omsi 56 + num_strings * (2.0*gamma - square(s inc(q_edge/2.0)))56 + num_strings * (2.0*gamma - square(sas_sinx_x(q_edge/2.0))) 57 57 ); 58 58 -
sasmodels/models/pearl_necklace.py
r2126131 r4b541ac 92 92 ] 93 93 94 source = ["lib/ Si.c", "lib/sph_j1c.c", "pearl_necklace.c"]94 source = ["lib/sas_Si.c", "lib/sas_3j1x_x.c", "pearl_necklace.c"] 95 95 single = False # use double precision unless told otherwise 96 96 -
sasmodels/models/polymer_micelle.c
rc3ebc71 r925ad6e 31 31 32 32 // Self-correlation term of the core 33 const double bes_core = s ph_j1c(q*radius_core);33 const double bes_core = sas_3j1x_x(q*radius_core); 34 34 const double term1 = square(n_aggreg*beta_core*bes_core); 35 35 … … 41 41 // Interference cross-term between core and chains 42 42 const double chain_ampl = (qrg2 == 0.0) ? 1.0 : -expm1(-qrg2)/qrg2; 43 const double bes_corona = s inc(q*(radius_core + d_penetration * rg));43 const double bes_corona = sas_sinx_x(q*(radius_core + d_penetration * rg)); 44 44 const double term3 = 2.0 * n_aggreg * n_aggreg * beta_core * beta_corona * 45 45 bes_core * chain_ampl * bes_corona; -
sasmodels/models/polymer_micelle.py
rbba9361 r925ad6e 104 104 single = False 105 105 106 source = ["lib/s ph_j1c.c", "polymer_micelle.c"]106 source = ["lib/sas_3j1x_x.c", "polymer_micelle.c"] 107 107 108 108 demo = dict(scale=1, background=0, -
sasmodels/models/pringle.c
r30fbe2e r1e7b0db0 91 91 SINCOS(psi, sin_psi, cos_psi); 92 92 double bessel_term = _sum_bessel_orders(radius, alpha, beta, q*sin_psi, q*cos_psi); 93 double sinc_term = square(s inc(q * thickness * cos_psi / 2.0));93 double sinc_term = square(sas_sinx_x(q * thickness * cos_psi / 2.0)); 94 94 double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term; 95 95 sum += Gauss76Wt[i] * pringle_kernel; -
sasmodels/models/raspberry.c
r2c74c11 r925ad6e 51 51 52 52 //Form factors for each particle 53 psiL = s ph_j1c(q*rL);54 psiS = s ph_j1c(q*rS);53 psiL = sas_3j1x_x(q*rL); 54 psiS = sas_3j1x_x(q*rS); 55 55 56 56 //Cross term between large and small particles 57 sfLS = psiL*psiS*s inc(q*(rL+deltaS*rS));57 sfLS = psiL*psiS*sas_sinx_x(q*(rL+deltaS*rS)); 58 58 //Cross term between small particles at the surface 59 sfSS = psiS*psiS*s inc(q*(rL+deltaS*rS))*sinc(q*(rL+deltaS*rS));59 sfSS = psiS*psiS*sas_sinx_x(q*(rL+deltaS*rS))*sas_sinx_x(q*(rL+deltaS*rS)); 60 60 61 61 //Large sphere form factor term -
sasmodels/models/raspberry.py
r40a87fa r8e68ea0 129 129 Ref: J. coll. inter. sci. (2010) vol. 343 (1) pp. 36-41.""" 130 130 category = "shape:sphere" 131 #single = False 131 132 132 133 133 # [ "name", "units", default, [lower, upper], "type", "description"], … … 152 152 ] 153 153 154 source = ["lib/s ph_j1c.c", "raspberry.c"]154 source = ["lib/sas_3j1x_x.c", "raspberry.c"] 155 155 156 156 # parameters for demo -
sasmodels/models/rectangular_prism.c
rab2aea8 r1e7b0db0 33 33 SINCOS(theta, sin_theta, cos_theta); 34 34 35 const double termC = s inc(q * c_half * cos_theta);35 const double termC = sas_sinx_x(q * c_half * cos_theta); 36 36 37 37 double inner_sum = 0.0; … … 42 42 43 43 // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0 44 const double termA = s inc(q * a_half * sin_theta * sin_phi);45 const double termB = s inc(q * b_half * sin_theta * cos_phi);44 const double termA = sas_sinx_x(q * a_half * sin_theta * sin_phi); 45 const double termB = sas_sinx_x(q * b_half * sin_theta * cos_phi); 46 46 const double AP = termA * termB * termC; 47 47 inner_sum += Gauss76Wt[j] * AP * AP; -
sasmodels/models/sc_paracrystal.py
r0bef47b r925ad6e 133 133 # pylint: enable=bad-whitespace, line-too-long 134 134 135 source = ["lib/s ph_j1c.c", "lib/sphere_form.c", "lib/gauss150.c", "sc_paracrystal.c"]135 source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c", "lib/gauss150.c", "sc_paracrystal.c"] 136 136 137 137 demo = dict(scale=1, background=0, -
sasmodels/models/sphere.py
r7e6bea81 r925ad6e 66 66 ] 67 67 68 source = ["lib/s ph_j1c.c", "lib/sphere_form.c"]68 source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c"] 69 69 70 70 # No volume normalization despite having a volume parameter -
sasmodels/models/spherical_sld.c
r54bcd4a r925ad6e 34 34 const double qr = q * r; 35 35 const double qrsq = qr * qr; 36 const double bes = s ph_j1c(qr);36 const double bes = sas_3j1x_x(qr); 37 37 double sinqr, cosqr; 38 38 SINCOS(qr, sinqr, cosqr); … … 60 60 61 61 // uniform shell; r=0 => r^3=0 => f=0, so works for core as well. 62 f -= M_4PI_3 * cube(r) * sld_l * s ph_j1c(q*r);62 f -= M_4PI_3 * cube(r) * sld_l * sas_3j1x_x(q*r); 63 63 r += thickness[shell]; 64 f += M_4PI_3 * cube(r) * sld_l * s ph_j1c(q*r);64 f += M_4PI_3 * cube(r) * sld_l * sas_3j1x_x(q*r); 65 65 66 66 // iterate over sub_shells in the interface … … 92 92 } 93 93 // add in solvent effect 94 f -= M_4PI_3 * cube(r) * sld_solvent * s ph_j1c(q*r);94 f -= M_4PI_3 * cube(r) * sld_solvent * sas_3j1x_x(q*r); 95 95 96 96 const double f2 = f * f * 1.0e-4; -
sasmodels/models/spherical_sld.py
r2d65d51 r925ad6e 214 214 ] 215 215 # pylint: enable=bad-whitespace, line-too-long 216 source = ["lib/polevl.c", "lib/sas_erf.c", "lib/s ph_j1c.c", "spherical_sld.c"]216 source = ["lib/polevl.c", "lib/sas_erf.c", "lib/sas_3j1x_x.c", "spherical_sld.c"] 217 217 single = False # TODO: fix low q behaviour 218 218 -
sasmodels/models/stacked_disks.c
r98ce141 r6c3e266 53 53 const double sinarg2 = q*(halfheight+thick_layer)*cos_alpha; 54 54 55 const double be1 = sas_ J1c(besarg1);56 //const double be2 = sas_ J1c(besarg2);55 const double be1 = sas_2J1x_x(besarg1); 56 //const double be2 = sas_2J1x_x(besarg2); 57 57 const double be2 = be1; 58 const double si1 = s inc(sinarg1);59 const double si2 = s inc(sinarg2);58 const double si1 = sas_sinx_x(sinarg1); 59 const double si2 = sas_sinx_x(sinarg2); 60 60 61 61 const double dr1 = core_sld - solvent_sld; -
sasmodels/models/stacked_disks.py
r07300ea rb7e8b94 148 148 sld_layer=0.0, 149 149 sld_solvent=5.0, 150 theta= 0,150 theta=90, 151 151 phi=0) 152 #After redefinition to spherical coordinates find new reasonable test values 153 #tests = [ 154 # # Accuracy tests based on content in test/utest_extra_models.py. 155 # # Added 2 tests with n_stacked = 5 using SasView 3.1.2 - PDB 156 # [{'thick_core': 10.0, 157 # 'thick_layer': 15.0, 158 # 'radius': 3000.0, 159 # 'n_stacking': 1.0, 160 # 'sigma_d': 0.0, 161 # 'sld_core': 4.0, 162 # 'sld_layer': -0.4, 163 # 'solvent_sd': 5.0, 164 # 'theta': 0.0, 165 # 'phi': 0.0, 166 # 'scale': 0.01, 167 # 'background': 0.001, 168 # }, 0.001, 5075.12], 169 170 # [{'thick_core': 10.0, 171 # 'thick_layer': 15.0, 172 # 'radius': 3000.0, 173 # 'n_stacking': 5.0, 174 # 'sigma_d': 0.0, 175 # 'sld_core': 4.0, 176 # 'sld_layer': -0.4, 177 # 'solvent_sd': 5.0, 178 # 'theta': 0.0, 179 # 'phi': 0.0, 180 # 'scale': 0.01, 181 # 'background': 0.001, 182 # }, 0.001, 5065.12793824], 183 184 # [{'thick_core': 10.0, 185 # 'thick_layer': 15.0, 186 # 'radius': 3000.0, 187 # 'n_stacking': 5.0, 188 # 'sigma_d': 0.0, 189 # 'sld_core': 4.0, 190 # 'sld_layer': -0.4, 191 # 'solvent_sd': 5.0, 192 # 'theta': 0.0, 193 # 'phi': 0.0, 194 # 'scale': 0.01, 195 # 'background': 0.001, 196 # }, 0.164, 0.0127673597265], 197 198 # [{'thick_core': 10.0, 199 # 'thick_layer': 15.0, 200 # 'radius': 3000.0, 201 # 'n_stacking': 1.0, 202 # 'sigma_d': 0.0, 203 # 'sld_core': 4.0, 204 # 'sld_layer': -0.4, 205 # 'solvent_sd': 5.0, 206 # 'theta': 0.0, 207 # 'phi': 0.0, 208 # 'scale': 0.01, 209 # 'background': 0.001, 210 # }, [0.001, 90.0], [5075.12, 0.001]], 211 212 # [{'thick_core': 10.0, 213 # 'thick_layer': 15.0, 214 # 'radius': 3000.0, 215 # 'n_stacking': 1.0, 216 # 'sigma_d': 0.0, 217 # 'sld_core': 4.0, 218 # 'sld_layer': -0.4, 219 # 'solvent_sd': 5.0, 220 # 'theta': 0.0, 221 # 'phi': 0.0, 222 # 'scale': 0.01, 223 # 'background': 0.001, 224 # }, ([0.4, 0.5]), [0.00105074, 0.00121761]], 225 226 # [{'thick_core': 10.0, 227 # 'thick_layer': 15.0, 228 # 'radius': 3000.0, 229 # 'n_stacking': 1.0, 230 # 'sigma_d': 0.0, 231 # 'sld_core': 4.0, 232 # 'sld_layer': -0.4, 233 # 'solvent_sd': 5.0, 234 # 'theta': 0.0, 235 # 'phi': 0.0, 236 # 'scale': 0.01, 237 # 'background': 0.001, 238 # }, ([1.3, 1.57]), [0.0010039, 0.0010038]], 239 # ] 240 # 21Mar2016 RKH notes that unit tests all have n_stacking=1, ought to test other values 241 152 # After redefinition of spherical coordinates - 153 # tests had in old coords theta=0, phi=0; new coords theta=90, phi=0 154 # but should not matter here as so far all the tests are 1D not 2D 155 tests = [ 156 # Accuracy tests based on content in test/utest_extra_models.py. 157 # Added 2 tests with n_stacked = 5 using SasView 3.1.2 - PDB; for which alas q=0.001 values seem closer to n_stacked=1 not 5, changed assuming my 4.1 code OK, RKH 158 [{'thick_core': 10.0, 159 'thick_layer': 15.0, 160 'radius': 3000.0, 161 'n_stacking': 1.0, 162 'sigma_d': 0.0, 163 'sld_core': 4.0, 164 'sld_layer': -0.4, 165 'solvent_sd': 5.0, 166 'theta': 90.0, 167 'phi': 0.0, 168 'scale': 0.01, 169 'background': 0.001, 170 }, 0.001, 5075.12], 171 [{'thick_core': 10.0, 172 'thick_layer': 15.0, 173 'radius': 3000.0, 174 'n_stacking': 5, 175 'sigma_d': 0.0, 176 'sld_core': 4.0, 177 'sld_layer': -0.4, 178 'solvent_sd': 5.0, 179 'theta': 90.0, 180 'phi': 0.0, 181 'scale': 0.01, 182 'background': 0.001, 183 # }, 0.001, 5065.12793824], n_stacking=1 not 5 ? slight change in value here 11jan2017, check other cpu types 184 # }, 0.001, 5075.11570493], 185 }, 0.001, 25325.635693], 186 [{'thick_core': 10.0, 187 'thick_layer': 15.0, 188 'radius': 3000.0, 189 'n_stacking': 5, 190 'sigma_d': 0.0, 191 'sld_core': 4.0, 192 'sld_layer': -0.4, 193 'solvent_sd': 5.0, 194 'theta': 90.0, 195 'phi': 0.0, 196 'scale': 0.01, 197 'background': 0.001, 198 # }, 0.164, 0.0127673597265], n_stacking=1 not 5 ? slight change in value here 11jan2017, check other cpu types 199 # }, 0.164, 0.01480812968], 200 }, 0.164, 0.0598367986509], 201 202 [{'thick_core': 10.0, 203 'thick_layer': 15.0, 204 'radius': 3000.0, 205 'n_stacking': 1.0, 206 'sigma_d': 0.0, 207 'sld_core': 4.0, 208 'sld_layer': -0.4, 209 'solvent_sd': 5.0, 210 'theta': 90.0, 211 'phi': 0.0, 212 'scale': 0.01, 213 'background': 0.001, 214 # second test here was at q=90, changed it to q=5, note I(q) is then just value of flat background 215 }, [0.001, 5.0], [5075.12, 0.001]], 216 217 [{'thick_core': 10.0, 218 'thick_layer': 15.0, 219 'radius': 3000.0, 220 'n_stacking': 1.0, 221 'sigma_d': 0.0, 222 'sld_core': 4.0, 223 'sld_layer': -0.4, 224 'solvent_sd': 5.0, 225 'theta': 90.0, 226 'phi': 0.0, 227 'scale': 0.01, 228 'background': 0.001, 229 }, ([0.4, 0.5]), [0.00105074, 0.00121761]], 230 231 [{'thick_core': 10.0, 232 'thick_layer': 15.0, 233 'radius': 3000.0, 234 'n_stacking': 1.0, 235 'sigma_d': 0.0, 236 'sld_core': 4.0, 237 'sld_layer': -0.4, 238 'solvent_sd': 5.0, 239 'theta': 90.0, 240 'phi': 0.0, 241 'scale': 0.01, 242 'background': 0.001, 243 }, ([1.3, 1.57]), [0.0010039, 0.0010038]], 244 ] 245 # 11Jan2017 RKH checking unit test agai, note they are all 1D, no 2D 246 -
sasmodels/models/surface_fractal.c
rb716cc6 r925ad6e 15 15 { 16 16 // calculate P(q) 17 const double pq = square(s ph_j1c(q*radius));17 const double pq = square(sas_3j1x_x(q*radius)); 18 18 19 19 // calculate S(q) -
sasmodels/models/surface_fractal.py
r5c94f41 r925ad6e 72 72 # pylint: enable=bad-whitespace, line-too-long 73 73 74 source = ["lib/s ph_j1c.c", "lib/sas_gamma.c", "surface_fractal.c"]74 source = ["lib/sas_3j1x_x.c", "lib/sas_gamma.c", "surface_fractal.c"] 75 75 76 76 demo = dict(scale=1, background=1e-5, -
sasmodels/models/triaxial_ellipsoid.c
r3a48772 r925ad6e 37 37 const double ysq = square(Gauss76Z[j]*zm + zb); 38 38 const double t = q*sqrt(acosx2 + bsinx2*(1.0-ysq) + c2*ysq); 39 const double fq = s ph_j1c(t);39 const double fq = sas_3j1x_x(t); 40 40 inner += Gauss76Wt[j] * fq * fq ; 41 41 } … … 64 64 + radius_equat_major*radius_equat_major*cmu*cmu 65 65 + radius_polar*radius_polar*calpha*calpha); 66 const double fq = s ph_j1c(t);66 const double fq = sas_3j1x_x(t); 67 67 const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 68 68 -
sasmodels/models/triaxial_ellipsoid.py
r416f5c7 r925ad6e 104 104 ] 105 105 106 source = ["lib/s ph_j1c.c", "lib/gauss76.c", "triaxial_ellipsoid.c"]106 source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "triaxial_ellipsoid.c"] 107 107 108 108 def ER(radius_equat_minor, radius_equat_major, radius_polar): -
sasmodels/models/unified_power_Rg.py
rb3f2a24 r66ca2a6 3 3 ---------- 4 4 5 Th e Beaucage model employs the empirical multiple level unified6 Exponential/Power-law fit method developed by G. Beaucage. Four functions 7 are included so that 1, 2, 3, or 4 levels can be used. In addition a 0 level 8 has been added which simplycalculates5 This model employs the empirical multiple level unified Exponential/Power-law 6 fit method developed by Beaucage. Four functions are included so that 1, 2, 3, 7 or 4 levels can be used. In addition a 0 level has been added which simply 8 calculates 9 9 10 10 .. math:: … … 15 15 many different types of particles, including fractal clusters, random coils 16 16 (Debye equation), ellipsoidal particles, etc. 17 18 The model works best for mass fractal systems characterized by Porod exponents 19 between 5/3 and 3. It should not be used for surface fractal systems. Hammouda 20 (2010) has pointed out a deficiency in the way this model handles the 21 transitioning between the Guinier and Porod regimes and which can create 22 artefacts that appear as kinks in the fitted model function. 23 24 Also see the Guinier_Porod model. 17 25 18 26 The empirical fit function is: … … 30 38 .. math:: 31 39 32 q_i^* = \frac{q}{\operatorname{erf}^3(q R_{gi}/\sqrt{6}} 40 q_i^* = q \left[\operatorname{erf} 41 \left(\frac{q R_{gi}}{\sqrt{6}}\right) 42 \right]^{-3} 33 43 34 44 … … 56 66 57 67 G Beaucage, *J. Appl. Cryst.*, 29 (1996) 134-146 68 69 B Hammouda, *Analysis of the Beaucage model, J. Appl. Cryst.*, (2010), 43, 1474-1478 58 70 59 71 """ -
sasmodels/models/vesicle.c
r3a48772 r925ad6e 30 30 contrast = sld_solvent-sld; 31 31 vol = M_4PI_3*cube(radius); 32 f = vol * s ph_j1c(q*radius) * contrast;32 f = vol * sas_3j1x_x(q*radius) * contrast; 33 33 34 34 //now the shell. No volume normalization as this is done by the caller 35 35 contrast = sld-sld_solvent; 36 36 vol = M_4PI_3*cube(radius+thickness); 37 f += vol * s ph_j1c(q*(radius+thickness)) * contrast;37 f += vol * sas_3j1x_x(q*(radius+thickness)) * contrast; 38 38 39 39 //rescale to [cm-1]. -
sasmodels/models/vesicle.py
r3a48772 r925ad6e 94 94 ] 95 95 96 source = ["lib/s ph_j1c.c", "vesicle.c"]96 source = ["lib/sas_3j1x_x.c", "vesicle.c"] 97 97 98 98 def ER(radius, thickness): -
sasmodels/resolution2d.py
r40a87fa r7e94989 64 64 # just need q_calc and weights 65 65 self.data = data 66 self.index = index 67 68 self.qx_data = data.qx_data[ index]69 self.qy_data = data.qy_data[ index]70 self.q_data = data.q_data[ index]66 self.index = index if index is not None else slice(None) 67 68 self.qx_data = data.qx_data[self.index] 69 self.qy_data = data.qy_data[self.index] 70 self.q_data = data.q_data[self.index] 71 71 72 72 dqx = getattr(data, 'dqx_data', None) … … 74 74 if dqx is not None and dqy is not None: 75 75 # Here dqx and dqy mean dq_parr and dq_perp 76 self.dqx_data = dqx[ index]77 self.dqy_data = dqy[ index]76 self.dqx_data = dqx[self.index] 77 self.dqy_data = dqy[self.index] 78 78 ## Remove singular points if exists 79 79 self.dqx_data[self.dqx_data < SIGMA_ZERO] = SIGMA_ZERO … … 126 126 # The angle (phi) of the original q point 127 127 q_phi = np.arctan(q_phi).repeat(nbins)\ 128 .reshape( nq, nbins).transpose().flatten()128 .reshape([nq, nbins]).transpose().flatten() 129 129 ## Find Gaussian weight for each dq bins: The weight depends only 130 130 # on r-direction (The integration may not need) -
sasmodels/sasview_model.py
rbb584b3 r64614ad 583 583 % type(qdist)) 584 584 585 def get_composition_models(self): 586 """ 587 Returns usable models that compose this model 588 """ 589 s_model = None 590 p_model = None 591 if hasattr(self._model_info, "composition") \ 592 and self._model_info.composition is not None: 593 p_model = _make_model_from_info(self._model_info.composition[1][0])() 594 s_model = _make_model_from_info(self._model_info.composition[1][1])() 595 return p_model, s_model 596 585 597 def calculate_Iq(self, qx, qy=None): 586 598 # type: (Sequence[float], Optional[Sequence[float]]) -> np.ndarray
Note: See TracChangeset
for help on using the changeset viewer.