Changeset ab2aea8 in sasmodels
- Timestamp:
- Oct 15, 2016 5:49:53 PM (8 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- e30d645
- Parents:
- ed0827a
- Location:
- sasmodels/models
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/hollow_rectangular_prism.c
r3a48772 rab2aea8 5 5 double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness) 6 6 { 7 double b_side= length_a * b2a_ratio;8 double c_side= length_a * c2a_ratio;7 double length_b = length_a * b2a_ratio; 8 double length_c = length_a * c2a_ratio; 9 9 double a_core = length_a - 2.0*thickness; 10 double b_core = b_side- 2.0*thickness;11 double c_core = c_side- 2.0*thickness;10 double b_core = length_b - 2.0*thickness; 11 double c_core = length_c - 2.0*thickness; 12 12 double vol_core = a_core * b_core * c_core; 13 double vol_total = length_a * b_side * c_side;13 double vol_total = length_a * length_b * length_c; 14 14 double vol_shell = vol_total - vol_core; 15 15 return vol_shell; … … 26 26 double termA1, termA2, termB1, termB2, termC1, termC2; 27 27 28 double b_side= length_a * b2a_ratio;29 double c_side= length_a * c2a_ratio;28 double length_b = length_a * b2a_ratio; 29 double length_c = length_a * c2a_ratio; 30 30 double a_half = 0.5 * length_a; 31 double b_half = 0.5 * b_side; 32 double c_half = 0.5 * c_side; 31 double b_half = 0.5 * length_b; 32 double c_half = 0.5 * length_c; 33 double vol_total = length_a * length_b * length_c; 34 double vol_core = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness); 33 35 34 //Integration limits to use in Gaussian quadrature36 //Integration limits to use in Gaussian quadrature 35 37 double v1a = 0.0; 36 38 double v1b = M_PI_2; //theta integration limits … … 38 40 double v2b = M_PI_2; //phi integration limits 39 41 40 //Order of integration41 int nordi=76;42 int nordj=76;42 double outer_sum = 0.0; 43 44 for(int i=0; i<76; i++) { 43 45 44 double sumi = 0.0; 45 46 for(int i=0; i<nordi; i++) { 46 double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 47 47 48 double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 48 double termC1 = sinc(q * c_half * cos(theta)); 49 double termC2 = sinc(q * (c_half-thickness)*cos(theta)); 49 50 50 double arg = q * c_half * cos(theta); 51 if (fabs(arg) > 1.e-16) {termC1 = sin(arg)/arg;} else {termC1 = 1.0;} 52 arg = q * (c_half-thickness)*cos(theta); 53 if (fabs(arg) > 1.e-16) {termC2 = sin(arg)/arg;} else {termC2 = 1.0;} 54 55 double sumj = 0.0; 51 double inner_sum = 0.0; 56 52 57 for(int j=0; j<nordj; j++) {53 for(int j=0; j<76; j++) { 58 54 59 55 double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); … … 61 57 // Amplitude AP from eqn. (13), rewritten to avoid round-off effects when arg=0 62 58 63 arg = q * a_half * sin(theta) * sin(phi); 64 if (fabs(arg) > 1.e-16) {termA1 = sin(arg)/arg;} else {termA1 = 1.0;} 65 arg = q * (a_half-thickness) * sin(theta) * sin(phi); 66 if (fabs(arg) > 1.e-16) {termA2 = sin(arg)/arg;} else {termA2 = 1.0;} 59 termA1 = sinc(q * a_half * sin(theta) * sin(phi)); 60 termA2 = sinc(q * (a_half-thickness) * sin(theta) * sin(phi)); 67 61 68 arg = q * b_half * sin(theta) * cos(phi); 69 if (fabs(arg) > 1.e-16) {termB1 = sin(arg)/arg;} else {termB1 = 1.0;} 70 arg = q * (b_half-thickness) * sin(theta) * cos(phi); 71 if (fabs(arg) > 1.e-16) {termB2 = sin(arg)/arg;} else {termB2 = 1.0;} 62 termB1 = sinc(q * b_half * sin(theta) * cos(phi)); 63 termB2 = sinc(q * (b_half-thickness) * sin(theta) * cos(phi)); 72 64 73 double AP1 = (length_a*b_side*c_side) * termA1 * termB1 * termC1; 74 double AP2 = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness) * termA2 * termB2 * termC2; 75 double AP = AP1 - AP2; 65 double AP1 = vol_total * termA1 * termB1 * termC1; 66 double AP2 = vol_core * termA2 * termB2 * termC2; 76 67 77 sumj += Gauss76Wt[j] * (AP*AP);68 inner_sum += Gauss76Wt[j] * square(AP1-AP2); 78 69 79 70 } 80 71 81 sumj = 0.5 * (v2b-v2a) * sumj;82 sumi += Gauss76Wt[i] * sumj* sin(theta);72 inner_sum = 0.5 * (v2b-v2a) * inner_sum; 73 outer_sum += Gauss76Wt[i] * inner_sum * sin(theta); 83 74 84 75 } 85 76 86 double answer = 0.5*(v1b-v1a)* sumi;77 double answer = 0.5*(v1b-v1a)*outer_sum; 87 78 88 79 // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization) … … 91 82 92 83 // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 93 answer *= (sld-solvent_sld)*(sld-solvent_sld);84 answer *= square(sld-solvent_sld); 94 85 95 86 // Convert from [1e-12 A-1] to [cm-1] -
sasmodels/models/hollow_rectangular_prism.py
ra807206 rab2aea8 13 13 the difference of the amplitudes of two massive parallelepipeds 14 14 differing in their outermost dimensions in each direction by the 15 same length increment :math:`2\Delta`(Nayuk, 2012).15 same length increment $2\Delta$ (Nayuk, 2012). 16 16 17 17 As in the case of the massive parallelepiped model (:ref:`rectangular-prism`), … … 58 58 59 59 .. math:: 60 I(q) = \text{scale} \times V \times (\rho_ {\text{p}} -61 \rho_ {\text{solvent}})^2 \times P(q) + \text{background}60 I(q) = \text{scale} \times V \times (\rho_\text{p} - 61 \rho_\text{solvent})^2 \times P(q) + \text{background} 62 62 63 where $\rho_ {\text{p}}$ is the scattering length of the parallelepiped,64 $\rho_ {\text{solvent}}$ is the scattering length of the solvent,63 where $\rho_\text{p}$ is the scattering length of the parallelepiped, 64 $\rho_\text{solvent}$ is the scattering length of the solvent, 65 65 and (if the data are in absolute units) *scale* represents the volume fraction 66 66 (which is unitless). … … 148 148 # parameters for demo 149 149 demo = dict(scale=1, background=0, 150 sld=6.3 e-6, sld_solvent=1.0e-6,150 sld=6.3, sld_solvent=1.0, 151 151 length_a=35, b2a_ratio=1, c2a_ratio=1, thickness=1, 152 152 length_a_pd=0.1, length_a_pd_n=10, -
sasmodels/models/hollow_rectangular_prism_thin_walls.c
r3a48772 rab2aea8 5 5 double form_volume(double length_a, double b2a_ratio, double c2a_ratio) 6 6 { 7 double b_side= length_a * b2a_ratio;8 double c_side= length_a * c2a_ratio;9 double vol_shell = 2.0 * (length_a* b_side + length_a*c_side + b_side*c_side);7 double length_b = length_a * b2a_ratio; 8 double length_c = length_a * c2a_ratio; 9 double vol_shell = 2.0 * (length_a*length_b + length_a*length_c + length_b*length_c); 10 10 return vol_shell; 11 11 } … … 18 18 double c2a_ratio) 19 19 { 20 double b_side= length_a * b2a_ratio;21 double c_side= length_a * c2a_ratio;22 double a_half = 0.5 * length_a;23 double b_half = 0.5 * b_side;24 double c_half = 0.5 * c_side;20 const double length_b = length_a * b2a_ratio; 21 const double length_c = length_a * c2a_ratio; 22 const double a_half = 0.5 * length_a; 23 const double b_half = 0.5 * length_b; 24 const double c_half = 0.5 * length_c; 25 25 26 26 //Integration limits to use in Gaussian quadrature 27 double v1a = 0.0;28 double v1b = M_PI_2; //theta integration limits29 double v2a = 0.0;30 double v2b = M_PI_2; //phi integration limits27 const double v1a = 0.0; 28 const double v1b = M_PI_2; //theta integration limits 29 const double v2a = 0.0; 30 const double v2b = M_PI_2; //phi integration limits 31 31 32 //Order of integration33 int nordi=76;34 int nordj=76;32 double outer_sum = 0.0; 33 for(int i=0; i<76; i++) { 34 const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 35 35 36 double sumi = 0.0; 37 38 for(int i=0; i<nordi; i++) { 36 double sin_theta, cos_theta; 37 double sin_c, cos_c; 38 SINCOS(theta, sin_theta, cos_theta); 39 SINCOS(q*c_half*cos_theta, sin_c, cos_c); 39 40 40 double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b );41 42 41 // To check potential problems if denominator goes to zero here !!! 43 double termAL_theta = 8.0*cos(q*c_half*cos(theta)) / (q*q*sin(theta)*sin(theta));44 double termAT_theta = 8.0*sin(q*c_half*cos(theta)) / (q*q*sin(theta)*cos(theta));42 const double termAL_theta = 8.0 * cos_c / (q*q*sin_theta*sin_theta); 43 const double termAT_theta = 8.0 * sin_c / (q*q*sin_theta*cos_theta); 45 44 46 double sumj= 0.0;47 48 for(int j=0; j<nordj; j++) { 45 double inner_sum = 0.0; 46 for(int j=0; j<76; j++) { 47 const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 49 48 50 double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 51 49 double sin_phi, cos_phi; 50 double sin_a, cos_a; 51 double sin_b, cos_b; 52 SINCOS(phi, sin_phi, cos_phi); 53 SINCOS(q*a_half*sin_theta*sin_phi, sin_a, cos_a); 54 SINCOS(q*b_half*sin_theta*cos_phi, sin_b, cos_b); 55 52 56 // Amplitude AL from eqn. (7c) 53 double AL = termAL_theta * sin(q*a_half*sin(theta)*sin(phi)) *54 sin(q*b_half*sin(theta)*cos(phi)) / (sin(phi)*cos(phi));57 const double AL = termAL_theta 58 * sin_a*sin_b / (sin_phi*cos_phi); 55 59 56 60 // Amplitude AT from eqn. (9) 57 double AT = termAT_theta * ( (cos(q*a_half*sin(theta)*sin(phi))*sin(q*b_half*sin(theta)*cos(phi))/cos(phi))58 + (cos(q*b_half*sin(theta)*cos(phi))*sin(q*a_half*sin(theta)*sin(phi))/sin(phi)));61 const double AT = termAT_theta 62 * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi ); 59 63 60 sumj += Gauss76Wt[j] * (AL+AT)*(AL+AT); 64 inner_sum += Gauss76Wt[j] * square(AL+AT); 65 } 61 66 62 } 63 64 sumj = 0.5 * (v2b-v2a) * sumj; 65 sumi += Gauss76Wt[i] * sumj * sin(theta); 66 67 inner_sum *= 0.5 * (v2b-v2a); 68 outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 67 69 } 68 70 69 double answer = 0.5*(v1b-v1a)*sumi;71 outer_sum *= 0.5*(v1b-v1a); 70 72 71 73 // Normalize as in Eqn. (15) without the volume factor (as cancels with (V*DelRho)^2 normalization) 72 74 // The factor 2 is due to the different theta integration limit (pi/2 instead of pi) 73 answer /=M_PI_2;75 double answer = outer_sum/M_PI_2; 74 76 75 77 // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 76 answer *= (sld-solvent_sld)*(sld-solvent_sld);78 answer *= square(sld-solvent_sld); 77 79 78 80 // Convert from [1e-12 A-1] to [cm-1] … … 80 82 81 83 return answer; 82 83 84 } -
sasmodels/models/hollow_rectangular_prism_thin_walls.py
ra807206 rab2aea8 3 3 r""" 4 4 5 This model provides the form factor, *P(q)*, for a hollow rectangular5 This model provides the form factor, $P(q)$, for a hollow rectangular 6 6 prism with infinitely thin walls. It computes only the 1D scattering, not the 2D. 7 7 … … 14 14 15 15 Assuming a hollow parallelepiped with infinitely thin walls, edge lengths 16 :math:`A \le B \le C`and presenting an orientation with respect to the17 scattering vector given by |theta| and |phi|, where |theta|is the angle18 between the *z* axis and the longest axis of the parallelepiped *C*, and19 |phi| is the angle between the scattering vector (lying in the *xy*plane)20 and the *y*axis, the form factor is given by16 $A \le B \le C$ and presenting an orientation with respect to the 17 scattering vector given by $\theta$ and $\phi$, where $\theta$ is the angle 18 between the $z$ axis and the longest axis of the parallelepiped $C$, and 19 $\phi$ is the angle between the scattering vector (lying in the $xy$ plane) 20 and the $y$ axis, the form factor is given by 21 21 22 22 .. math:: 23 P(q) = \frac{1}{V^2} \frac{2}{\pi} \int_0^{\frac{\pi}{2}} 24 \int_0^{\frac{\pi}{2}} [A_L(q)+A_T(q)]^2 \sin\theta d\theta d\phi 23 24 P(q) = \frac{1}{V^2} \frac{2}{\pi} \int_0^{\frac{\pi}{2}} 25 \int_0^{\frac{\pi}{2}} [A_L(q)+A_T(q)]^2 \sin\theta\,d\theta\,d\phi 25 26 26 27 where 27 28 28 29 .. math:: 29 V = 2AB + 2AC + 2BC30 30 31 .. math:: 32 A_L(q) = 8 \times \frac{ \sin \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr)33 \sin \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr)34 \cos \bigl( q \frac{C}{2} \cos\theta \bigr) }35 {q^2 \, \sin^2\theta \, \sin\phi \cos\phi}36 37 .. math:: 38 A_T(q) = A_F(q) \times \frac{2 \, \sin \bigl( q \frac{C}{2} \cos\theta \bigr)}{q \,\cos\theta}31 V &= 2AB + 2AC + 2BC \\ 32 A_L(q) &= 8 \times \frac{ 33 \sin \left( \tfrac{1}{2} q A \sin\phi \sin\theta \right) 34 \sin \left( \tfrac{1}{2} q B \cos\phi \sin\theta \right) 35 \cos \left( \tfrac{1}{2} q C \cos\theta \right) 36 }{q^2 \, \sin^2\theta \, \sin\phi \cos\phi} \\ 37 A_T(q) &= A_F(q) \times 38 \frac{2\,\sin \left( \tfrac{1}{2} q C \cos\theta \right)}{q\,\cos\theta} 39 39 40 40 and 41 41 42 42 .. math:: 43 A_F(q) = 4 \frac{ \cos \bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr) 44 \sin \bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr) } 43 44 A_F(q) = 4 \frac{ \cos \left( \tfrac{1}{2} q A \sin\phi \sin\theta \right) 45 \sin \left( \tfrac{1}{2} q B \cos\phi \sin\theta \right) } 45 46 {q \, \cos\phi \, \sin\theta} + 46 4 \frac{ \sin \ bigl( q \frac{A}{2} \sin\phi \sin\theta \bigr)47 \cos \ bigl( q \frac{B}{2} \cos\phi \sin\theta \bigr) }47 4 \frac{ \sin \left( \tfrac{1}{2} q A \sin\phi \sin\theta \right) 48 \cos \left( \tfrac{1}{2} q B \cos\phi \sin\theta \right) } 48 49 {q \, \sin\phi \, \sin\theta} 49 50 … … 51 52 52 53 .. math:: 53 I(q) = \mbox{scale} \times V \times (\rho_{\mbox{p}} - \rho_{\mbox{solvent}})^2 \times P(q)54 54 55 where *V* is the volume of the rectangular prism, :math:`\rho_{\mbox{p}}` 56 is the scattering length of the parallelepiped, :math:`\rho_{\mbox{solvent}}` 55 I(q) = \text{scale} \times V \times (\rho_\text{p} - \rho_\text{solvent})^2 \times P(q) 56 57 where $V$ is the volume of the rectangular prism, $\rho_\text{p}$ 58 is the scattering length of the parallelepiped, $\rho_\text{solvent}$ 57 59 is the scattering length of the solvent, and (if the data are in absolute 58 60 units) *scale* represents the volume fraction (which is unitless). … … 127 129 # parameters for demo 128 130 demo = dict(scale=1, background=0, 129 sld=6.3 e-6, sld_solvent=1.0e-6,131 sld=6.3, sld_solvent=1.0, 130 132 length_a=35, b2a_ratio=1, c2a_ratio=1, 131 133 length_a_pd=0.1, length_a_pd_n=10, -
sasmodels/models/rectangular_prism.c
r3a48772 rab2aea8 15 15 double c2a_ratio) 16 16 { 17 double termA, termB, termC; 18 19 double b_side = length_a * b2a_ratio; 20 double c_side = length_a * c2a_ratio; 21 double volume = length_a * b_side * c_side; 22 double a_half = 0.5 * length_a; 23 double b_half = 0.5 * b_side; 24 double c_half = 0.5 * c_side; 17 const double length_b = length_a * b2a_ratio; 18 const double length_c = length_a * c2a_ratio; 19 const double a_half = 0.5 * length_a; 20 const double b_half = 0.5 * length_b; 21 const double c_half = 0.5 * length_c; 25 22 26 23 //Integration limits to use in Gaussian quadrature 27 double v1a = 0.0;28 double v1b = M_PI_2; //theta integration limits29 double v2a = 0.0;30 double v2b = M_PI_2; //phi integration limits24 const double v1a = 0.0; 25 const double v1b = M_PI_2; //theta integration limits 26 const double v2a = 0.0; 27 const double v2b = M_PI_2; //phi integration limits 31 28 32 //Order of integration 33 int nordi=76; 34 int nordj=76; 29 double outer_sum = 0.0; 30 for(int i=0; i<76; i++) { 31 const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 32 double sin_theta, cos_theta; 33 SINCOS(theta, sin_theta, cos_theta); 35 34 36 double sumi = 0.0; 37 38 for(int i=0; i<nordi; i++) { 35 const double termC = sinc(q * c_half * cos_theta); 39 36 40 double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 37 double inner_sum = 0.0; 38 for(int j=0; j<76; j++) { 39 double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 40 double sin_phi, cos_phi; 41 SINCOS(phi, sin_phi, cos_phi); 41 42 42 double arg = q * c_half * cos(theta); 43 if (fabs(arg) > 1.e-16) {termC = sin(arg)/arg;} else {termC = 1.0;} 44 45 double sumj = 0.0; 46 47 for(int j=0; j<nordj; j++) { 48 49 double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 50 51 // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0 52 53 arg = q * a_half * sin(theta) * sin(phi); 54 if (fabs(arg) > 1.e-16) {termA = sin(arg)/arg;} else {termA = 1.0;} 55 56 arg = q * b_half * sin(theta) * cos(phi); 57 if (fabs(arg) > 1.e-16) {termB = sin(arg)/arg;} else {termB = 1.0;} 58 59 double AP = termA * termB * termC; 60 61 sumj += Gauss76Wt[j] * (AP*AP); 62 63 } 64 65 sumj = 0.5 * (v2b-v2a) * sumj; 66 sumi += Gauss76Wt[i] * sumj * sin(theta); 67 43 // Amplitude AP from eqn. (12), rewritten to avoid round-off effects when arg=0 44 const double termA = sinc(q * a_half * sin_theta * sin_phi); 45 const double termB = sinc(q * b_half * sin_theta * cos_phi); 46 const double AP = termA * termB * termC; 47 inner_sum += Gauss76Wt[j] * AP * AP; 48 } 49 inner_sum = 0.5 * (v2b-v2a) * inner_sum; 50 outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 68 51 } 69 52 70 double answer = 0.5*(v1b-v1a)* sumi;53 double answer = 0.5*(v1b-v1a)*outer_sum; 71 54 72 55 // Normalize by Pi (Eqn. 16). … … 78 61 79 62 // Multiply by contrast^2 and volume^2 80 answer *= (sld-solvent_sld)*(sld-solvent_sld)*volume*volume; 63 const double volume = length_a * length_b * length_c; 64 answer *= square((sld-solvent_sld)*volume); 81 65 82 66 // Convert from [1e-12 A-1] to [cm-1] … … 84 68 85 69 return answer; 86 87 70 } -
sasmodels/models/rectangular_prism.py
r2941abf rab2aea8 128 128 # parameters for demo 129 129 demo = dict(scale=1, background=0, 130 sld=6.3 e-6, sld_solvent=1.0e-6,130 sld=6.3, sld_solvent=1.0, 131 131 length_a=35, b2a_ratio=1, c2a_ratio=1, 132 132 length_a_pd=0.1, length_a_pd_n=10,
Note: See TracChangeset
for help on using the changeset viewer.