Changeset a807206 in sasmodels for sasmodels/models/ellipsoid.c
- Timestamp:
- Sep 30, 2016 10:42:06 PM (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:
- caddb14, 5031ca3
- Parents:
- 2222134
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/ellipsoid.c
r50e1e40 ra807206 1 double form_volume(double r polar, double requatorial);2 double Iq(double q, double sld, double s olvent_sld, double rpolar, double requatorial);3 double Iqxy(double qx, double qy, double sld, double s olvent_sld,4 double r polar, double requatorial, double theta, double phi);1 double form_volume(double radius_polar, double radius_equatorial); 2 double Iq(double q, double sld, double sld_solvent, double radius_polar, double radius_equatorial); 3 double Iqxy(double qx, double qy, double sld, double sld_solvent, 4 double radius_polar, double radius_equatorial, double theta, double phi); 5 5 6 double _ellipsoid_kernel(double q, double r polar, double requatorial, double sin_alpha);7 double _ellipsoid_kernel(double q, double r polar, double requatorial, double sin_alpha)6 double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha); 7 double _ellipsoid_kernel(double q, double radius_polar, double radius_equatorial, double sin_alpha) 8 8 { 9 double ratio = r polar/requatorial;10 const double u = q*r equatorial*sqrt(1.09 double ratio = radius_polar/radius_equatorial; 10 const double u = q*radius_equatorial*sqrt(1.0 11 11 + sin_alpha*sin_alpha*(ratio*ratio - 1.0)); 12 12 const double f = sph_j1c(u); … … 15 15 } 16 16 17 double form_volume(double r polar, double requatorial)17 double form_volume(double radius_polar, double radius_equatorial) 18 18 { 19 return M_4PI_3*r polar*requatorial*requatorial;19 return M_4PI_3*radius_polar*radius_equatorial*radius_equatorial; 20 20 } 21 21 22 22 double Iq(double q, 23 23 double sld, 24 double s olvent_sld,25 double r polar,26 double r equatorial)24 double sld_solvent, 25 double radius_polar, 26 double radius_equatorial) 27 27 { 28 28 // translate a point in [-1,1] to a point in [0, 1] … … 33 33 //const double sin_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 34 34 const double sin_alpha = Gauss76Z[i]*zm + zb; 35 total += Gauss76Wt[i] * _ellipsoid_kernel(q, r polar, requatorial, sin_alpha);35 total += Gauss76Wt[i] * _ellipsoid_kernel(q, radius_polar, radius_equatorial, sin_alpha); 36 36 } 37 37 // translate dx in [-1,1] to dx in [lower,upper] 38 38 const double form = total*zm; 39 const double s = (sld - s olvent_sld) * form_volume(rpolar, requatorial);39 const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 40 40 return 1.0e-4 * s * s * form; 41 41 } … … 43 43 double Iqxy(double qx, double qy, 44 44 double sld, 45 double s olvent_sld,46 double r polar,47 double r equatorial,45 double sld_solvent, 46 double radius_polar, 47 double radius_equatorial, 48 48 double theta, 49 49 double phi) … … 55 55 // TODO: check if this is actually sin(alpha), not cos(alpha) 56 56 const double cos_alpha = cn*cos(phi*M_PI_180)*(qx/q) + sn*(qy/q); 57 const double form = _ellipsoid_kernel(q, r polar, requatorial, cos_alpha);58 const double s = (sld - s olvent_sld) * form_volume(rpolar, requatorial);57 const double form = _ellipsoid_kernel(q, radius_polar, radius_equatorial, cos_alpha); 58 const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 59 59 60 60 return 1.0e-4 * form * s * s;
Note: See TracChangeset
for help on using the changeset viewer.