Changeset 2222134 in sasmodels for sasmodels/models
- Timestamp:
- Sep 30, 2016 11:07:16 AM (8 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- a807206
- Parents:
- 6e5b2a7
- Location:
- sasmodels/models
- Files:
-
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/barbell.c
r2f5c6d4 r2222134 1 double form_volume(double bell_radius, double radius, double length);1 double form_volume(double radius_bell, double radius, double length); 2 2 double Iq(double q, double sld, double solvent_sld, 3 double bell_radius, double radius, double length);3 double radius_bell, double radius, double length); 4 4 double Iqxy(double qx, double qy, double sld, double solvent_sld, 5 double bell_radius, double radius, double length,5 double radius_bell, double radius, double length, 6 6 double theta, double phi); 7 7 8 #define INVALID(v) (v. bell_radius< v.radius)8 #define INVALID(v) (v.radius_bell < v.radius) 9 9 10 10 //barbell kernel - same as dumbell 11 11 static double 12 _bell_kernel(double q, double h, double bell_radius,12 _bell_kernel(double q, double h, double radius_bell, 13 13 double half_length, double sin_alpha, double cos_alpha) 14 14 { 15 15 // translate a point in [-1,1] to a point in [lower,upper] 16 16 const double upper = 1.0; 17 const double lower = h/ bell_radius;17 const double lower = h/radius_bell; 18 18 const double zm = 0.5*(upper-lower); 19 19 const double zb = 0.5*(upper+lower); … … 26 26 // m = q R cos(alpha) 27 27 // b = q(L/2-h) cos(alpha) 28 const double m = q* bell_radius*cos_alpha; // cos argument slope28 const double m = q*radius_bell*cos_alpha; // cos argument slope 29 29 const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 30 const double qrst = q* bell_radius*sin_alpha; // Q*R*sin(theta)30 const double qrst = q*radius_bell*sin_alpha; // Q*R*sin(theta) 31 31 double total = 0.0; 32 32 for (int i = 0; i < 76; i++){ … … 39 39 // translate dx in [-1,1] to dx in [lower,upper] 40 40 const double integral = total*zm; 41 const double bell_Fq = 2*M_PI*cube( bell_radius)*integral;41 const double bell_Fq = 2*M_PI*cube(radius_bell)*integral; 42 42 return bell_Fq; 43 43 } 44 44 45 double form_volume(double bell_radius,45 double form_volume(double radius_bell, 46 46 double radius, 47 47 double length) … … 49 49 50 50 // bell radius should never be less than radius when this is called 51 const double hdist = sqrt(square( bell_radius) - square(radius));52 const double p1 = 2.0/3.0*cube( bell_radius);53 const double p2 = square( bell_radius)*hdist;51 const double hdist = sqrt(square(radius_bell) - square(radius)); 52 const double p1 = 2.0/3.0*cube(radius_bell); 53 const double p2 = square(radius_bell)*hdist; 54 54 const double p3 = cube(hdist)/3.0; 55 55 … … 58 58 59 59 double Iq(double q, double sld, double solvent_sld, 60 double bell_radius, double radius, double length)60 double radius_bell, double radius, double length) 61 61 { 62 const double h = -sqrt( bell_radius*bell_radius- radius*radius);62 const double h = -sqrt(radius_bell*radius_bell - radius*radius); 63 63 const double half_length = 0.5*length; 64 64 … … 72 72 SINCOS(alpha, sin_alpha, cos_alpha); 73 73 74 const double bell_Fq = _bell_kernel(q, h, bell_radius, half_length, sin_alpha, cos_alpha);74 const double bell_Fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 75 75 const double bj = sas_J1c(q*radius*sin_alpha); 76 76 const double si = sinc(q*half_length*cos_alpha); … … 90 90 double Iqxy(double qx, double qy, 91 91 double sld, double solvent_sld, 92 double bell_radius, double radius, double length,92 double radius_bell, double radius, double length, 93 93 double theta, double phi) 94 94 { … … 100 100 const double alpha = acos(cos_val); // rod angle relative to q 101 101 102 const double h = -sqrt(square( bell_radius) - square(radius));102 const double h = -sqrt(square(radius_bell) - square(radius)); 103 103 const double half_length = 0.5*length; 104 104 105 105 double sin_alpha, cos_alpha; // slots to hold sincos function output 106 106 SINCOS(alpha, sin_alpha, cos_alpha); 107 const double bell_Fq = _bell_kernel(q, h, bell_radius, half_length, sin_alpha, cos_alpha);107 const double bell_Fq = _bell_kernel(q, h, radius_bell, half_length, sin_alpha, cos_alpha); 108 108 const double bj = sas_J1c(q*radius*sin_alpha); 109 109 const double si = sinc(q*half_length*cos_alpha); -
sasmodels/models/barbell.py
r42356c8 r2222134 13 13 .. figure:: img/barbell_geometry.jpg 14 14 15 Barbell geometry, where $r$ is *radius*, $R$ is * bell_radius* and15 Barbell geometry, where $r$ is *radius*, $R$ is *radius_bell* and 16 16 $L$ is *length*. Since the end cap radius $R \geq r$ and by definition 17 17 for this geometry $h < 0$, $h$ is then defined by $r$ and $R$ as … … 103 103 parameters = [["sld", "1e-6/Ang^2", 4, [-inf, inf], "sld", "Barbell scattering length density"], 104 104 ["sld_solvent", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Solvent scattering length density"], 105 [" bell_radius", "Ang", 40, [0, inf], "volume", "Spherical bell radius"],105 ["radius_bell", "Ang", 40, [0, inf], "volume", "Spherical bell radius"], 106 106 ["radius", "Ang", 20, [0, inf], "volume", "Cylindrical bar radius"], 107 107 ["length", "Ang", 400, [0, inf], "volume", "Cylinder bar length"], … … 116 116 demo = dict(scale=1, background=0, 117 117 sld=6, sld_solvent=1, 118 bell_radius=40, radius=20, length=400,118 radius_bell=40, radius=20, length=400, 119 119 theta=60, phi=60, 120 120 radius_pd=.2, radius_pd_n=5, -
sasmodels/models/capped_cylinder.c
r2f5c6d4 r2222134 1 double form_volume(double radius, double cap_radius, double length);1 double form_volume(double radius, double radius_cap, double length); 2 2 double Iq(double q, double sld, double solvent_sld, 3 double radius, double cap_radius, double length);3 double radius, double radius_cap, double length); 4 4 double Iqxy(double qx, double qy, double sld, double solvent_sld, 5 double radius, double cap_radius, double length, double theta, double phi);5 double radius, double radius_cap, double length, double theta, double phi); 6 6 7 #define INVALID(v) (v. cap_radius< v.radius)7 #define INVALID(v) (v.radius_cap < v.radius) 8 8 9 9 // Integral over a convex lens kernel for t in [h/R,1]. See the docs for … … 12 12 // h is the length of the lens "inside" the cylinder. This negative wrt the 13 13 // definition of h in the docs. 14 // cap_radiusis the radius of the lens14 // radius_cap is the radius of the lens 15 15 // length is the cylinder length, or the separation between the lens halves 16 16 // alpha is the angle of the cylinder wrt q. 17 17 static double 18 _cap_kernel(double q, double h, double cap_radius,18 _cap_kernel(double q, double h, double radius_cap, 19 19 double half_length, double sin_alpha, double cos_alpha) 20 20 { 21 21 // translate a point in [-1,1] to a point in [lower,upper] 22 22 const double upper = 1.0; 23 const double lower = h/ cap_radius; // integral lower bound23 const double lower = h/radius_cap; // integral lower bound 24 24 const double zm = 0.5*(upper-lower); 25 25 const double zb = 0.5*(upper+lower); … … 32 32 // m = q R cos(alpha) 33 33 // b = q(L/2-h) cos(alpha) 34 const double m = q* cap_radius*cos_alpha; // cos argument slope34 const double m = q*radius_cap*cos_alpha; // cos argument slope 35 35 const double b = q*(half_length-h)*cos_alpha; // cos argument intercept 36 const double qrst = q* cap_radius*sin_alpha; // Q*R*sin(theta)36 const double qrst = q*radius_cap*sin_alpha; // Q*R*sin(theta) 37 37 double total = 0.0; 38 38 for (int i=0; i<76 ;i++) { … … 45 45 // translate dx in [-1,1] to dx in [lower,upper] 46 46 const double integral = total*zm; 47 const double cap_Fq = 2*M_PI*cube( cap_radius)*integral;47 const double cap_Fq = 2*M_PI*cube(radius_cap)*integral; 48 48 return cap_Fq; 49 49 } 50 50 51 double form_volume(double radius, double cap_radius, double length)51 double form_volume(double radius, double radius_cap, double length) 52 52 { 53 53 // cap radius should never be less than radius when this is called … … 74 74 // = pi r^2 L + pi hc (r^2 + hc^2/3) 75 75 // = pi (r^2 (L+hc) + hc^3/3) 76 const double hc = cap_radius - sqrt(cap_radius*cap_radius- radius*radius);76 const double hc = radius_cap - sqrt(radius_cap*radius_cap - radius*radius); 77 77 return M_PI*(radius*radius*(length+hc) + hc*hc*hc/3.0); 78 78 } 79 79 80 80 double Iq(double q, double sld, double solvent_sld, 81 double radius, double cap_radius, double length)81 double radius, double radius_cap, double length) 82 82 { 83 const double h = sqrt( cap_radius*cap_radius- radius*radius);83 const double h = sqrt(radius_cap*radius_cap - radius*radius); 84 84 const double half_length = 0.5*length; 85 85 … … 93 93 SINCOS(alpha, sin_alpha, cos_alpha); 94 94 95 const double cap_Fq = _cap_kernel(q, h, cap_radius, half_length, sin_alpha, cos_alpha);95 const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 96 96 const double bj = sas_J1c(q*radius*sin_alpha); 97 97 const double si = sinc(q*half_length*cos_alpha); … … 111 111 double Iqxy(double qx, double qy, 112 112 double sld, double solvent_sld, double radius, 113 double cap_radius, double length,113 double radius_cap, double length, 114 114 double theta, double phi) 115 115 { … … 121 121 const double alpha = acos(cos_val); // rod angle relative to q 122 122 123 const double h = sqrt( cap_radius*cap_radius- radius*radius);123 const double h = sqrt(radius_cap*radius_cap - radius*radius); 124 124 const double half_length = 0.5*length; 125 125 126 126 double sin_alpha, cos_alpha; // slots to hold sincos function output 127 127 SINCOS(alpha, sin_alpha, cos_alpha); 128 const double cap_Fq = _cap_kernel(q, h, cap_radius, half_length, sin_alpha, cos_alpha);128 const double cap_Fq = _cap_kernel(q, h, radius_cap, half_length, sin_alpha, cos_alpha); 129 129 const double bj = sas_J1c(q*radius*sin_alpha); 130 130 const double si = sinc(q*half_length*cos_alpha); -
sasmodels/models/capped_cylinder.py
r42356c8 r2222134 96 96 Note: As the length of cylinder -->0, 97 97 it becomes a Convex Lens. 98 It must be that radius <(=) cap_radius.98 It must be that radius <(=) radius_cap. 99 99 [Parameters]; 100 100 scale: volume fraction of spheres, … … 102 102 radius: radius of the cylinder, 103 103 length: length of the cylinder, 104 cap_radius: radius of the semi-spherical cap,104 radius_cap: radius of the semi-spherical cap, 105 105 sld: SLD of the capped cylinder, 106 106 sld_solvent: SLD of the solvent. … … 122 122 # in the capped cylinder, and zero for no bar in the barbell model. In 123 123 # both models, one would be a pill. 124 [" cap_radius", "Ang", 20, [0, inf], "volume", "Cap radius"],124 ["radius_cap", "Ang", 20, [0, inf], "volume", "Cap radius"], 125 125 ["length", "Ang", 400, [0, inf], "volume", "Cylinder length"], 126 126 ["theta", "degrees", 60, [-inf, inf], "orientation", "In plane angle"], … … 133 133 demo = dict(scale=1, background=0, 134 134 sld=6, sld_solvent=1, 135 radius=260, cap_radius=290, length=290,135 radius=260, radius_cap=290, length=290, 136 136 theta=30, phi=15, 137 137 radius_pd=.2, radius_pd_n=1, 138 cap_radius_pd=.2, cap_radius_pd_n=1,138 radius_cap_pd=.2, radius_cap_pd_n=1, 139 139 length_pd=.2, length_pd_n=10, 140 140 theta_pd=15, theta_pd_n=45, -
sasmodels/models/core_shell_bicelle.c
r43b7eea r2222134 1 double form_volume(double radius, double rim_thickness, double face_thickness, double length);1 double form_volume(double radius, double thick_rim, double thick_face, double length); 2 2 double Iq(double q, 3 3 double radius, 4 double rim_thickness,5 double face_thickness,4 double thick_rim, 5 double thick_face, 6 6 double length, 7 7 double core_sld, … … 13 13 double Iqxy(double qx, double qy, 14 14 double radius, 15 double rim_thickness,16 double face_thickness,15 double thick_rim, 16 double thick_face, 17 17 double length, 18 18 double core_sld, … … 24 24 25 25 26 double form_volume(double radius, double rim_thickness, double face_thickness, double length)26 double form_volume(double radius, double thick_rim, double thick_face, double length) 27 27 { 28 return M_PI*(radius+ rim_thickness)*(radius+rim_thickness)*(length+2*face_thickness);28 return M_PI*(radius+thick_rim)*(radius+thick_rim)*(length+2*thick_face); 29 29 } 30 30 … … 102 102 bicelle_kernel_2d(double q, double q_x, double q_y, 103 103 double radius, 104 double rim_thickness,105 double face_thickness,104 double thick_rim, 105 double thick_face, 106 106 double length, 107 107 double core_sld, … … 125 125 126 126 // Get the kernel 127 double answer = bicelle_kernel(q, radius, rim_thickness, face_thickness,127 double answer = bicelle_kernel(q, radius, thick_rim, thick_face, 128 128 length/2.0, core_sld, face_sld, rim_sld, 129 129 solvent_sld, alpha) / fabs(sin(alpha)); … … 136 136 double Iq(double q, 137 137 double radius, 138 double rim_thickness,139 double face_thickness,138 double thick_rim, 139 double thick_face, 140 140 double length, 141 141 double core_sld, … … 144 144 double solvent_sld) 145 145 { 146 double intensity = bicelle_integration(q, radius, rim_thickness, face_thickness,146 double intensity = bicelle_integration(q, radius, thick_rim, thick_face, 147 147 length, core_sld, face_sld, rim_sld, solvent_sld); 148 148 return intensity*1.0e-4; … … 152 152 double Iqxy(double qx, double qy, 153 153 double radius, 154 double rim_thickness,155 double face_thickness,154 double thick_rim, 155 double thick_face, 156 156 double length, 157 157 double core_sld, … … 166 166 double intensity = bicelle_kernel_2d(q, qx/q, qy/q, 167 167 radius, 168 rim_thickness,169 face_thickness,168 thick_rim, 169 thick_face, 170 170 length, 171 171 core_sld, -
sasmodels/models/core_shell_bicelle.py
r40a87fa r2222134 76 76 parameters = [ 77 77 ["radius", "Ang", 20, [0, inf], "volume", "Cylinder core radius"], 78 [" rim_thickness", "Ang", 10, [0, inf], "volume", "Rim shell thickness"],79 [" face_thickness", "Ang", 10, [0, inf], "volume", "Cylinder face thickness"],78 ["thick_rim", "Ang", 10, [0, inf], "volume", "Rim shell thickness"], 79 ["thick_face", "Ang", 10, [0, inf], "volume", "Cylinder face thickness"], 80 80 ["length", "Ang", 400, [0, inf], "volume", "Cylinder length"], 81 81 ["sld_core", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Cylinder core scattering length density"], … … 94 94 demo = dict(scale=1, background=0, 95 95 radius=20.0, 96 rim_thickness=10.0,97 face_thickness=10.0,96 thick_rim=10.0, 97 thick_face=10.0, 98 98 length=400.0, 99 99 sld_core=1.0, … … 108 108 # Accuracy tests based on content in test/utest_other_models.py 109 109 [{'radius': 20.0, 110 ' rim_thickness': 10.0,111 ' face_thickness': 10.0,110 'thick_rim': 10.0, 111 'thick_face': 10.0, 112 112 'length': 400.0, 113 113 'sld_core': 1.0, … … 119 119 120 120 [{'radius': 20.0, 121 ' rim_thickness': 10.0,122 ' face_thickness': 10.0,121 'thick_rim': 10.0, 122 'thick_face': 10.0, 123 123 'length': 400.0, 124 124 'sld_core': 1.0, … … 133 133 # Additional tests with larger range of parameters 134 134 [{'radius': 3.0, 135 ' rim_thickness': 100.0,136 ' face_thickness': 100.0,135 'thick_rim': 100.0, 136 'thick_face': 100.0, 137 137 'length': 1200.0, 138 138 'sld_core': 5.0, -
sasmodels/models/core_shell_ellipsoid.c
r29172aa r2222134 1 double form_volume(double equat_core,2 double polar_core,3 double equat_shell,4 double polar_shell);1 double form_volume(double radius_equat_core, 2 double radius_polar_core, 3 double radius_equat_shell, 4 double radius_polar_shell); 5 5 double Iq(double q, 6 double equat_core,7 double polar_core,8 double equat_shell,9 double polar_shell,6 double radius_equat_core, 7 double radius_polar_core, 8 double radius_equat_shell, 9 double radius_polar_shell, 10 10 double sld_core, 11 11 double sld_shell, … … 14 14 15 15 double Iqxy(double qx, double qy, 16 double equat_core,17 double polar_core,18 double equat_shell,19 double polar_shell,16 double radius_equat_core, 17 double radius_polar_core, 18 double radius_equat_shell, 19 double radius_polar_shell, 20 20 double sld_core, 21 21 double sld_shell, … … 25 25 26 26 27 double form_volume(double equat_core,28 double polar_core,29 double equat_shell,30 double polar_shell)27 double form_volume(double radius_equat_core, 28 double radius_polar_core, 29 double radius_equat_shell, 30 double radius_polar_shell) 31 31 { 32 double vol = 4.0*M_PI/3.0* equat_shell*equat_shell*polar_shell;32 double vol = 4.0*M_PI/3.0*radius_equat_shell*radius_equat_shell*radius_polar_shell; 33 33 return vol; 34 34 } … … 36 36 static double 37 37 core_shell_ellipsoid_kernel(double q, 38 double equat_core,39 double polar_core,40 double equat_shell,41 double polar_shell,38 double radius_equat_core, 39 double radius_polar_core, 40 double radius_equat_shell, 41 double radius_polar_shell, 42 42 double sld_core, 43 43 double sld_shell, … … 57 57 double zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 58 58 double yyy = Gauss76Wt[i] * gfn4(zi, 59 equat_core,60 polar_core,61 equat_shell,62 polar_shell,59 radius_equat_core, 60 radius_polar_core, 61 radius_equat_shell, 62 radius_polar_shell, 63 63 delpc, 64 64 delps, … … 77 77 static double 78 78 core_shell_ellipsoid_kernel_2d(double q, double q_x, double q_y, 79 double equat_core,80 double polar_core,81 double equat_shell,82 double polar_shell,79 double radius_equat_core, 80 double radius_polar_core, 81 double radius_equat_shell, 82 double radius_polar_shell, 83 83 double sld_core, 84 84 double sld_shell, … … 105 105 // Call the IGOR library function to get the kernel: MUST use gfn4 not gf2 because of the def of params. 106 106 double answer = gfn4(cos_val, 107 equat_core,108 polar_core,109 equat_shell,110 polar_shell,107 radius_equat_core, 108 radius_polar_core, 109 radius_equat_shell, 110 radius_polar_shell, 111 111 sldcs, 112 112 sldss, … … 120 120 121 121 double Iq(double q, 122 double equat_core,123 double polar_core,124 double equat_shell,125 double polar_shell,122 double radius_equat_core, 123 double radius_polar_core, 124 double radius_equat_shell, 125 double radius_polar_shell, 126 126 double sld_core, 127 127 double sld_shell, … … 129 129 { 130 130 double intensity = core_shell_ellipsoid_kernel(q, 131 equat_core,132 polar_core,133 equat_shell,134 polar_shell,131 radius_equat_core, 132 radius_polar_core, 133 radius_equat_shell, 134 radius_polar_shell, 135 135 sld_core, 136 136 sld_shell, … … 142 142 143 143 double Iqxy(double qx, double qy, 144 double equat_core,145 double polar_core,146 double equat_shell,147 double polar_shell,144 double radius_equat_core, 145 double radius_polar_core, 146 double radius_equat_shell, 147 double radius_polar_shell, 148 148 double sld_core, 149 149 double sld_shell, … … 155 155 q = sqrt(qx*qx+qy*qy); 156 156 double intensity = core_shell_ellipsoid_kernel_2d(q, qx/q, qy/q, 157 equat_core,158 polar_core,159 equat_shell,160 polar_shell,157 radius_equat_core, 158 radius_polar_core, 159 radius_equat_shell, 160 radius_polar_shell, 161 161 sld_core, 162 162 sld_shell, -
sasmodels/models/core_shell_ellipsoid.py
r40a87fa r2222134 52 52 .. note:: 53 53 The 2nd virial coefficient of the solid ellipsoid is calculated based on 54 the *radius_a* (= * polar_shell)* and *radius_b (=equat_shell)* values,54 the *radius_a* (= *radius_polar_shell)* and *radius_b (= radius_equat_shell)* values, 55 55 and used as the effective radius for *S(Q)* when $P(Q) * S(Q)$ is applied. 56 56 … … 82 82 single particle scattering amplitude. 83 83 [Parameters]: 84 equat_core = equatorial radius of core, Rminor_core,85 polar_core = polar radius of core, Rmajor_core,86 equat_shell = equatorial radius of shell, Rminor_outer,87 polar_shell = polar radius of shell, Rmajor_outer,84 radius_equat_core = equatorial radius of core, Rminor_core, 85 radius_polar_core = polar radius of core, Rmajor_core, 86 radius_equat_shell = equatorial radius of shell, Rminor_outer, 87 radius_polar_shell = polar radius of shell, Rmajor_outer, 88 88 sld_core = scattering length density of core, 89 89 sld_shell = scattering length density of shell, … … 102 102 # ["name", "units", default, [lower, upper], "type", "description"], 103 103 parameters = [ 104 [" equat_core", "Ang", 200, [0, inf], "volume", "Equatorial radius of core, r minor core"],105 [" polar_core", "Ang", 10, [0, inf], "volume", "Polar radius of core, r major core"],106 [" equat_shell", "Ang", 250, [0, inf], "volume", "Equatorial radius of shell, r minor outer"],107 [" polar_shell", "Ang", 30, [0, inf], "volume", "Polar radius of shell, r major outer"],104 ["radius_equat_core", "Ang", 200, [0, inf], "volume", "Equatorial radius of core, r minor core"], 105 ["radius_polar_core", "Ang", 10, [0, inf], "volume", "Polar radius of core, r major core"], 106 ["radius_equat_shell", "Ang", 250, [0, inf], "volume", "Equatorial radius of shell, r minor outer"], 107 ["radius_polar_shell", "Ang", 30, [0, inf], "volume", "Polar radius of shell, r major outer"], 108 108 ["sld_core", "1e-6/Ang^2", 2, [-inf, inf], "sld", "Core scattering length density"], 109 109 ["sld_shell", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Shell scattering length density"], … … 116 116 source = ["lib/sph_j1c.c", "lib/gfn.c", "lib/gauss76.c", "core_shell_ellipsoid.c"] 117 117 118 def ER( equat_core, polar_core, equat_shell,polar_shell):118 def ER(radius_equat_core, radius_polar_core, radius_equat_shell, radius_polar_shell): 119 119 """ 120 120 Returns the effective radius used in the S*P calculation 121 121 """ 122 122 from .ellipsoid import ER as ellipsoid_ER 123 return ellipsoid_ER( polar_shell,equat_shell)123 return ellipsoid_ER(radius_polar_shell, radius_equat_shell) 124 124 125 125 126 126 demo = dict(scale=1, background=0.001, 127 equat_core=200.0,128 polar_core=10.0,129 equat_shell=250.0,130 polar_shell=30.0,127 radius_equat_core=200.0, 128 radius_polar_core=10.0, 129 radius_equat_shell=250.0, 130 radius_polar_shell=30.0, 131 131 sld_core=2.0, 132 132 sld_shell=1.0, … … 142 142 tests = [ 143 143 # Accuracy tests based on content in test/utest_other_models.py 144 [{' equat_core': 200.0,145 ' polar_core': 20.0,146 ' equat_shell': 250.0,147 ' polar_shell': 30.0,144 [{'radius_equat_core': 200.0, 145 'radius_polar_core': 20.0, 146 'radius_equat_shell': 250.0, 147 'radius_polar_shell': 30.0, 148 148 'sld_core': 2.0, 149 149 'sld_shell': 1.0, … … 156 156 [{'background': 0.01}, 0.1, 8.86741], 157 157 158 [{' equat_core': 20.0,159 ' polar_core': 200.0,160 ' equat_shell': 54.0,161 ' polar_shell': 3.0,158 [{'radius_equat_core': 20.0, 159 'radius_polar_core': 200.0, 160 'radius_equat_shell': 54.0, 161 'radius_polar_shell': 3.0, 162 162 'sld_core': 20.0, 163 163 'sld_shell': 10.0, … … 169 169 [{'background': 0.001}, (0.4, 0.5), 0.00170471], 170 170 171 [{' equat_core': 20.0,172 ' polar_core': 200.0,173 ' equat_shell': 54.0,174 ' polar_shell': 3.0,171 [{'radius_equat_core': 20.0, 172 'radius_polar_core': 200.0, 173 'radius_equat_shell': 54.0, 174 'radius_polar_shell': 3.0, 175 175 'sld_core': 20.0, 176 176 'sld_shell': 10.0, -
sasmodels/models/core_shell_ellipsoid_xt.c
re7678b2 r2222134 1 double form_volume(double equat_core,1 double form_volume(double radius_equat_core, 2 2 double polar_core, 3 3 double equat_shell, 4 4 double polar_shell); 5 5 double Iq(double q, 6 double equat_core,6 double radius_equat_core, 7 7 double x_core, 8 double t _shell,8 double thick_shell, 9 9 double x_polar_shell, 10 10 double core_sld, … … 14 14 15 15 double Iqxy(double qx, double qy, 16 double equat_core,16 double radius_equat_core, 17 17 double x_core, 18 double t _shell,18 double thick_shell, 19 19 double x_polar_shell, 20 20 double core_sld, … … 25 25 26 26 27 double form_volume(double equat_core,27 double form_volume(double radius_equat_core, 28 28 double x_core, 29 double t _shell,29 double thick_shell, 30 30 double x_polar_shell) 31 31 { 32 const double equat_shell = equat_core + t_shell;33 const double polar_shell = equat_core*x_core + t_shell*x_polar_shell;32 const double equat_shell = radius_equat_core + thick_shell; 33 const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 34 34 double vol = 4.0*M_PI/3.0*equat_shell*equat_shell*polar_shell; 35 35 return vol; … … 38 38 static double 39 39 core_shell_ellipsoid_xt_kernel(double q, 40 double equat_core,40 double radius_equat_core, 41 41 double x_core, 42 double t _shell,42 double thick_shell, 43 43 double x_polar_shell, 44 44 double core_sld, … … 55 55 56 56 57 const double polar_core = equat_core*x_core;58 const double equat_shell = equat_core + t_shell;59 const double polar_shell = equat_core*x_core + t_shell*x_polar_shell;57 const double polar_core = radius_equat_core*x_core; 58 const double equat_shell = radius_equat_core + thick_shell; 59 const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 60 60 61 61 for(int i=0;i<N_POINTS_76;i++) { 62 62 double zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 63 63 double yyy = Gauss76Wt[i] * gfn4(zi, 64 equat_core,64 radius_equat_core, 65 65 polar_core, 66 66 equat_shell, … … 81 81 static double 82 82 core_shell_ellipsoid_xt_kernel_2d(double q, double q_x, double q_y, 83 double equat_core,83 double radius_equat_core, 84 84 double x_core, 85 double t _shell,85 double thick_shell, 86 86 double x_polar_shell, 87 87 double core_sld, … … 106 106 const double cos_val = cyl_x*q_x + cyl_y*q_y; 107 107 108 const double polar_core = equat_core*x_core;109 const double equat_shell = equat_core + t_shell;110 const double polar_shell = equat_core*x_core + t_shell*x_polar_shell;108 const double polar_core = radius_equat_core*x_core; 109 const double equat_shell = radius_equat_core + thick_shell; 110 const double polar_shell = radius_equat_core*x_core + thick_shell*x_polar_shell; 111 111 112 112 // Call the IGOR library function to get the kernel: 113 113 // MUST use gfn4 not gf2 because of the def of params. 114 114 double answer = gfn4(cos_val, 115 equat_core,115 radius_equat_core, 116 116 polar_core, 117 117 equat_shell, … … 128 128 129 129 double Iq(double q, 130 double equat_core,130 double radius_equat_core, 131 131 double x_core, 132 double t _shell,132 double thick_shell, 133 133 double x_polar_shell, 134 134 double core_sld, … … 137 137 { 138 138 double intensity = core_shell_ellipsoid_xt_kernel(q, 139 equat_core,139 radius_equat_core, 140 140 x_core, 141 t _shell,141 thick_shell, 142 142 x_polar_shell, 143 143 core_sld, … … 150 150 151 151 double Iqxy(double qx, double qy, 152 double equat_core,152 double radius_equat_core, 153 153 double x_core, 154 double t _shell,154 double thick_shell, 155 155 double x_polar_shell, 156 156 double core_sld, … … 163 163 q = sqrt(qx*qx+qy*qy); 164 164 double intensity = core_shell_ellipsoid_xt_kernel_2d(q, qx/q, qy/q, 165 equat_core,165 radius_equat_core, 166 166 x_core, 167 t _shell,167 thick_shell, 168 168 x_polar_shell, 169 169 core_sld, -
sasmodels/models/core_shell_ellipsoid_xt.py
r40a87fa r2222134 14 14 The geometric parameters of this model are 15 15 16 * equat_core =* equatorial core radius *= R_minor_core*16 *radius_equat_core =* equatorial core radius *= R_minor_core* 17 17 18 *X_core = polar_core / equat_core = Rmajor_core / Rminor_core*18 *X_core = polar_core / radius_equat_core = Rmajor_core / Rminor_core* 19 19 20 *T _shell = equat_outer -equat_core = Rminor_outer - Rminor_core*20 *Thick_shell = equat_outer - radius_equat_core = Rminor_outer - Rminor_core* 21 21 22 *XpolarShell = Tpolar_shell / T _shell = (Rmajor_outer - Rmajor_core)/22 *XpolarShell = Tpolar_shell / Thick_shell = (Rmajor_outer - Rmajor_core)/ 23 23 (Rminor_outer - Rminor_core)* 24 24 25 25 In terms of the original radii 26 26 27 *polar_core = equat_core * X_core*27 *polar_core = radius_equat_core * X_core* 28 28 29 *equat_shell = equat_core + T_shell*29 *equat_shell = radius_equat_core + Thick_shell* 30 30 31 *polar_shell = equat_core * X_core + T_shell * XpolarShell*31 *polar_shell = radius_equat_core * X_core + Thick_shell * XpolarShell* 32 32 33 33 (where we note that "shell" perhaps confusingly, relates to the outer radius) … … 68 68 single particle scattering amplitude. 69 69 [Parameters]: 70 equat_core = equatorial radius of core,70 radius_equat_core = equatorial radius of core, 71 71 x_core = ratio of core polar/equatorial radii, 72 t _shell = equatorial radius of outer surface,72 thick_shell = equatorial radius of outer surface, 73 73 x_polar_shell = ratio of polar shell thickness to equatorial shell thickness, 74 74 sld_core = SLD_core … … 88 88 # ["name", "units", default, [lower, upper], "type", "description"], 89 89 parameters = [ 90 [" equat_core", "Ang",20, [0, inf], "volume", "Equatorial radius of core"],90 ["radius_equat_core","Ang", 20, [0, inf], "volume", "Equatorial radius of core"], 91 91 ["x_core", "None", 3, [0, inf], "volume", "axial ratio of core, X = r_polar/r_equatorial"], 92 ["t _shell","Ang", 30, [0, inf], "volume", "thickness of shell at equator"],92 ["thick_shell", "Ang", 30, [0, inf], "volume", "thickness of shell at equator"], 93 93 ["x_polar_shell", "", 1, [0, inf], "volume", "ratio of thickness of shell at pole to that at equator"], 94 94 ["sld_core", "1e-6/Ang^2", 2, [-inf, inf], "sld", "Core scattering length density"], … … 103 103 "core_shell_ellipsoid_xt.c"] 104 104 105 def ER( equat_core, x_core, t_shell, x_polar_shell):105 def ER(radius_equat_core, x_core, thick_shell, x_polar_shell): 106 106 """ 107 107 Returns the effective radius used in the S*P calculation 108 108 """ 109 109 from .ellipsoid import ER as ellipsoid_ER 110 polar_outer = equat_core*x_core + t_shell*x_polar_shell111 equat_outer = equat_core + t_shell110 polar_outer = radius_equat_core*x_core + thick_shell*x_polar_shell 111 equat_outer = radius_equat_core + thick_shell 112 112 return ellipsoid_ER(polar_outer, equat_outer) 113 113 114 114 115 115 demo = dict(scale=0.05, background=0.001, 116 equat_core=20.0,116 radius_equat_core=20.0, 117 117 x_core=3.0, 118 t _shell=30.0,118 thick_shell=30.0, 119 119 x_polar_shell=1.0, 120 120 sld_core=2.0, … … 131 131 tests = [ 132 132 # Accuracy tests based on content in test/utest_coreshellellipsoidXTmodel.py 133 [{' equat_core': 200.0,133 [{'radius_equat_core': 200.0, 134 134 'x_core': 0.1, 135 't _shell': 50.0,135 'thick_shell': 50.0, 136 136 'x_polar_shell': 0.2, 137 137 'sld_core': 2.0, … … 145 145 [{'background': 0.01}, 0.1, 11.6915], 146 146 147 [{' equat_core': 20.0,147 [{'radius_equat_core': 20.0, 148 148 'x_core': 200.0, 149 't _shell': 54.0,149 'thick_shell': 54.0, 150 150 'x_polar_shell': 3.0, 151 151 'sld_core': 20.0, … … 158 158 [{'background': 0.001}, (0.4, 0.5), 0.00690673], 159 159 160 [{' equat_core': 20.0,160 [{'radius_equat_core': 20.0, 161 161 'x_core': 200.0, 162 't _shell': 54.0,162 'thick_shell': 54.0, 163 163 'x_polar_shell': 3.0, 164 164 'sld_core': 20.0, -
sasmodels/models/core_shell_parallelepiped.c
rabdd01c r2222134 1 double form_volume(double a_side, double b_side, double c_side,2 double arim_thickness, double brim_thickness, double crim_thickness);1 double form_volume(double length_a, double length_b, double length_c, 2 double thick_rim_a, double thick_rim_b, double thick_rim_c); 3 3 double Iq(double q, double core_sld, double arim_sld, double brim_sld, double crim_sld, 4 double solvent_sld, double a_side, double b_side, double c_side,5 double arim_thickness, double brim_thickness, double crim_thickness);4 double solvent_sld, double length_a, double length_b, double length_c, 5 double thick_rim_a, double thick_rim_b, double thick_rim_c); 6 6 double Iqxy(double qx, double qy, double core_sld, double arim_sld, double brim_sld, 7 double crim_sld, double solvent_sld, double a_side, double b_side,8 double c_side, double arim_thickness, double brim_thickness,9 double crim_thickness, double theta, double phi, double psi);10 11 double form_volume(double a_side, double b_side, double c_side,12 double arim_thickness, double brim_thickness, double crim_thickness)7 double crim_sld, double solvent_sld, double length_a, double length_b, 8 double length_c, double thick_rim_a, double thick_rim_b, 9 double thick_rim_c, double theta, double phi, double psi); 10 11 double form_volume(double length_a, double length_b, double length_c, 12 double thick_rim_a, double thick_rim_b, double thick_rim_c) 13 13 { 14 //return a_side * b_side * c_side;15 return a_side * b_side * c_side+16 2.0 * arim_thickness * b_side * c_side+17 2.0 * brim_thickness * a_side * c_side+18 2.0 * crim_thickness * a_side * b_side;14 //return length_a * length_b * length_c; 15 return length_a * length_b * length_c + 16 2.0 * thick_rim_a * length_b * length_c + 17 2.0 * thick_rim_b * length_a * length_c + 18 2.0 * thick_rim_c * length_a * length_b; 19 19 } 20 20 … … 25 25 double crim_sld, 26 26 double solvent_sld, 27 double a_side,28 double b_side,29 double c_side,30 double arim_thickness,31 double brim_thickness,32 double crim_thickness)27 double length_a, 28 double length_b, 29 double length_c, 30 double thick_rim_a, 31 double thick_rim_b, 32 double thick_rim_c) 33 33 { 34 34 // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled … … 36 36 37 37 double t1, t2, t3, t4, tmp, answer; 38 double mu = q * b_side;38 double mu = q * length_b; 39 39 40 40 //calculate volume before rescaling (in original code, but not used) 41 //double vol = form_volume( a_side, b_side, c_side, arim_thickness, brim_thickness, crim_thickness);42 //double vol = a_side * b_side * c_side+43 // 2.0 * arim_thickness * b_side * c_side+44 // 2.0 * brim_thickness * a_side * c_side+45 // 2.0 * crim_thickness * a_side * b_side;41 //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 42 //double vol = length_a * length_b * length_c + 43 // 2.0 * thick_rim_a * length_b * length_c + 44 // 2.0 * thick_rim_b * length_a * length_c + 45 // 2.0 * thick_rim_c * length_a * length_b; 46 46 47 47 // Scale sides by B 48 double a_scaled = a_side / b_side;49 double c_scaled = c_side / b_side;48 double a_scaled = length_a / length_b; 49 double c_scaled = length_c / length_b; 50 50 51 51 // DelRho values (note that drC is not used later) … … 74 74 double mudum = mu * sqrt(1.0-sigma*sigma); 75 75 76 double Vin = a_side * b_side * c_side;77 //double Vot = ( a_side * b_side * c_side+78 // 2.0 * arim_thickness * b_side * c_side+79 // 2.0 * a_side * brim_thickness * c_side+80 // 2.0 * a_side * b_side * crim_thickness);81 double V1 = (2.0 * arim_thickness * b_side * c_side); // incorrect V1 (aa*bb*cc+2*ta*bb*cc)82 double V2 = (2.0 * a_side * brim_thickness * c_side); // incorrect V2(aa*bb*cc+2*aa*tb*cc)76 double Vin = length_a * length_b * length_c; 77 //double Vot = (length_a * length_b * length_c + 78 // 2.0 * thick_rim_a * length_b * length_c + 79 // 2.0 * length_a * thick_rim_b * length_c + 80 // 2.0 * length_a * length_b * thick_rim_c); 81 double V1 = (2.0 * thick_rim_a * length_b * length_c); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 82 double V2 = (2.0 * length_a * thick_rim_b * length_c); // incorrect V2(aa*bb*cc+2*aa*tb*cc) 83 83 84 84 // ta and tb correspond to the definitions in CSPPKernel, but they don't make sense to me (MG) … … 86 86 // while in cspkernel in csparallelepiped.cpp (used for the 2D), all the definitions 87 87 // for ta, tb, tc use also A + 2*rim_thickness (but not scaled by B!!!) 88 double ta = (a_scaled+2.0* arim_thickness)/b_side;89 double tb = (a_scaled+2.0* brim_thickness)/b_side;88 double ta = (a_scaled+2.0*thick_rim_a)/length_b; 89 double tb = (a_scaled+2.0*thick_rim_b)/length_b; 90 90 91 91 double arg1 = (0.5*mudum*a_scaled) * sin(0.5*M_PI*uu); … … 160 160 double crim_sld, 161 161 double solvent_sld, 162 double a_side,163 double b_side,164 double c_side,165 double arim_thickness,166 double brim_thickness,167 double crim_thickness,162 double length_a, 163 double length_b, 164 double length_c, 165 double thick_rim_a, 166 double thick_rim_b, 167 double thick_rim_c, 168 168 double theta, 169 169 double phi, … … 217 217 double drB = brim_sld-solvent_sld; 218 218 double drC = crim_sld-solvent_sld; 219 double Vin = a_side * b_side * c_side;219 double Vin = length_a * length_b * length_c; 220 220 // As for the 1D case, Vot is not used 221 //double Vot = ( a_side * b_side * c_side+222 // 2.0 * arim_thickness * b_side * c_side+223 // 2.0 * a_side * brim_thickness * c_side+224 // 2.0 * a_side * b_side * crim_thickness);225 double V1 = (2.0 * arim_thickness * b_side * c_side); // incorrect V1 (aa*bb*cc+2*ta*bb*cc)226 double V2 = (2.0 * a_side * brim_thickness * c_side); // incorrect V2(aa*bb*cc+2*aa*tb*cc)227 double V3 = (2.0 * a_side * b_side * crim_thickness);221 //double Vot = (length_a * length_b * length_c + 222 // 2.0 * thick_rim_a * length_b * length_c + 223 // 2.0 * length_a * thick_rim_b * length_c + 224 // 2.0 * length_a * length_b * thick_rim_c); 225 double V1 = (2.0 * thick_rim_a * length_b * length_c); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 226 double V2 = (2.0 * length_a * thick_rim_b * length_c); // incorrect V2(aa*bb*cc+2*aa*tb*cc) 227 double V3 = (2.0 * length_a * length_b * thick_rim_c); 228 228 // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 229 // the scaling by B. The use of a_sidefor the 3 of them seems clearly a mistake to me,229 // the scaling by B. The use of length_a for the 3 of them seems clearly a mistake to me, 230 230 // but for the moment I let it like this until understanding better the code. 231 double ta = a_side + 2.0*arim_thickness;232 double tb = a_side + 2.0*brim_thickness;233 double tc = a_side + 2.0*crim_thickness;231 double ta = length_a + 2.0*thick_rim_a; 232 double tb = length_a + 2.0*thick_rim_b; 233 double tc = length_a + 2.0*thick_rim_c; 234 234 //handle arg=0 separately, as sin(t)/t -> 1 as t->0 235 double argA = 0.5*q* a_side*cos_val_a;236 double argB = 0.5*q* b_side*cos_val_b;237 double argC = 0.5*q* c_side*cos_val_c;235 double argA = 0.5*q*length_a*cos_val_a; 236 double argB = 0.5*q*length_b*cos_val_b; 237 double argC = 0.5*q*length_c*cos_val_c; 238 238 double argtA = 0.5*q*ta*cos_val_a; 239 239 double argtB = 0.5*q*tb*cos_val_b; -
sasmodels/models/core_shell_parallelepiped.py
r42356c8 r2222134 124 124 ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "sld", 125 125 "Solvent scattering length density"], 126 [" a_side", "Ang", 35, [0, inf], "volume",126 ["length_a", "Ang", 35, [0, inf], "volume", 127 127 "Shorter side of the parallelepiped"], 128 [" b_side", "Ang", 75, [0, inf], "volume",128 ["length_b", "Ang", 75, [0, inf], "volume", 129 129 "Second side of the parallelepiped"], 130 [" c_side", "Ang", 400, [0, inf], "volume",130 ["length_c", "Ang", 400, [0, inf], "volume", 131 131 "Larger side of the parallelepiped"], 132 [" arim_thickness", "Ang", 10, [0, inf], "volume",132 ["thick_rim_a", "Ang", 10, [0, inf], "volume", 133 133 "Thickness of A rim"], 134 [" brim_thickness", "Ang", 10, [0, inf], "volume",134 ["thick_rim_b", "Ang", 10, [0, inf], "volume", 135 135 "Thickness of B rim"], 136 [" crim_thickness", "Ang", 10, [0, inf], "volume",136 ["thick_rim_c", "Ang", 10, [0, inf], "volume", 137 137 "Thickness of C rim"], 138 138 ["theta", "degrees", 0, [-inf, inf], "orientation", … … 147 147 148 148 149 def ER( a_side, b_side, c_side, arim_thickness, brim_thickness, crim_thickness):149 def ER(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c): 150 150 """ 151 151 Return equivalent radius (ER) … … 153 153 154 154 # surface average radius (rough approximation) 155 surf_rad = sqrt(( a_side + 2.0*arim_thickness) * (b_side + 2.0*brim_thickness) / pi)155 surf_rad = sqrt((length_a + 2.0*thick_rim_a) * (length_b + 2.0*thick_rim_b) / pi) 156 156 157 height = c_side + 2.0*crim_thickness157 height = length_c + 2.0*thick_rim_c 158 158 159 159 ddd = 0.75 * surf_rad * (2 * surf_rad * height + (height + surf_rad) * (height + pi * surf_rad)) … … 166 166 sld_core=1e-6, sld_a=2e-6, sld_b=4e-6, 167 167 sld_c=2e-6, sld_solvent=6e-6, 168 a_side=35, b_side=75, c_side=400,169 arim_thickness=10, brim_thickness=10, crim_thickness=10,168 length_a=35, length_b=75, length_c=400, 169 thick_rim_a=10, thick_rim_b=10, thick_rim_c=10, 170 170 theta=0, phi=0, psi=0, 171 a_side_pd=0.1, a_side_pd_n=1,172 b_side_pd=0.1, b_side_pd_n=1,173 c_side_pd=0.1, c_side_pd_n=1,174 arim_thickness_pd=0.1, arim_thickness_pd_n=1,175 brim_thickness_pd=0.1, brim_thickness_pd_n=1,176 crim_thickness_pd=0.1, crim_thickness_pd_n=1,171 length_a_pd=0.1, length_a_pd_n=1, 172 length_b_pd=0.1, length_b_pd_n=1, 173 length_c_pd=0.1, length_c_pd_n=1, 174 thick_rim_a_pd=0.1, thick_rim_a_pd_n=1, 175 thick_rim_b_pd=0.1, thick_rim_b_pd_n=1, 176 thick_rim_c_pd=0.1, thick_rim_c_pd_n=1, 177 177 theta_pd=10, theta_pd_n=1, 178 178 phi_pd=10, phi_pd_n=1,
Note: See TracChangeset
for help on using the changeset viewer.