Changeset 3b571ae in sasmodels for sasmodels/models/ellipsoid.c
- Timestamp:
- Mar 22, 2017 4:32:33 PM (7 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:
- 61104c8
- Parents:
- b00a646
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/ellipsoid.c
r130d4c7 r3b571ae 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 double7 _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 the11 // 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 be18 // slightly better for 2D which has already computed it, but it introduces19 // an extra sqrt and square for 1-D not required by the current form, so20 // leave it as is.21 const double r = radius_equatorial22 * 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 }27 5 28 6 double form_volume(double radius_polar, double radius_equatorial) … … 37 15 double radius_equatorial) 38 16 { 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 dT 19 // = int_0^pi/2 Phi^2(h a sqrt((1-sin^2) + v^2 sin^2) cos dT 20 // = int_0^pi/2 Phi^2(h a sqrt(1 + sin^2(v^2-1)) cos dT 21 // u-substitution of 22 // u = sin, du = cos dT 23 // i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du 24 const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0; 25 39 26 // translate a point in [-1,1] to a point in [0, 1] 27 // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2; 40 28 const double zm = 0.5; 41 29 const double zb = 0.5; 42 30 double total = 0.0; 43 31 for (int i=0;i<76;i++) { 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); 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; 47 36 } 48 37 // translate dx in [-1,1] to dx in [lower,upper] … … 62 51 double q, sin_alpha, cos_alpha; 63 52 ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sin_alpha, cos_alpha); 64 const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, 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); 65 56 const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 66 57 67 return 1.0e-4 * form * s * s;58 return 1.0e-4 * square(f * s); 68 59 } 69 60
Note: See TracChangeset
for help on using the changeset viewer.