Changeset 68dd6a9 in sasmodels
 Timestamp:
 Mar 22, 2017 11:46:16 PM (7 years ago)
 Branches:
 master, core_shell_microgels, costrafo411, magnetic_model, ticket1257vesicleproduct, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
 Children:
 d86967b
 Parents:
 01ea374
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sasmodels/models/triaxial_ellipsoid.c
r925ad6e r68dd6a9 20 20 double radius_polar) 21 21 { 22 double sn, cn; 23 // translate a point in [1,1] to a point in [0, 1] 24 const double zm = 0.5; 25 const double zb = 0.5; 22 const double pa = square(radius_equat_minor/radius_equat_major)  1.0; 23 const double pc = square(radius_polar/radius_equat_major)  1.0; 24 // translate a point in [1,1] to a point in [0, pi/2] 25 const double zm = M_PI_4; 26 const double zb = M_PI_4; 26 27 double outer = 0.0; 27 28 for (int i=0;i<76;i++) { 28 //const double cos_alpha = (Gauss76Z[i]*(upperlower) + upper + lower)/2; 29 const double x = 0.5*(Gauss76Z[i] + 1.0); 30 SINCOS(M_PI_2*x, sn, cn); 31 const double acosx2 = radius_equat_minor*radius_equat_minor*cn*cn; 32 const double bsinx2 = radius_equat_major*radius_equat_major*sn*sn; 33 const double c2 = radius_polar*radius_polar; 29 //const double u = Gauss76Z[i]*(upperlower)/2 + (upper + lower)/2; 30 const double phi = Gauss76Z[i]*zm + zb; 31 const double pa_sinsq_phi = pa*square(sin(phi)); 34 32 35 33 double inner = 0.0; 34 const double um = 0.5; 35 const double ub = 0.5; 36 36 for (int j=0;j<76;j++) { 37 const double ysq = square(Gauss76Z[j]*zm + zb); 38 const double t = q*sqrt(acosx2 + bsinx2*(1.0ysq) + c2*ysq); 39 const double fq = sas_3j1x_x(t); 40 inner += Gauss76Wt[j] * fq * fq ; 37 // translate a point in [1,1] to a point in [0, 1] 38 const double usq = square(Gauss76Z[j]*um + ub); 39 const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0usq) + 1.0 + pc*usq); 40 const double fq = sas_3j1x_x(q*r); 41 inner += Gauss76Wt[j] * fq * fq; 41 42 } 42 outer += Gauss76Wt[i] * 0.5 * inner;43 outer += Gauss76Wt[i] * inner; // correcting for dx later 43 44 } 44 // translate dx in [1,1] to dx in [lower,upper]45 const double fqsq = outer *zm;45 // translate integration ranges from [1,1] to [lower,upper] and normalize by 4 pi 46 const double fqsq = outer/4.0; // = outer*um*zm*8.0/(4.0*M_PI); 46 47 const double s = (sld  sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 47 48 return 1.0e4 * s * s * fqsq; … … 58 59 double psi) 59 60 { 60 double q, calpha, cmu, cnu;61 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, calpha, cmu, cnu);61 double q, xhat, yhat, zhat; 62 ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat); 62 63 63 const double t = q*sqrt(radius_equat_minor*radius_equat_minor*cnu*cnu64 + radius_equat_major*radius_equat_major*cmu*cmu65 + radius_polar*radius_polar*calpha*calpha);66 const double fq = sas_3j1x_x( t);64 const double r = sqrt(square(radius_equat_minor*xhat) 65 + square(radius_equat_major*yhat) 66 + square(radius_polar*zhat)); 67 const double fq = sas_3j1x_x(q*r); 67 68 const double s = (sld  sld_solvent) * form_volume(radius_equat_minor, radius_equat_major, radius_polar); 68 69
Note: See TracChangeset
for help on using the changeset viewer.