Changeset 2222134 in sasmodels for sasmodels/models/capped_cylinder.c
- Timestamp:
- Sep 30, 2016 9:07:16 AM (7 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
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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);
Note: See TracChangeset
for help on using the changeset viewer.