Changes in / [d77eca8:630156b] in sasmodels
- Files:
-
- 8 added
- 8 deleted
- 29 edited
Legend:
- Unmodified
- Added
- Removed
-
explore/angular_pd.py
r8267e0b r12eb36b 47 47 48 48 def draw_mesh_new(ax, theta, dtheta, phi, dphi, flow, radius=10., dist='gauss'): 49 theta_center = radians( 90-theta)49 theta_center = radians(theta) 50 50 phi_center = radians(phi) 51 51 flow_center = radians(flow) … … 137 137 radius=11., dist=dist) 138 138 if not axis.startswith('d'): 139 ax.view_init(elev= 90-theta if use_new elsetheta, azim=phi)139 ax.view_init(elev=theta, azim=phi) 140 140 plt.gcf().canvas.draw() 141 141 -
sasmodels/kernel_header.c
rb00a646 r1e7b0db0 148 148 inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; } 149 149 150 // To rotate from the canonical position to theta, phi, psi, first rotate by151 // psi about the major axis, oriented along z, which is a rotation in the152 // detector plane xy. Next rotate by theta about the y axis, aligning the major153 // axis in the xz plane. Finally, rotate by phi in the detector plane xy.154 // To compute the scattering, undo these rotations in reverse order:155 // rotate in xy by -phi, rotate in xz by -theta, rotate in xy by -psi156 // The returned q is the length of the q vector and (xhat, yhat, zhat) is a unit157 // vector in the q direction.158 // To change between counterclockwise and clockwise rotation, change the159 // sign of phi and psi.160 161 150 #if 1 162 151 //think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 … … 177 166 #endif 178 167 179 #if 1180 #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat) do { \181 q = sqrt(qx*qx + qy*qy); \182 const double qxhat = qx/q; \183 const double qyhat = qy/q; \184 double sin_theta, cos_theta; \185 double sin_phi, cos_phi; \186 double sin_psi, cos_psi; \187 SINCOS(theta*M_PI_180, sin_theta, cos_theta); \188 SINCOS(phi*M_PI_180, sin_phi, cos_phi); \189 SINCOS(psi*M_PI_180, sin_psi, cos_psi); \190 xhat = qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi) \191 + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi); \192 yhat = qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi) \193 + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi); \194 zhat = qxhat*(-sin_theta*cos_phi) \195 + qyhat*(-sin_theta*sin_phi); \196 } while (0)197 #else198 // SasView 3.x definition of orientation199 168 #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 200 169 q = sqrt(qx*qx + qy*qy); \ … … 211 180 cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 212 181 } while (0) 213 #endif -
sasmodels/models/barbell.py
r9802ab3 rfcb33e4 68 68 The 2D scattering intensity is calculated similar to the 2D cylinder model. 69 69 70 .. figure:: img/cylinder_angle_definition. png70 .. figure:: img/cylinder_angle_definition.jpg 71 71 72 72 Definition of the angles for oriented 2D barbells. … … 87 87 * **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 88 88 """ 89 from numpy import inf , sin, cos, pi89 from numpy import inf 90 90 91 91 name = "barbell" … … 108 108 ["radius", "Ang", 20, [0, inf], "volume", "Cylindrical bar radius"], 109 109 ["length", "Ang", 400, [0, inf], "volume", "Cylinder bar length"], 110 ["theta", "degrees", 60, [- 360, 360], "orientation", "Barbell axis to beamangle"],111 ["phi", "degrees", 60, [- 360, 360], "orientation", "Rotation about beam"],110 ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"], 111 ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"], 112 112 ] 113 113 # pylint: enable=bad-whitespace, line-too-long … … 125 125 phi_pd=15, phi_pd_n=0, 126 126 ) 127 q = 0.1128 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!129 qx = q*cos(pi/6.0)130 qy = q*sin(pi/6.0)131 tests = [[{}, 0.075, 25.5691260532],132 [{'theta':80., 'phi':10.}, (qx, qy), 3.04233067789],133 ]134 del qx, qy # not necessary to delete, but cleaner -
sasmodels/models/bcc_paracrystal.c
r50beefe r4962519 90 90 double theta, double phi, double psi) 91 91 { 92 double q, zhat, yhat, xhat;93 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);92 double q, cos_a1, cos_a2, cos_a3; 93 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 94 94 95 const double a1 = + xhat - zhat + yhat;96 const double a2 = + xhat + zhat - yhat;97 const double a3 = - xhat + zhat + yhat;95 const double a1 = +cos_a3 - cos_a1 + cos_a2; 96 const double a2 = +cos_a3 + cos_a1 - cos_a2; 97 const double a3 = -cos_a3 + cos_a1 + cos_a2; 98 98 99 99 const double qd = 0.5*q*dnn; -
sasmodels/models/bcc_paracrystal.py
r69e1afc r925ad6e 79 79 be accurate. 80 80 81 .. figure:: img/ parallelepiped_angle_definition.png81 .. figure:: img/bcc_angle_definition.png 82 82 83 Orientation of the crystal with respect to the scattering plane, when 84 $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis). 83 Orientation of the crystal with respect to the scattering plane. 85 84 86 85 References … … 100 99 """ 101 100 102 from numpy import inf , pi101 from numpy import inf 103 102 104 103 name = "bcc_paracrystal" … … 123 122 ["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Particle scattering length density"], 124 123 ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Solvent scattering length density"], 125 ["theta", "degrees", 60, [- 360, 360], "orientation", "c axis to beamangle"],126 ["phi", "degrees", 60, [- 360, 360], "orientation", "rotation about beam"],127 ["psi", "degrees", 60, [- 360, 360], "orientation", "rotation about c axis"]124 ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"], 125 ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"], 126 ["psi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"] 128 127 ] 129 128 # pylint: enable=bad-whitespace, line-too-long … … 142 141 psi_pd=15, psi_pd_n=0, 143 142 ) 144 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!145 # add 2d test later146 q =4.*pi/220.147 tests = [148 [{ },149 [0.001, q, 0.215268], [1.46601394721, 2.85851284174, 0.00866710287078]],150 [{'theta':20.0,'phi':30,'psi':40.0},(-0.017,0.035),2082.20264399 ],151 [{'theta':20.0,'phi':30,'psi':40.0},(-0.081,0.011),0.436323144781 ]152 ] -
sasmodels/models/capped_cylinder.py
r9802ab3 rfcb33e4 71 71 The 2D scattering intensity is calculated similar to the 2D cylinder model. 72 72 73 .. figure:: img/cylinder_angle_definition. png73 .. figure:: img/cylinder_angle_definition.jpg 74 74 75 75 Definition of the angles for oriented 2D cylinders. … … 91 91 92 92 """ 93 from numpy import inf , sin, cos, pi93 from numpy import inf 94 94 95 95 name = "capped_cylinder" … … 129 129 ["radius_cap", "Ang", 20, [0, inf], "volume", "Cap radius"], 130 130 ["length", "Ang", 400, [0, inf], "volume", "Cylinder length"], 131 ["theta", "degrees", 60, [- 360, 360], "orientation", "cylinder axis to beamangle"],132 ["phi", "degrees", 60, [- 360, 360], "orientation", "rotation about beam"],131 ["theta", "degrees", 60, [-inf, inf], "orientation", "inclination angle"], 132 ["phi", "degrees", 60, [-inf, inf], "orientation", "deflection angle"], 133 133 ] 134 134 # pylint: enable=bad-whitespace, line-too-long … … 145 145 theta_pd=15, theta_pd_n=45, 146 146 phi_pd=15, phi_pd_n=1) 147 q = 0.1148 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!149 qx = q*cos(pi/6.0)150 qy = q*sin(pi/6.0)151 tests = [[{}, 0.075, 26.0698570695],152 [{'theta':80., 'phi':10.}, (qx, qy), 0.561811990502],153 ]154 del qx, qy # not necessary to delete, but cleaner -
sasmodels/models/core_shell_bicelle.c
rb260926 r592343f 30 30 31 31 static double 32 bicelle_kernel(double q ,32 bicelle_kernel(double qq, 33 33 double rad, 34 34 double radthick, 35 35 double facthick, 36 double halflength,36 double length, 37 37 double rhoc, 38 38 double rhoh, … … 42 42 double cos_alpha) 43 43 { 44 double si1,si2,be1,be2; 45 44 46 const double dr1 = rhoc-rhoh; 45 47 const double dr2 = rhor-rhosolv; 46 48 const double dr3 = rhoh-rhor; 47 const double vol1 = M_PI*square(rad)*2.0*(halflength); 48 const double vol2 = M_PI*square(rad+radthick)*2.0*(halflength+facthick); 49 const double vol3 = M_PI*square(rad)*2.0*(halflength+facthick); 49 const double vol1 = M_PI*rad*rad*(2.0*length); 50 const double vol2 = M_PI*(rad+radthick)*(rad+radthick)*2.0*(length+facthick); 51 const double vol3 = M_PI*rad*rad*2.0*(length+facthick); 52 double besarg1 = qq*rad*sin_alpha; 53 double besarg2 = qq*(rad+radthick)*sin_alpha; 54 double sinarg1 = qq*length*cos_alpha; 55 double sinarg2 = qq*(length+facthick)*cos_alpha; 50 56 51 const double be1 = sas_2J1x_x(q*(rad)*sin_alpha);52 const double be2 = sas_2J1x_x(q*(rad+radthick)*sin_alpha);53 const double si1 = sas_sinx_x(q*(halflength)*cos_alpha);54 const double si2 = sas_sinx_x(q*(halflength+facthick)*cos_alpha);57 be1 = sas_2J1x_x(besarg1); 58 be2 = sas_2J1x_x(besarg2); 59 si1 = sas_sinx_x(sinarg1); 60 si2 = sas_sinx_x(sinarg2); 55 61 56 62 const double t = vol1*dr1*si1*be1 + … … 58 64 vol3*dr3*si2*be1; 59 65 60 const double retval = t*t ;66 const double retval = t*t*sin_alpha; 61 67 62 68 return retval; … … 65 71 66 72 static double 67 bicelle_integration(double q ,73 bicelle_integration(double qq, 68 74 double rad, 69 75 double radthick, … … 77 83 // set up the integration end points 78 84 const double uplim = M_PI_4; 79 const double half length= 0.5*length;85 const double halfheight = 0.5*length; 80 86 81 87 double summ = 0.0; … … 84 90 double sin_alpha, cos_alpha; // slots to hold sincos function output 85 91 SINCOS(alpha, sin_alpha, cos_alpha); 86 double yyy = Gauss76Wt[i] * bicelle_kernel(q , rad, radthick, facthick,87 half length, rhoc, rhoh, rhor, rhosolv,92 double yyy = Gauss76Wt[i] * bicelle_kernel(qq, rad, radthick, facthick, 93 halfheight, rhoc, rhoh, rhor, rhosolv, 88 94 sin_alpha, cos_alpha); 89 summ += yyy *sin_alpha;95 summ += yyy; 90 96 } 91 97 … … 113 119 double answer = bicelle_kernel(q, radius, thick_rim, thick_face, 114 120 0.5*length, core_sld, face_sld, rim_sld, 115 solvent_sld, sin_alpha, cos_alpha); 116 return 1.0e-4*answer; 121 solvent_sld, sin_alpha, cos_alpha) / fabs(sin_alpha); 122 123 answer *= 1.0e-4; 124 125 return answer; 117 126 } 118 127 -
sasmodels/models/core_shell_bicelle.py
r9802ab3 r8afefae 47 47 .. math:: 48 48 49 \begin{align}49 \begin{align} 50 50 F(Q,\alpha) = &\bigg[ 51 51 (\rho_c - \rho_f) V_c \frac{2J_1(QRsin \alpha)}{QRsin\alpha}\frac{sin(QLcos\alpha/2)}{Q(L/2)cos\alpha} \\ … … 63 63 cylinders is then given by integrating over all possible $\theta$ and $\phi$. 64 64 65 For oriented bicelles the *theta*, and *phi* orientation parameters will appear when fitting 2D data, 66 see the :ref:`cylinder` model for further information. 65 The *theta* and *phi* parameters are not used for the 1D output. 67 66 Our implementation of the scattering kernel and the 1D scattering intensity 68 67 use the c-library from NIST. 69 68 70 .. figure:: img/cylinder_angle_definition. png69 .. figure:: img/cylinder_angle_definition.jpg 71 70 72 Definition of the angles for the oriented core shell bicelle model, 73 note that the cylinder axis of the bicelle starts along the beam direction 74 when $\theta = \phi = 0$. 71 Definition of the angles for the oriented core shell bicelle tmodel. 75 72 76 73 … … 91 88 """ 92 89 93 from numpy import inf, sin, cos , pi90 from numpy import inf, sin, cos 94 91 95 92 name = "core_shell_bicelle" … … 138 135 ["sld_rim", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Cylinder rim scattering length density"], 139 136 ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Solvent scattering length density"], 140 ["theta", "degrees", 90, [- 360, 360], "orientation", "cylinder axis to beamangle"],141 ["phi", "degrees", 0, [- 360, 360], "orientation", "rotation about beam"]137 ["theta", "degrees", 90, [-inf, inf], "orientation", "In plane angle"], 138 ["phi", "degrees", 0, [-inf, inf], "orientation", "Out of plane angle"], 142 139 ] 143 140 … … 158 155 theta=90, 159 156 phi=0) 160 q = 0.1161 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!162 qx = q*cos(pi/6.0)163 qy = q*sin(pi/6.0)164 tests = [[{}, 0.05, 7.4883545957],165 [{'theta':80., 'phi':10.}, (qx, qy), 2.81048892474 ]166 ]167 del qx, qy # not necessary to delete, but cleaner168 157 158 #qx, qy = 0.4 * cos(pi/2.0), 0.5 * sin(0) -
sasmodels/models/core_shell_bicelle_elliptical.c
rdedcf34 r592343f 1 double form_volume(double radius, double x_core, double thick_rim, double thick_face, double length); 2 double Iq(double q, 3 double radius, 4 double x_core, 5 double thick_rim, 6 double thick_face, 7 double length, 8 double core_sld, 9 double face_sld, 10 double rim_sld, 11 double solvent_sld); 12 13 14 double Iqxy(double qx, double qy, 15 double radius, 16 double x_core, 17 double thick_rim, 18 double thick_face, 19 double length, 20 double core_sld, 21 double face_sld, 22 double rim_sld, 23 double solvent_sld, 24 double theta, 25 double phi, 26 double psi); 27 1 28 // NOTE that "length" here is the full height of the core! 2 static double 3 form_volume(double r_minor, 4 double x_core, 5 double thick_rim, 6 double thick_face, 7 double length) 29 double form_volume(double radius, double x_core, double thick_rim, double thick_face, double length) 8 30 { 9 return M_PI*(r _minor+thick_rim)*(r_minor*x_core+thick_rim)*(length+2.0*thick_face);31 return M_PI*(radius+thick_rim)*(radius*x_core+thick_rim)*(length+2.0*thick_face); 10 32 } 11 33 12 static double 13 Iq(doubleq,14 double r_minor,15 double x_core,16 double thick_rim,17 double thick_face,18 double length,19 double rhoc,20 double rhoh,21 double rhor,22 double rhosolv)34 double 35 Iq(double qq, 36 double rad, 37 double x_core, 38 double radthick, 39 double facthick, 40 double length, 41 double rhoc, 42 double rhoh, 43 double rhor, 44 double rhosolv) 23 45 { 24 46 double si1,si2,be1,be2; 25 47 // core_shell_bicelle_elliptical, RKH Dec 2016, based on elliptical_cylinder and core_shell_bicelle 26 // tested against limiting cases of cylinder, elliptical_cylinder , stacked_discs,and core_shell_bicelle48 // tested against limiting cases of cylinder, elliptical_cylinder and core_shell_bicelle 27 49 // const double uplim = M_PI_4; 28 50 const double halfheight = 0.5*length; … … 33 55 //const double vbj=M_PI; 34 56 35 const double r _major = r_minor* x_core;36 const double r 2A = 0.5*(square(r_major) + square(r_minor));37 const double r 2B = 0.5*(square(r_major) - square(r_minor));38 const double dr1 = (rhoc-rhoh) *M_PI*r _minor*r_major*(2.0*halfheight);;39 const double dr2 = (rhor-rhosolv)*M_PI*(r _minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);40 const double dr3 = (rhoh-rhor) *M_PI*r _minor*r_major*2.0*(halfheight+thick_face);41 //const double vol1 = M_PI*r _minor*r_major*(2.0*halfheight);42 //const double vol2 = M_PI*(r _minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);43 //const double vol3 = M_PI*r _minor*r_major*2.0*(halfheight+thick_face);57 const double radius_major = rad * x_core; 58 const double rA = 0.5*(square(radius_major) + square(rad)); 59 const double rB = 0.5*(square(radius_major) - square(rad)); 60 const double dr1 = (rhoc-rhoh) *M_PI*rad*radius_major*(2.0*halfheight);; 61 const double dr2 = (rhor-rhosolv)*M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick); 62 const double dr3 = (rhoh-rhor) *M_PI*rad*radius_major*2.0*(halfheight+facthick); 63 //const double vol1 = M_PI*rad*radius_major*(2.0*halfheight); 64 //const double vol2 = M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick); 65 //const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick); 44 66 45 67 //initialize integral … … 52 74 const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 53 75 double inner_sum=0; 54 double sinarg1 = q *halfheight*cos_alpha;55 double sinarg2 = q *(halfheight+thick_face)*cos_alpha;76 double sinarg1 = qq*halfheight*cos_alpha; 77 double sinarg2 = qq*(halfheight+facthick)*cos_alpha; 56 78 si1 = sas_sinx_x(sinarg1); 57 79 si2 = sas_sinx_x(sinarg2); … … 60 82 //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 61 83 const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 62 const double rr = sqrt(r 2A - r2B*cos(beta));63 double besarg1 = q *rr*sin_alpha;64 double besarg2 = q *(rr+thick_rim)*sin_alpha;84 const double rr = sqrt(rA - rB*cos(beta)); 85 double besarg1 = qq*rr*sin_alpha; 86 double besarg2 = qq*(rr+radthick)*sin_alpha; 65 87 be1 = sas_2J1x_x(besarg1); 66 88 be2 = sas_2J1x_x(besarg2); … … 76 98 } 77 99 78 static double 100 double 79 101 Iqxy(double qx, double qy, 80 double r _minor,102 double rad, 81 103 double x_core, 82 double thick_rim,83 double thick_face,104 double radthick, 105 double facthick, 84 106 double length, 85 107 double rhoc, … … 92 114 { 93 115 // THIS NEEDS TESTING 94 double q , xhat, yhat, zhat;95 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q , xhat, yhat, zhat);116 double qq, cos_val, cos_mu, cos_nu; 117 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, qq, cos_val, cos_mu, cos_nu); 96 118 const double dr1 = rhoc-rhoh; 97 119 const double dr2 = rhor-rhosolv; 98 120 const double dr3 = rhoh-rhor; 99 const double r _major = r_minor*x_core;121 const double radius_major = rad*x_core; 100 122 const double halfheight = 0.5*length; 101 const double vol1 = M_PI*r _minor*r_major*length;102 const double vol2 = M_PI*(r _minor+thick_rim)*(r_major+thick_rim)*2.0*(halfheight+thick_face);103 const double vol3 = M_PI*r _minor*r_major*2.0*(halfheight+thick_face);123 const double vol1 = M_PI*rad*radius_major*length; 124 const double vol2 = M_PI*(rad+radthick)*(radius_major+radthick)*2.0*(halfheight+facthick); 125 const double vol3 = M_PI*rad*radius_major*2.0*(halfheight+facthick); 104 126 105 // Compute effective radius in rotated coordinates106 const double r_hat = sqrt(square(r_major*xhat) + square(r_minor*yhat));107 const double rshell_hat = sqrt(square((r_major+thick_rim)*xhat)108 + square((r_minor+thick_rim)*yhat));109 const double be1 = sas_2J1x_x( q *r_hat);110 const double be2 = sas_2J1x_x( q *rshell_hat);111 const double si1 = sas_sinx_x( q *halfheight*zhat);112 const double si2 = sas_sinx_x( q *(halfheight + thick_face)*zhat);127 // Compute: r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 128 // Given: radius_major = r_ratio * radius_minor 129 // ASSUME the sin_alpha is included in the separate integration over orientation of rod angle 130 const double r = rad*sqrt(square(x_core*cos_nu) + cos_mu*cos_mu); 131 const double be1 = sas_2J1x_x( qq*r ); 132 const double be2 = sas_2J1x_x( qq*(r + radthick ) ); 133 const double si1 = sas_sinx_x( qq*halfheight*cos_val ); 134 const double si2 = sas_sinx_x( qq*(halfheight + facthick)*cos_val ); 113 135 const double Aq = square( vol1*dr1*si1*be1 + vol2*dr2*si2*be2 + vol3*dr3*si2*be1); 136 //const double vol = form_volume(radius_minor, r_ratio, length); 114 137 return 1.0e-4 * Aq; 115 138 } -
sasmodels/models/core_shell_bicelle_elliptical.py
r9802ab3 r8afefae 76 76 bicelles is then given by integrating over all possible $\alpha$ and $\psi$. 77 77 78 For oriented bicell es the *theta*, *phi* and *psi* orientation parameters willappear when fitting 2D data,78 For oriented bicellles the *theta*, *phi* and *psi* orientation parameters only appear when fitting 2D data, 79 79 see the :ref:`elliptical-cylinder` model for further information. 80 80 81 81 82 .. figure:: img/elliptical_cylinder_angle_definition. png82 .. figure:: img/elliptical_cylinder_angle_definition.jpg 83 83 84 Definition of the angles for the oriented core_shell_bicelle_elliptical particles.85 84 Definition of the angles for the oriented core_shell_bicelle_elliptical model. 85 Note that *theta* and *phi* are currently defined differently to those for the core_shell_bicelle model. 86 86 87 87 … … 99 99 """ 100 100 101 from numpy import inf, sin, cos , pi101 from numpy import inf, sin, cos 102 102 103 103 name = "core_shell_bicelle_elliptical" … … 119 119 ["radius", "Ang", 30, [0, inf], "volume", "Cylinder core radius"], 120 120 ["x_core", "None", 3, [0, inf], "volume", "axial ratio of core, X = r_polar/r_equatorial"], 121 ["thick_rim", "Ang", 122 ["thick_face", "Ang", 123 ["length", "Ang", 121 ["thick_rim", "Ang", 8, [0, inf], "volume", "Rim shell thickness"], 122 ["thick_face", "Ang", 14, [0, inf], "volume", "Cylinder face thickness"], 123 ["length", "Ang", 50, [0, inf], "volume", "Cylinder length"], 124 124 ["sld_core", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Cylinder core scattering length density"], 125 125 ["sld_face", "1e-6/Ang^2", 7, [-inf, inf], "sld", "Cylinder face scattering length density"], 126 126 ["sld_rim", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Cylinder rim scattering length density"], 127 127 ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "sld", "Solvent scattering length density"], 128 ["theta", "degrees", 90.0, [-360, 360], "orientation", "cylinder axis to beamangle"],129 ["phi", "degrees", 0, [-360, 360], "orientation", "rotation about beam"],130 ["psi", "degrees", 0, [-360, 360], "orientation", "rotation about cylinder axis"]128 ["theta", "degrees", 90, [-360, 360], "orientation", "In plane angle"], 129 ["phi", "degrees", 0, [-360, 360], "orientation", "Out of plane angle"], 130 ["psi", "degrees", 0, [-360, 360], "orientation", "Major axis angle relative to Q"], 131 131 ] 132 132 … … 150 150 psi=0) 151 151 152 q = 0.1 153 # april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct! 154 qx = q*cos(pi/6.0) 155 qy = q*sin(pi/6.0) 152 #qx, qy = 0.4 * cos(pi/2.0), 0.5 * sin(0) 156 153 157 154 tests = [ … … 162 159 'sld_core':4.0, 'sld_face':7.0, 'sld_rim':1.0, 'sld_solvent':6.0, 'background':0.0}, 163 160 0.015, 286.540286], 164 # [{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ], 165 ] 166 167 del qx, qy # not necessary to delete, but cleaner 161 ] -
sasmodels/models/core_shell_cylinder.py
r9b79f29 r8e68ea0 73 73 """ 74 74 75 from numpy import pi, inf , sin, cos75 from numpy import pi, inf 76 76 77 77 name = "core_shell_cylinder" … … 117 117 ["length", "Ang", 400, [0, inf], "volume", 118 118 "Cylinder length"], 119 ["theta", "degrees", 60, [- 360, 360], "orientation",120 " cylinder axis to beamangle"],121 ["phi", "degrees", 60, [-360, 360], "orientation",122 " rotation about beam"],119 ["theta", "degrees", 60, [-inf, inf], "orientation", 120 "In plane angle"], 121 ["phi", "degrees", 60, [-inf, inf], "orientation", 122 "Out of plane angle"], 123 123 ] 124 124 … … 151 151 theta_pd=15, theta_pd_n=45, 152 152 phi_pd=15, phi_pd_n=1) 153 q = 0.1 154 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct! 155 qx = q*cos(pi/6.0) 156 qy = q*sin(pi/6.0) 157 tests = [[{}, 0.075, 10.8552692237], 158 [{}, (qx, qy), 0.444618752741 ], 159 ] 160 del qx, qy # not necessary to delete, but cleaner 153 -
sasmodels/models/core_shell_ellipsoid.py
r9802ab3 r8e68ea0 77 77 F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha} 78 78 79 For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data,80 see the :ref:`elliptical-cylinder` model for further information.81 79 82 80 References … … 134 132 ["sld_shell", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Shell scattering length density"], 135 133 ["sld_solvent", "1e-6/Ang^2", 6.3, [-inf, inf], "sld", "Solvent scattering length density"], 136 ["theta", "degrees", 0, [- 360, 360], "orientation", "elipsoid axis to beam angle"],137 ["phi", "degrees", 0, [- 360, 360], "orientation", "rotation about beam"],134 ["theta", "degrees", 0, [-inf, inf], "orientation", "Oblate orientation wrt incoming beam"], 135 ["phi", "degrees", 0, [-inf, inf], "orientation", "Oblate orientation in the plane of the detector"], 138 136 ] 139 137 # pylint: enable=bad-whitespace, line-too-long -
sasmodels/models/core_shell_parallelepiped.c
r92dfe0c r1e7b0db0 134 134 double psi) 135 135 { 136 double q, zhat, yhat, xhat;137 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);136 double q, cos_val_a, cos_val_b, cos_val_c; 137 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a); 138 138 139 139 // cspkernel in csparallelepiped recoded here … … 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 = sas_sinx_x(0.5*q*length_a* xhat);163 double siB = sas_sinx_x(0.5*q*length_b* yhat);164 double siC = sas_sinx_x(0.5*q*length_c* zhat);165 double siAt = sas_sinx_x(0.5*q*ta* xhat);166 double siBt = sas_sinx_x(0.5*q*tb* yhat);167 double siCt = sas_sinx_x(0.5*q*tc* zhat);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_parallelepiped.py
r9b79f29 rcb0dc22 6 6 The thickness and the scattering length density of the shell or 7 7 "rim" can be different on each (pair) of faces. However at this time 8 the 1D calculationdoes **NOT** actually calculate a c face rim despite the presence of9 the parameter. Some other aspects of the 1D calculation may be wrong.8 the model does **NOT** actually calculate a c face rim despite the presence of 9 the parameter. 10 10 11 11 .. note:: 12 This model was originally ported from NIST IGOR macros. However, it is not13 yet fully understood by the SasView developers and is currently underreview.12 This model was originally ported from NIST IGOR macros. However,t is not 13 yet fully understood by the SasView developers and is currently review. 14 14 15 15 The form factor is normalized by the particle volume $V$ such that … … 48 48 amplitudes of the core and shell, in the same manner as a core-shell model. 49 49 50 .. math::51 52 F_{a}(Q,\alpha,\beta)=53 \left[\frac{\sin(\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha \sin\beta)}{\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha\sin\beta}54 - \frac{\sin(\tfrac{1}{2}QL_A\sin\alpha \sin\beta)}{\tfrac{1}{2}QL_A\sin\alpha \sin\beta} \right]55 \left[\frac{\sin(\tfrac{1}{2}QL_B\sin\alpha \sin\beta)}{\tfrac{1}{2}QL_B\sin\alpha \sin\beta} \right]56 \left[\frac{\sin(\tfrac{1}{2}QL_C\sin\alpha \sin\beta)}{\tfrac{1}{2}QL_C\sin\alpha \sin\beta} \right]57 50 58 51 .. note:: 59 52 60 Why does t_B not appear in the above equation?61 53 For the calculation of the form factor to be valid, the sides of the solid 62 MUST (perhaps not any more?)be chosen such that** $A < B < C$.54 MUST be chosen such that** $A < B < C$. 63 55 If this inequality is not satisfied, the model will not report an error, 64 56 but the calculation will not be correct and thus the result wrong. 65 57 66 58 FITTING NOTES 67 If the scale is set equal to the particle volume fraction, $\phi$, the returned59 If the scale is set equal to the particle volume fraction, |phi|, the returned 68 60 value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 69 61 However, **no interparticle interference effects are included in this … … 81 73 NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 82 74 based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 83 and length $(C+2t_C)$ values, after appropriately 84 sorting the three dimensions to give an oblate or prolate particle, to give an 85 effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied. 75 and length $(C+2t_C)$ values, and used as the effective radius 76 for $S(Q)$ when $P(Q) * S(Q)$ is applied. 86 77 87 78 To provide easy access to the orientation of the parallelepiped, we define the … … 92 83 *x*-axis of the detector. 93 84 94 .. figure:: img/parallelepiped_angle_definition. png85 .. figure:: img/parallelepiped_angle_definition.jpg 95 86 96 87 Definition of the angles for oriented core-shell parallelepipeds. 97 88 98 .. figure:: img/parallelepiped_angle_projection. png89 .. figure:: img/parallelepiped_angle_projection.jpg 99 90 100 91 Examples of the angles for oriented core-shell parallelepipeds against the … … 121 112 122 113 import numpy as np 123 from numpy import pi, inf, sqrt , cos, sin114 from numpy import pi, inf, sqrt 124 115 125 116 name = "core_shell_parallelepiped" … … 153 144 ["thick_rim_c", "Ang", 10, [0, inf], "volume", 154 145 "Thickness of C rim"], 155 ["theta", "degrees", 0, [- 360, 360], "orientation",156 " c axis to beamangle"],157 ["phi", "degrees", 0, [- 360, 360], "orientation",158 " rotation about beam"],159 ["psi", "degrees", 0, [- 360, 360], "orientation",160 " rotation about c axis"],146 ["theta", "degrees", 0, [-inf, inf], "orientation", 147 "In plane angle"], 148 ["phi", "degrees", 0, [-inf, inf], "orientation", 149 "Out of plane angle"], 150 ["psi", "degrees", 0, [-inf, inf], "orientation", 151 "Rotation angle around its own c axis against q plane"], 161 152 ] 162 153 … … 195 186 psi_pd=10, psi_pd_n=1) 196 187 197 # rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models 198 qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 188 qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 199 189 tests = [[{}, 0.2, 0.533149288477], 200 190 [{}, [0.2], [0.533149288477]], 201 [{'theta':10.0, 'phi': 20.0}, (qx, qy), 0.0853299803222],202 [{'theta':10.0, 'phi': 20.0}, [(qx, qy)], [0.0853299803222]],191 [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.032102135569], 192 [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.032102135569]], 203 193 ] 204 194 del qx, qy # not necessary to delete, but cleaner -
sasmodels/models/cylinder.py
r9802ab3 r3330bb4 38 38 39 39 40 Numerical integration is simplified by a change of variable to $u = cos(\alpha)$ with 41 $sin(\alpha)=\sqrt{1-u^2}$. 40 Numerical integration is simplified by a change of variable to $u = cos(\alpha)$ with 41 $sin(\alpha)=\sqrt{1-u^2}$. 42 42 43 43 The output of the 1D scattering intensity function for randomly oriented … … 61 61 .. _cylinder-angle-definition: 62 62 63 .. figure:: img/cylinder_angle_definition. png63 .. figure:: img/cylinder_angle_definition.jpg 64 64 65 Definition of the $\theta$ and $\phi$ orientation angles for a cylinder relative 66 to the beam line coordinates, plus an indication of their orientation distributions 67 which are described as rotations about each of the perpendicular axes $\delta_1$ and $\delta_2$ 68 in the frame of the cylinder itself, which when $\theta = \phi = 0$ are parallel to the $Y$ and $X$ axes. 65 Definition of the angles for oriented cylinders. 69 66 70 .. figure:: img/cylinder_angle_projection.png 71 72 Examples for oriented cylinders. 73 74 The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data. 75 On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will 76 appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, which when $\theta = \phi = 0$ are parallel 77 to the $Y$ and $X$ axes of the instrument respectively. Some experimentation may be required to understand the 2d patterns fully. 78 (Earlier implementations had numerical integration issues in some circumstances when orientation distributions passed through 90 degrees, such 79 situations, with very broad distributions, should still be approached with care.) 67 The $\theta$ and $\phi$ parameters only appear in the model when fitting 2d data. 80 68 81 69 Validation … … 135 123 ["length", "Ang", 400, [0, inf], "volume", 136 124 "Cylinder length"], 137 ["theta", "degrees", 60, [- 360, 360], "orientation",138 " cylinder axis to beam angle"],139 ["phi", "degrees", 60, [-360, 360], "orientation",140 " rotation about beam"],125 ["theta", "degrees", 60, [-inf, inf], "orientation", 126 "latitude"], 127 ["phi", "degrees", 60, [-inf, inf], "orientation", 128 "longitude"], 141 129 ] 142 130 … … 164 152 tests = [[{}, 0.2, 0.042761386790780453], 165 153 [{}, [0.2], [0.042761386790780453]], 166 # new coords 154 # new coords 167 155 [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852], 168 156 [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]], -
sasmodels/models/ellipsoid.c
r3b571ae r130d4c7 3 3 double Iqxy(double qx, double qy, double sld, double sld_solvent, 4 4 double radius_polar, double radius_equatorial, double theta, double phi); 5 6 static double 7 _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double cos_alpha) 8 { 9 double ratio = radius_polar/radius_equatorial; 10 // Using ratio v = Rp/Re, we can expand the following to match the 11 // form given in Guinier (1955) 12 // r = Re * sqrt(1 + cos^2(T) (v^2 - 1)) 13 // = Re * sqrt( (1 - cos^2(T)) + v^2 cos^2(T) ) 14 // = Re * sqrt( sin^2(T) + v^2 cos^2(T) ) 15 // = sqrt( Re^2 sin^2(T) + Rp^2 cos^2(T) ) 16 // 17 // Instead of using pythagoras we could pass in sin and cos; this may be 18 // slightly better for 2D which has already computed it, but it introduces 19 // an extra sqrt and square for 1-D not required by the current form, so 20 // leave it as is. 21 const double r = radius_equatorial 22 * sqrt(1.0 + cos_alpha*cos_alpha*(ratio*ratio - 1.0)); 23 const double f = sas_3j1x_x(q*r); 24 25 return f*f; 26 } 5 27 6 28 double form_volume(double radius_polar, double radius_equatorial) … … 15 37 double radius_equatorial) 16 38 { 17 // Using ratio v = Rp/Re, we can implement the form given in Guinier (1955)18 // i(h) = int_0^pi/2 Phi^2(h a sqrt(cos^2 + v^2 sin^2) cos dT19 // = int_0^pi/2 Phi^2(h a sqrt((1-sin^2) + v^2 sin^2) cos dT20 // = int_0^pi/2 Phi^2(h a sqrt(1 + sin^2(v^2-1)) cos dT21 // u-substitution of22 // u = sin, du = cos dT23 // i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du24 const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0;25 26 39 // translate a point in [-1,1] to a point in [0, 1] 27 // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2;28 40 const double zm = 0.5; 29 41 const double zb = 0.5; 30 42 double total = 0.0; 31 43 for (int i=0;i<76;i++) { 32 const double u = Gauss76Z[i]*zm + zb; 33 const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 34 const double f = sas_3j1x_x(q*r); 35 total += Gauss76Wt[i] * f * f; 44 //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 45 const double cos_alpha = Gauss76Z[i]*zm + zb; 46 total += Gauss76Wt[i] * _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 36 47 } 37 48 // translate dx in [-1,1] to dx in [lower,upper] … … 51 62 double q, sin_alpha, cos_alpha; 52 63 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 53 const double r = sqrt(square(radius_equatorial*sin_alpha) 54 + square(radius_polar*cos_alpha)); 55 const double f = sas_3j1x_x(q*r); 64 const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 56 65 const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 57 66 58 return 1.0e-4 * square(f * s);67 return 1.0e-4 * form * s * s; 59 68 } 60 69 -
sasmodels/models/ellipsoid.py
r9b79f29 r925ad6e 18 18 .. math:: 19 19 20 F(q,\alpha) = \Delta \rho V \frac{3(\sin qr - qr \cos qr)}{(qr)^3} 20 F(q,\alpha) = \frac{3 \Delta \rho V (\sin[qr(R_p,R_e,\alpha)] 21 - \cos[qr(R_p,R_e,\alpha)])} 22 {[qr(R_p,R_e,\alpha)]^3} 21 23 22 for 24 and 23 25 24 26 .. math:: 25 27 26 r = \left[ R_e^2 \sin^2 \alpha + R_p^2 \cos^2 \alpha \right]^{1/2} 28 r(R_p,R_e,\alpha) = \left[ R_e^2 \sin^2 \alpha 29 + R_p^2 \cos^2 \alpha \right]^{1/2} 27 30 28 31 29 32 $\alpha$ is the angle between the axis of the ellipsoid and $\vec q$, 30 $V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid, $R_p$ is the polar 31 radius along the rotational axis of the ellipsoid, $R_e$ is the equatorial 32 radius perpendicular to the rotational axis of the ellipsoid and 33 $\Delta \rho$ (contrast) is the scattering length density difference between 34 the scatterer and the solvent. 33 $V = (4/3)\pi R_pR_e^2$ is the volume of the ellipsoid , $R_p$ is the polar radius along the 34 rotational axis of the ellipsoid, $R_e$ is the equatorial radius perpendicular 35 to the rotational axis of the ellipsoid and $\Delta \rho$ (contrast) is the 36 scattering length density difference between the scatterer and the solvent. 35 37 36 For randomly oriented particles use the orientational average,38 For randomly oriented particles: 37 39 38 40 .. math:: 39 41 40 \langle F^2(q) \rangle = \int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)\,d\alpha}42 F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha} 41 43 42 43 computed via substitution of $u=\sin(\alpha)$, $du=\cos(\alpha)\,d\alpha$ as44 45 .. math::46 47 \langle F^2(q) \rangle = \int_0^1{F^2(q, u)\,du}48 49 with50 51 .. math::52 53 r = R_e \left[ 1 + u^2\left(R_p^2/R_e^2 - 1\right)\right]^{1/2}54 44 55 45 To provide easy access to the orientation of the ellipsoid, we define … … 58 48 :ref:`cylinder orientation figure <cylinder-angle-definition>`. 59 49 For the ellipsoid, $\theta$ is the angle between the rotational axis 60 and the $z$ -axis in the $xz$ plane followed by a rotation by $\phi$ 61 in the $xy$ plane. 50 and the $z$ -axis. 62 51 63 52 NB: The 2nd virial coefficient of the solid ellipsoid is calculated based … … 101 90 than 500. 102 91 103 Model was also tested against the triaxial ellipsoid model with equal major104 and minor equatorial radii. It is also consistent with the cyclinder model105 with polar radius equal to length and equatorial radius equal to radius.106 107 92 References 108 93 ---------- … … 111 96 *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, 112 97 Plenum Press, New York, 1987. 113 114 Authorship and Verification115 ----------------------------116 117 * **Author:** NIST IGOR/DANSE **Date:** pre 2010118 * **Converted to sasmodels by:** Helen Park **Date:** July 9, 2014119 * **Last Modified by:** Paul Kienzle **Date:** March 22, 2017120 98 """ 121 99 122 from numpy import inf , sin, cos, pi100 from numpy import inf 123 101 124 102 name = "ellipsoid" … … 151 129 ["radius_equatorial", "Ang", 400, [0, inf], "volume", 152 130 "Equatorial radius"], 153 ["theta", "degrees", 60, [- 360, 360], "orientation",154 " ellipsoid axis to beamangle"],155 ["phi", "degrees", 60, [- 360, 360], "orientation",156 " rotation about beam"],131 ["theta", "degrees", 60, [-inf, inf], "orientation", 132 "In plane angle"], 133 ["phi", "degrees", 60, [-inf, inf], "orientation", 134 "Out of plane angle"], 157 135 ] 158 136 … … 161 139 def ER(radius_polar, radius_equatorial): 162 140 import numpy as np 163 # see equation (26) in A.Isihara, J.Chem.Phys. 18(1950)1446-1449 141 164 142 ee = np.empty_like(radius_polar) 165 143 idx = radius_polar > radius_equatorial … … 190 168 theta_pd=15, theta_pd_n=45, 191 169 phi_pd=15, phi_pd_n=1) 192 q = 0.1193 # april 6 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!194 qx = q*cos(pi/6.0)195 qy = q*sin(pi/6.0)196 tests = [[{}, 0.05, 54.8525847025],197 [{'theta':80., 'phi':10.}, (qx, qy), 1.74134670026 ],198 ]199 del qx, qy # not necessary to delete, but cleaner -
sasmodels/models/elliptical_cylinder.c
r61104c8 r592343f 67 67 double theta, double phi, double psi) 68 68 { 69 double q, xhat, yhat, zhat;70 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);69 double q, cos_val, cos_mu, cos_nu; 70 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val, cos_mu, cos_nu); 71 71 72 72 // Compute: r = sqrt((radius_major*cos_nu)^2 + (radius_minor*cos_mu)^2) 73 73 // Given: radius_major = r_ratio * radius_minor 74 const double r = radius_minor*sqrt(square(r_ratio* xhat) + square(yhat));74 const double r = radius_minor*sqrt(square(r_ratio*cos_nu) + cos_mu*cos_mu); 75 75 const double be = sas_2J1x_x(q*r); 76 const double si = sas_sinx_x(q* zhat*0.5*length);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
r9802ab3 rfcb33e4 57 57 define the axis of the cylinder using two angles $\theta$, $\phi$ and $\Psi$ 58 58 (see :ref:`cylinder orientation <cylinder-angle-definition>`). The angle 59 $\Psi$ is the rotational angle around its own long_c axis. 59 $\Psi$ is the rotational angle around its own long_c axis against the $q$ plane. 60 For example, $\Psi = 0$ when the $r_\text{minor}$ axis is parallel to the 61 $x$ axis of the detector. 60 62 61 63 All angle parameters are valid and given only for 2D calculation; ie, an 62 64 oriented system. 63 65 64 .. figure:: img/elliptical_cylinder_angle_definition. png66 .. figure:: img/elliptical_cylinder_angle_definition.jpg 65 67 66 Definition of angles for oriented elliptical cylinder, where axis_ratio is drawn >1, 67 and angle $\Psi$ is now a rotation around the axis of the cylinder. 68 Definition of angles for 2D 68 69 69 .. figure:: img/ elliptical_cylinder_angle_projection.png70 .. figure:: img/cylinder_angle_projection.jpg 70 71 71 72 Examples of the angles for oriented elliptical cylinders against the 72 detector plane, with $\Psi$ = 0. 73 74 The $\theta$ and $\phi$ parameters to orient the cylinder only appear in the model when fitting 2d data. 75 On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will 76 appear. These are actually rotations about the axes $\delta_1$ and $\delta_2$ of the cylinder, the $b$ and $a$ axes of the 77 cylinder cross section. (When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.) 78 The third orientation distribution, in $\psi$, is about the $c$ axis of the particle. Some experimentation may be required to 79 understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation 80 distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.) 73 detector plane. 81 74 82 75 NB: The 2nd virial coefficient of the cylinder is calculated based on the … … 115 108 """ 116 109 117 from numpy import pi, inf, sqrt , sin, cos110 from numpy import pi, inf, sqrt 118 111 119 112 name = "elliptical_cylinder" … … 132 125 ["sld", "1e-6/Ang^2", 4.0, [-inf, inf], "sld", "Cylinder scattering length density"], 133 126 ["sld_solvent", "1e-6/Ang^2", 1.0, [-inf, inf], "sld", "Solvent scattering length density"], 134 ["theta", "degrees", 90.0, [-360, 360], "orientation", " cylinder axis to beamangle"],135 ["phi", "degrees", 0, [-360, 360], "orientation", " rotation about beam"],136 ["psi", "degrees", 0, [-360, 360], "orientation", " rotation about cylinder axis"]]127 ["theta", "degrees", 90.0, [-360, 360], "orientation", "In plane angle"], 128 ["phi", "degrees", 0, [-360, 360], "orientation", "Out of plane angle"], 129 ["psi", "degrees", 0, [-360, 360], "orientation", "Major axis angle relative to Q"]] 137 130 138 131 # pylint: enable=bad-whitespace, line-too-long … … 156 149 + (length + radius) * (length + pi * radius)) 157 150 return 0.5 * (ddd) ** (1. / 3.) 158 q = 0.1159 # april 6 2017, rkh added a 2d unit test, NOT READY YET pull #890 branch assume correct!160 qx = q*cos(pi/6.0)161 qy = q*sin(pi/6.0)162 151 163 152 tests = [ … … 169 158 'sld_solvent':1.0, 'background':0.0}, 170 159 0.001, 675.504402], 171 # [{'theta':80., 'phi':10.}, (qx, qy), 7.88866563001 ],172 160 ] -
sasmodels/models/fcc_paracrystal.c
r50beefe r4962519 90 90 double theta, double phi, double psi) 91 91 { 92 double q, zhat, yhat, xhat;93 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);92 double q, cos_a1, cos_a2, cos_a3; 93 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 94 94 95 const double a1 = yhat + xhat;96 const double a2 = xhat + zhat;97 const double a3 = yhat + zhat;95 const double a1 = cos_a2 + cos_a3; 96 const double a2 = cos_a3 + cos_a1; 97 const double a3 = cos_a2 + cos_a1; 98 98 const double qd = 0.5*q*dnn; 99 99 const double arg = 0.5*square(qd*d_factor)*(a1*a1 + a2*a2 + a3*a3); -
sasmodels/models/fcc_paracrystal.py
r69e1afc r925ad6e 76 76 2D model computation. 77 77 78 .. figure:: img/ parallelepiped_angle_definition.png78 .. figure:: img/bcc_angle_definition.png 79 79 80 Orientation of the crystal with respect to the scattering plane, when 81 $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis). 80 Orientation of the crystal with respect to the scattering plane. 82 81 83 82 References … … 91 90 """ 92 91 93 from numpy import inf , pi92 from numpy import inf 94 93 95 94 name = "fcc_paracrystal" … … 111 110 ["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Particle scattering length density"], 112 111 ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Solvent scattering length density"], 113 ["theta", "degrees", 60, [-360, 360], "orientation", "c axis to beamangle"],114 ["phi", "degrees", 60, [-360, 360], "orientation", "rotation about beam"],115 ["psi", "degrees", 60, [-360, 360], "orientation", "rotation about c axis"]112 ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"], 113 ["phi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"], 114 ["psi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"] 116 115 ] 117 116 # pylint: enable=bad-whitespace, line-too-long … … 129 128 psi_pd=15, psi_pd_n=0, 130 129 ) 131 # april 10 2017, rkh add unit tests, NOT compared with any other calc method, assume correct!132 q =4.*pi/220.133 tests = [134 [{ },135 [0.001, q, 0.215268], [0.275164706668, 5.7776842567, 0.00958167119232]],136 [{}, (-0.047,-0.007), 238.103096286],137 [{}, (0.053,0.063), 0.863609587796 ],138 ] -
sasmodels/models/hollow_cylinder.py
r9b79f29 raea2e2a 60 60 """ 61 61 62 from numpy import pi, inf , sin, cos62 from numpy import pi, inf 63 63 64 64 name = "hollow_cylinder" … … 82 82 ["sld", "1/Ang^2", 6.3, [-inf, inf], "sld", "Cylinder sld"], 83 83 ["sld_solvent", "1/Ang^2", 1, [-inf, inf], "sld", "Solvent sld"], 84 ["theta", "degrees", 90, [-360, 360], "orientation", " Cylinder axis to beamangle"],85 ["phi", "degrees", 0, [-360, 360], "orientation", " Rotation about beam"],84 ["theta", "degrees", 90, [-360, 360], "orientation", "Theta angle"], 85 ["phi", "degrees", 0, [-360, 360], "orientation", "Phi angle"], 86 86 ] 87 87 # pylint: enable=bad-whitespace, line-too-long … … 129 129 theta_pd=10, theta_pd_n=5, 130 130 ) 131 q = 0.1 132 # april 6 2017, rkh added a 2d unit test, assume correct! 133 qx = q*cos(pi/6.0) 134 qy = q*sin(pi/6.0) 131 135 132 # Parameters for unit tests 136 133 tests = [ 137 134 [{}, 0.00005, 1764.926], 138 135 [{}, 'VR', 1.8], 139 [{}, 0.001, 1756.76], 140 [{}, (qx, qy), 2.36885476192 ], 141 ] 142 del qx, qy # not necessary to delete, but cleaner 136 [{}, 0.001, 1756.76] 137 ] -
sasmodels/models/parallelepiped.c
rd605080 r1e7b0db0 67 67 double psi) 68 68 { 69 double q, xhat, yhat, zhat;70 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);69 double q, cos_val_a, cos_val_b, cos_val_c; 70 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_val_c, cos_val_b, cos_val_a); 71 71 72 const double siA = sas_sinx_x(0.5* length_a*q*xhat);73 const double siB = sas_sinx_x(0.5* length_b*q*yhat);74 const double siC = sas_sinx_x(0.5* length_c*q*zhat);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/parallelepiped.py
r34a9e4e r3330bb4 9 9 ---------- 10 10 11 This model calculates the scattering from a rectangular parallelepiped12 (\:numref:`parallelepiped-image`\).13 If you need to apply polydispersity, see also :ref:`rectangular-prism`.11 | This model calculates the scattering from a rectangular parallelepiped 12 | (\:numref:`parallelepiped-image`\). 13 | If you need to apply polydispersity, see also :ref:`rectangular-prism`. 14 14 15 15 .. _parallelepiped-image: 16 16 17 18 17 .. figure:: img/parallelepiped_geometry.jpg 19 18 … … 22 21 .. note:: 23 22 24 The three dimensions of the parallelepiped (strictly here a cuboid) may be given in 25 $any$ size order. To avoid multiple fit solutions, especially 26 with Monte-Carlo fit methods, it may be advisable to restrict their ranges. There may 27 be a number of closely similar "best fits", so some trial and error, or fixing of some 28 dimensions at expected values, may help. 23 The edge of the solid must satisfy the condition that $A < B < C$. 24 This requirement is not enforced in the model, so it is up to the 25 user to check this during the analysis. 29 26 30 27 The 1D scattering intensity $I(q)$ is calculated as: … … 70 67 \mu &= qB 71 68 69 72 70 The scattering intensity per unit volume is returned in units of |cm^-1|. 73 71 74 72 NB: The 2nd virial coefficient of the parallelepiped is calculated based on 75 the averaged effective radius, after appropriately sorting the three 76 dimensions, to give an oblate or prolate particle, $(=\sqrt{AB/\pi})$ and 73 the averaged effective radius $(=\sqrt{A B / \pi})$ and 77 74 length $(= C)$ values, and used as the effective radius for 78 75 $S(q)$ when $P(q) \cdot S(q)$ is applied. … … 105 102 .. _parallelepiped-orientation: 106 103 107 .. figure:: img/parallelepiped_angle_definition. png108 109 Definition of the angles for oriented parallelepiped , shown with $A<B<C$.110 111 .. figure:: img/parallelepiped_angle_projection. png112 113 Examples of the angles for an oriented parallelepipedagainst the104 .. figure:: img/parallelepiped_angle_definition.jpg 105 106 Definition of the angles for oriented parallelepipeds. 107 108 .. figure:: img/parallelepiped_angle_projection.jpg 109 110 Examples of the angles for oriented parallelepipeds against the 114 111 detector plane. 115 112 116 On introducing "Orientational Distribution" in the angles, "distribution of theta" and "distribution of phi" parameters will117 appear. These are actually rotations about axes $\delta_1$ and $\delta_2$ of the parallelepiped, perpendicular to the $a$ x $c$ and $b$ x $c$ faces.118 (When $\theta = \phi = 0$ these are parallel to the $Y$ and $X$ axes of the instrument.) The third orientation distribution, in $\psi$, is119 about the $c$ axis of the particle, perpendicular to the $a$ x $b$ face. Some experimentation may be required to120 understand the 2d patterns fully. (Earlier implementations had numerical integration issues in some circumstances when orientation121 distributions passed through 90 degrees, such situations, with very broad distributions, should still be approached with care.)122 123 124 113 For a given orientation of the parallelepiped, the 2D form factor is 125 114 calculated as … … 127 116 .. math:: 128 117 129 P(q_x, q_y) = \left[\frac{\sin( \tfrac{1}{2}qA\cos\alpha)}{(\tfrac{1}{2}qA\cos\alpha)}\right]^2130 \left[\frac{\sin( \tfrac{1}{2}qB\cos\beta)}{(\tfrac{1}{2}qB\cos\beta)}\right]^2131 \left[\frac{\sin( \tfrac{1}{2}qC\cos\gamma)}{(\tfrac{1}{2}qC\cos\gamma)}\right]^2118 P(q_x, q_y) = \left[\frac{\sin(qA\cos\alpha/2)}{(qA\cos\alpha/2)}\right]^2 119 \left[\frac{\sin(qB\cos\beta/2)}{(qB\cos\beta/2)}\right]^2 120 \left[\frac{\sin(qC\cos\gamma/2)}{(qC\cos\gamma/2)}\right]^2 132 121 133 122 with … … 165 154 angles. 166 155 156 This model is based on form factor calculations implemented in a c-library 157 provided by the NIST Center for Neutron Research (Kline, 2006). 167 158 168 159 References … … 172 163 173 164 R Nayuk and K Huber, *Z. Phys. Chem.*, 226 (2012) 837-854 174 175 Authorship and Verification176 ----------------------------177 178 * **Author:** This model is based on form factor calculations implemented179 in a c-library provided by the NIST Center for Neutron Research (Kline, 2006).180 * **Last Modified by:** Paul Kienzle **Date:** April 05, 2017181 * **Last Reviewed by:** Richard Heenan **Date:** April 06, 2017182 183 165 """ 184 166 185 167 import numpy as np 186 from numpy import pi, inf, sqrt , sin, cos168 from numpy import pi, inf, sqrt 187 169 188 170 name = "parallelepiped" … … 198 180 mu = q*B 199 181 V: Volume of the rectangular parallelepiped 200 alpha: angle between the long axis of the 182 alpha: angle between the long axis of the 201 183 parallelepiped and the q-vector for 1D 202 184 """ … … 214 196 ["length_c", "Ang", 400, [0, inf], "volume", 215 197 "Larger side of the parallelepiped"], 216 ["theta", "degrees", 60, [- 360, 360], "orientation",217 " c axis to beamangle"],218 ["phi", "degrees", 60, [- 360, 360], "orientation",219 " rotation about beam"],220 ["psi", "degrees", 60, [- 360, 360], "orientation",221 " rotation about c axis"],198 ["theta", "degrees", 60, [-inf, inf], "orientation", 199 "In plane angle"], 200 ["phi", "degrees", 60, [-inf, inf], "orientation", 201 "Out of plane angle"], 202 ["psi", "degrees", 60, [-inf, inf], "orientation", 203 "Rotation angle around its own c axis against q plane"], 222 204 ] 223 205 … … 226 208 def ER(length_a, length_b, length_c): 227 209 """ 228 Return effective radius (ER) for P(q)*S(q)210 Return effective radius (ER) for P(q)*S(q) 229 211 """ 230 # now that axes can be in any size order, need to sort a,b,c where a~b and c is either much smaller 231 # or much larger 232 abc = np.vstack((length_a, length_b, length_c)) 233 abc = np.sort(abc, axis=0) 234 selector = (abc[1] - abc[0]) > (abc[2] - abc[1]) 235 length = np.where(selector, abc[0], abc[2]) 212 236 213 # surface average radius (rough approximation) 237 radius = np.sqrt(np.where(~selector, abc[0]*abc[1], abc[1]*abc[2])/ pi)238 239 ddd = 0.75 * radius * (2*radius*length + (length + radius)*(length + pi*radius))214 surf_rad = sqrt(length_a * length_b / pi) 215 216 ddd = 0.75 * surf_rad * (2 * surf_rad * length_c + (length_c + surf_rad) * (length_c + pi * surf_rad)) 240 217 return 0.5 * (ddd) ** (1. / 3.) 241 218 … … 253 230 phi_pd=10, phi_pd_n=1, 254 231 psi_pd=10, psi_pd_n=10) 255 # rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models 256 qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.)232 233 qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 257 234 tests = [[{}, 0.2, 0.17758004974], 258 235 [{}, [0.2], [0.17758004974]], 259 [{'theta':10.0, 'phi': 20.0}, (qx, qy), 0.0089517140475],260 [{'theta':10.0, 'phi': 20.0}, [(qx, qy)], [0.0089517140475]],236 [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.00560296014], 237 [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.00560296014]], 261 238 ] 262 239 del qx, qy # not necessary to delete, but cleaner -
sasmodels/models/sc_paracrystal.c
r50beefe r4962519 111 111 double psi) 112 112 { 113 double q, zhat, yhat, xhat;114 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);113 double q, cos_a1, cos_a2, cos_a3; 114 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_a3, cos_a2, cos_a1); 115 115 116 116 const double qd = q*dnn; … … 118 118 const double tanh_qd = tanh(arg); 119 119 const double cosh_qd = cosh(arg); 120 const double Zq = tanh_qd/(1. - cos(qd* zhat)/cosh_qd)121 * tanh_qd/(1. - cos(qd* yhat)/cosh_qd)122 * tanh_qd/(1. - cos(qd* xhat)/cosh_qd);120 const double Zq = tanh_qd/(1. - cos(qd*cos_a1)/cosh_qd) 121 * tanh_qd/(1. - cos(qd*cos_a2)/cosh_qd) 122 * tanh_qd/(1. - cos(qd*cos_a3)/cosh_qd); 123 123 124 124 const double Fq = sphere_form(q, radius, sphere_sld, solvent_sld)*Zq; -
sasmodels/models/sc_paracrystal.py
r69e1afc r3330bb4 83 83 model computation. 84 84 85 .. figure:: img/parallelepiped_angle_definition.png 86 87 Orientation of the crystal with respect to the scattering plane, when 88 $\theta = \phi = 0$ the $c$ axis is along the beam direction (the $z$ axis). 85 .. figure:: img/sc_crystal_angle_definition.jpg 89 86 90 87 Reference … … 130 127 ["sld", "1e-6/Ang^2", 3.0, [0.0, inf], "sld", "Sphere scattering length density"], 131 128 ["sld_solvent", "1e-6/Ang^2", 6.3, [0.0, inf], "sld", "Solvent scattering length density"], 132 ["theta", "degrees", 0, [-360, 360], "orientation", "c axis to beam angle"],133 ["phi", "degrees", 0, [-360, 360], "orientation", "rotation about beam"],134 ["psi", "degrees", 0, [-360, 360], "orientation", "rotation about c axis"]129 ["theta", "degrees", 0.0, [-inf, inf], "orientation", "Orientation of the a1 axis w/respect incoming beam"], 130 ["phi", "degrees", 0.0, [-inf, inf], "orientation", "Orientation of the a2 in the plane of the detector"], 131 ["psi", "degrees", 0.0, [-inf, inf], "orientation", "Orientation of the a3 in the plane of the detector"], 135 132 ] 136 133 # pylint: enable=bad-whitespace, line-too-long … … 149 146 150 147 tests = [ 151 # Accuracy tests based on content in test/utest_extra_models.py , 2d tests added April 10, 2017148 # Accuracy tests based on content in test/utest_extra_models.py 152 149 [{}, 0.001, 10.3048], 153 150 [{}, 0.215268, 0.00814889], 154 [{}, (0.414467), 0.001313289], 155 [{'theta':10.0,'phi':20,'psi':30.0},(0.045,-0.035),18.0397138402 ], 156 [{'theta':10.0,'phi':20,'psi':30.0},(0.023,0.045),0.0177333171285 ] 151 [{}, (0.414467), 0.001313289] 157 152 ] 158 153 -
sasmodels/models/stacked_disks.py
r9802ab3 rb7e8b94 58 58 and $\sigma_d$ = the Gaussian standard deviation of the d-spacing (*sigma_d*). 59 59 Note that $D\cos(\alpha)$ is the component of $D$ parallel to $q$ and the last 60 term in the equation above is effectively a Debye-Waller factor term. 60 term in the equation above is effectively a Debye-Waller factor term. 61 61 62 62 .. note:: … … 77 77 the axis of the cylinder using two angles $\theta$ and $\varphi$. 78 78 79 .. figure:: img/cylinder_angle_definition. png79 .. figure:: img/cylinder_angle_definition.jpg 80 80 81 81 Examples of the angles against the detector plane. … … 103 103 """ 104 104 105 from numpy import inf , sin, cos, pi105 from numpy import inf 106 106 107 107 name = "stacked_disks" … … 131 131 ["sld_layer", "1e-6/Ang^2", 0.0, [-inf, inf], "sld", "Layer scattering length density"], 132 132 ["sld_solvent", "1e-6/Ang^2", 5.0, [-inf, inf], "sld", "Solvent scattering length density"], 133 ["theta", "degrees", 0, [- 360, 360], "orientation", "Orientation of the stacked disk axis w/respect incoming beam"],134 ["phi", "degrees", 0, [- 360, 360], "orientation", "Rotation about beam"],133 ["theta", "degrees", 0, [-inf, inf], "orientation", "Orientation of the stacked disk axis w/respect incoming beam"], 134 ["phi", "degrees", 0, [-inf, inf], "orientation", "Orientation of the stacked disk in the plane of the detector"], 135 135 ] 136 136 # pylint: enable=bad-whitespace, line-too-long … … 152 152 # After redefinition of spherical coordinates - 153 153 # tests had in old coords theta=0, phi=0; new coords theta=90, phi=0 154 q = 0.1 155 # april 6 2017, rkh added a 2d unit test, assume correct! 156 qx = q*cos(pi/6.0) 157 qy = q*sin(pi/6.0) 154 # but should not matter here as so far all the tests are 1D not 2D 158 155 tests = [ 159 156 # Accuracy tests based on content in test/utest_extra_models.py. 160 # Added 2 tests with n_stacked = 5 using SasView 3.1.2 - PDB; 161 # for which alas q=0.001 values seem closer to n_stacked=1 not 5, 162 # changed assuming my 4.1 code OK, RKH 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 163 158 [{'thick_core': 10.0, 164 159 'thick_layer': 15.0, … … 191 186 [{'thick_core': 10.0, 192 187 'thick_layer': 15.0, 193 'radius': 100.0,194 'n_stacking': 5,195 'sigma_d': 0.0,196 'sld_core': 4.0,197 'sld_layer': -0.4,198 'solvent_sd': 5.0,199 'theta': 90.0,200 'phi': 20.0,201 'scale': 0.01,202 'background': 0.001},203 (qx, qy), 0.0491167089952 ],204 [{'thick_core': 10.0,205 'thick_layer': 15.0,206 188 'radius': 3000.0, 207 189 'n_stacking': 5, … … 246 228 'background': 0.001, 247 229 }, ([0.4, 0.5]), [0.00105074, 0.00121761]], 248 [{'thick_core': 10.0,249 'thick_layer': 15.0,250 'radius': 3000.0,251 'n_stacking': 1.0,252 'sigma_d': 0.0,253 'sld_core': 4.0,254 'sld_layer': -0.4,255 'solvent_sd': 5.0,256 'theta': 90.0,257 'phi': 20.0,258 'scale': 0.01,259 'background': 0.001,260 }, (qx, qy), 0.0341738733124 ],261 230 262 231 [{'thick_core': 10.0, -
sasmodels/models/triaxial_ellipsoid.c
r68dd6a9 r925ad6e 20 20 double radius_polar) 21 21 { 22 const double pa = square(radius_equat_minor/radius_equat_major) - 1.0; 23 const double pc = square(radius_polar/radius_equat_major) - 1.0; 24 // translate a point in [-1,1] to a point in [0, pi/2] 25 const double zm = M_PI_4; 26 const double zb = M_PI_4; 22 double sn, cn; 23 // translate a point in [-1,1] to a point in [0, 1] 24 const double zm = 0.5; 25 const double zb = 0.5; 27 26 double outer = 0.0; 28 27 for (int i=0;i<76;i++) { 29 //const double u = Gauss76Z[i]*(upper-lower)/2 + (upper + lower)/2; 30 const double phi = Gauss76Z[i]*zm + zb; 31 const double pa_sinsq_phi = pa*square(sin(phi)); 28 //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 29 const double x = 0.5*(Gauss76Z[i] + 1.0); 30 SINCOS(M_PI_2*x, sn, cn); 31 const double acosx2 = radius_equat_minor*radius_equat_minor*cn*cn; 32 const double bsinx2 = radius_equat_major*radius_equat_major*sn*sn; 33 const double c2 = radius_polar*radius_polar; 32 34 33 35 double inner = 0.0; 34 const double um = 0.5;35 const double ub = 0.5;36 36 for (int j=0;j<76;j++) { 37 // translate a point in [-1,1] to a point in [0, 1] 38 const double usq = square(Gauss76Z[j]*um + ub); 39 const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0-usq) + 1.0 + pc*usq); 40 const double fq = sas_3j1x_x(q*r); 41 inner += Gauss76Wt[j] * fq * fq; 37 const double ysq = square(Gauss76Z[j]*zm + zb); 38 const double t = q*sqrt(acosx2 + bsinx2*(1.0-ysq) + c2*ysq); 39 const double fq = sas_3j1x_x(t); 40 inner += Gauss76Wt[j] * fq * fq ; 42 41 } 43 outer += Gauss76Wt[i] * inner; // correcting for dx later42 outer += Gauss76Wt[i] * 0.5 * inner; 44 43 } 45 // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi46 const double fqsq = outer /4.0; // = outer*um*zm*8.0/(4.0*M_PI);44 // translate dx in [-1,1] to dx in [lower,upper] 45 const double fqsq = outer*zm; 47 46 const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 48 47 return 1.0e-4 * s * s * fqsq; … … 59 58 double psi) 60 59 { 61 double q, xhat, yhat, zhat;62 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat);60 double q, calpha, cmu, cnu; 61 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, calpha, cmu, cnu); 63 62 64 const double r = sqrt(square(radius_equat_minor*xhat)65 + square(radius_equat_major*yhat)66 + square(radius_polar*zhat));67 const double fq = sas_3j1x_x( q*r);63 const double t = q*sqrt(radius_equat_minor*radius_equat_minor*cnu*cnu 64 + radius_equat_major*radius_equat_major*cmu*cmu 65 + radius_polar*radius_polar*calpha*calpha); 66 const double fq = sas_3j1x_x(t); 68 67 const double s = (sld - sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 69 68 -
sasmodels/models/triaxial_ellipsoid.py
r34a9e4e r3330bb4 2 2 # Note: model title and parameter table are inserted automatically 3 3 r""" 4 All three axes are of different lengths with $R_a \leq R_b \leq R_c$ 5 **Users should maintain this inequality for all calculations**. 6 7 .. math:: 8 9 P(q) = \text{scale} V \left< F^2(q) \right> + \text{background} 10 11 where the volume $V = 4/3 \pi R_a R_b R_c$, and the averaging 12 $\left<\ldots\right>$ is applied over all orientations for 1D. 13 14 .. figure:: img/triaxial_ellipsoid_geometry.jpg 15 16 Ellipsoid schematic. 17 4 18 Definition 5 19 ---------- 6 20 7 .. figure:: img/triaxial_ellipsoid_geometry.jpg 8 9 Ellipsoid with $R_a$ as *radius_equat_minor*, $R_b$ as *radius_equat_major* 10 and $R_c$ as *radius_polar*. 11 12 Given an ellipsoid 21 The form factor calculated is 13 22 14 23 .. math:: 15 24 16 \frac{X^2}{R_a^2} + \frac{Y^2}{R_b^2} + \frac{Z^2}{R_c^2} = 1 17 18 the scattering for randomly oriented particles is defined by the average over 19 all orientations $\Omega$ of: 20 21 .. math:: 22 23 P(q) = \text{scale}(\Delta\rho)^2\frac{V}{4 \pi}\int_\Omega\Phi^2(qr)\,d\Omega 24 + \text{background} 25 P(q) = \frac{\text{scale}}{V}\int_0^1\int_0^1 26 \Phi^2(qR_a^2\cos^2( \pi x/2) + qR_b^2\sin^2(\pi y/2)(1-y^2) + R_c^2y^2) 27 dx dy 25 28 26 29 where … … 28 31 .. math:: 29 32 30 \Phi(qr) &= 3 j_1(qr)/qr = 3 (\sin qr - qr \cos qr)/(qr)^3 \\ 31 r^2 &= R_a^2e^2 + R_b^2f^2 + R_c^2g^2 \\ 32 V &= \tfrac{4}{3} \pi R_a R_b R_c 33 \Phi(u) = 3 u^{-3} (\sin u - u \cos u) 33 34 34 The $e$, $f$ and $g$ terms are the projections of the orientation vector on $X$,35 $Y$ and $Z$ respectively. Keeping the orientation fixed at the canonical36 axes, we can integrate over the incident direction using polar angle37 $-\pi/2 \le \gamma \le \pi/2$ and equatorial angle $0 \le \phi \le 2\pi$38 (as defined in ref [1]),39 40 .. math::41 42 \langle\Phi^2\rangle = \int_0^{2\pi} \int_{-\pi/2}^{\pi/2} \Phi^2(qr)43 \cos \gamma\,d\gamma d\phi44 45 with $e = \cos\gamma \sin\phi$, $f = \cos\gamma \cos\phi$ and $g = \sin\gamma$.46 A little algebra yields47 48 .. math::49 50 r^2 = b^2(p_a \sin^2 \phi \cos^2 \gamma + 1 + p_c \sin^2 \gamma)51 52 for53 54 .. math::55 56 p_a = \frac{a^2}{b^2} - 1 \text{ and } p_c = \frac{c^2}{b^2} - 157 58 Due to symmetry, the ranges can be restricted to a single quadrant59 $0 \le \gamma \le \pi/2$ and $0 \le \phi \le \pi/2$, scaling the resulting60 integral by 8. The computation is done using the substitution $u = \sin\gamma$,61 $du = \cos\gamma\,d\gamma$, giving62 63 .. math::64 65 \langle\Phi^2\rangle &= 8 \int_0^{\pi/2} \int_0^1 \Phi^2(qr) du d\phi \\66 r^2 &= b^2(p_a \sin^2(\phi)(1 - u^2) + 1 + p_c u^2)67 68 Though for convenience we describe the three radii of the ellipsoid as equatorial69 and polar, they may be given in $any$ size order. To avoid multiple solutions, especially70 with Monte-Carlo fit methods, it may be advisable to restrict their ranges. For typical71 small angle diffraction situations there may be a number of closely similar "best fits",72 so some trial and error, or fixing of some radii at expected values, may help.73 74 35 To provide easy access to the orientation of the triaxial ellipsoid, 75 36 we define the axis of the cylinder using the angles $\theta$, $\phi$ 76 and $\psi$. These angles are defined analogously to the elliptical_cylinder below, note that 77 angle $\phi$ is now NOT the same as in the equations above. 78 79 .. figure:: img/elliptical_cylinder_angle_definition.png 80 81 Definition of angles for oriented triaxial ellipsoid, where radii are for illustration here 82 $a < b << c$ and angle $\Psi$ is a rotation around the axis of the particle. 83 84 For oriented ellipsoids the *theta*, *phi* and *psi* orientation parameters will appear when fitting 2D data, 85 see the :ref:`elliptical-cylinder` model for further information. 37 and $\psi$. These angles are defined on 38 :numref:`triaxial-ellipsoid-angles` . 39 The angle $\psi$ is the rotational angle around its own $c$ axis 40 against the $q$ plane. For example, $\psi = 0$ when the 41 $a$ axis is parallel to the $x$ axis of the detector. 86 42 87 43 .. _triaxial-ellipsoid-angles: 88 44 89 .. figure:: img/triaxial_ellipsoid_angle_projection. png45 .. figure:: img/triaxial_ellipsoid_angle_projection.jpg 90 46 91 Some examples for an oriented triaxialellipsoid.47 The angles for oriented ellipsoid. 92 48 93 49 The radius-of-gyration for this system is $R_g^2 = (R_a R_b R_c)^2/5$. 94 50 95 The contrast $\Delta\rho$is defined as SLD(ellipsoid) - SLD(solvent). In the51 The contrast is defined as SLD(ellipsoid) - SLD(solvent). In the 96 52 parameters, $R_a$ is the minor equatorial radius, $R_b$ is the major 97 53 equatorial radius, and $R_c$ is the polar radius of the ellipsoid. 98 54 99 55 NB: The 2nd virial coefficient of the triaxial solid ellipsoid is 100 calculated after sorting the three radii to give the most appropriate 101 prolate or oblate form, from the new polar radius $R_p = R_c$ and effective equatorial 102 radius, $R_e = \sqrt{R_a R_b}$, to then be used as the effective radius for 56 calculated based on the polar radius $R_p = R_c$ and equatorial 57 radius $R_e = \sqrt{R_a R_b}$, and used as the effective radius for 103 58 $S(q)$ when $P(q) \cdot S(q)$ is applied. 104 59 … … 114 69 ---------- 115 70 116 [1] Finnigan, J.A., Jacobs, D.J., 1971. 117 *Light scattering by ellipsoidal particles in solution*, 118 J. Phys. D: Appl. Phys. 4, 72-77. doi:10.1088/0022-3727/4/1/310 119 120 Authorship and Verification 121 ---------------------------- 122 123 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 124 * **Last Modified by:** Paul Kienzle (improved calculation) **Date:** April 4, 2017 125 * **Last Reviewed by:** Paul Kienzle & Richard Heenan **Date:** April 4, 2017 71 L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray 72 and Neutron Scattering*, Plenum, New York, 1987. 126 73 """ 127 74 128 from numpy import inf , sin, cos, pi75 from numpy import inf 129 76 130 77 name = "triaxial_ellipsoid" 131 78 title = "Ellipsoid of uniform scattering length density with three independent axes." 132 79 133 description = """ 134 Triaxial ellipsoid - see main documentation. 80 description = """\ 81 Note: During fitting ensure that the inequality ra<rb<rc is not 82 violated. Otherwise the calculation will 83 not be correct. 135 84 """ 136 85 category = "shape:ellipsoid" … … 142 91 "Solvent scattering length density"], 143 92 ["radius_equat_minor", "Ang", 20, [0, inf], "volume", 144 "Minor equatorial radius , Ra"],93 "Minor equatorial radius"], 145 94 ["radius_equat_major", "Ang", 400, [0, inf], "volume", 146 "Major equatorial radius , Rb"],95 "Major equatorial radius"], 147 96 ["radius_polar", "Ang", 10, [0, inf], "volume", 148 "Polar radius , Rc"],149 ["theta", "degrees", 60, [- 360, 360], "orientation",150 " polar axis to beamangle"],151 ["phi", "degrees", 60, [- 360, 360], "orientation",152 " rotation about beam"],153 ["psi", "degrees", 60, [- 360, 360], "orientation",154 " rotation about polar axis"],97 "Polar radius"], 98 ["theta", "degrees", 60, [-inf, inf], "orientation", 99 "In plane angle"], 100 ["phi", "degrees", 60, [-inf, inf], "orientation", 101 "Out of plane angle"], 102 ["psi", "degrees", 60, [-inf, inf], "orientation", 103 "Out of plane angle"], 155 104 ] 156 105 … … 159 108 def ER(radius_equat_minor, radius_equat_major, radius_polar): 160 109 """ 161 Returns the effective radius used in the S*P calculation110 Returns the effective radius used in the S*P calculation 162 111 """ 163 112 import numpy as np 164 113 from .ellipsoid import ER as ellipsoid_ER 165 166 # now that radii can be in any size order, radii need sorting a,b,c 167 # where a~b and c is either much smaller or much larger 168 radii = np.vstack((radius_equat_major, radius_equat_minor, radius_polar)) 169 radii = np.sort(radii, axis=0) 170 selector = (radii[1] - radii[0]) > (radii[2] - radii[1]) 171 polar = np.where(selector, radii[0], radii[2]) 172 equatorial = np.sqrt(np.where(~selector, radii[0]*radii[1], radii[1]*radii[2])) 173 return ellipsoid_ER(polar, equatorial) 114 return ellipsoid_ER(radius_polar, np.sqrt(radius_equat_minor * radius_equat_major)) 174 115 175 116 demo = dict(scale=1, background=0, … … 183 124 phi_pd=15, phi_pd_n=1, 184 125 psi_pd=15, psi_pd_n=1) 185 186 q = 0.1187 # april 6 2017, rkh add unit tests188 # NOT compared with any other calc method, assume correct!189 # check 2d test after pull #890190 qx = q*cos(pi/6.0)191 qy = q*sin(pi/6.0)192 tests = [[{}, 0.05, 24.8839548033],193 [{'theta':80., 'phi':10.}, (qx, qy), 166.712060266 ],194 ]195 del qx, qy # not necessary to delete, but cleaner
Note: See TracChangeset
for help on using the changeset viewer.