Changeset 994d77f in sasmodels for sasmodels/models/triaxial_ellipsoid.c
- Timestamp:
- Oct 30, 2014 10:33:53 AM (9 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:
- ef2861b
- Parents:
- d087487b
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/triaxial_ellipsoid.c
r5d4777d r994d77f 1 real form_volume(real req_minor, real req_major, realrpolar);2 real Iq(real q, real sld, realsolvent_sld,3 real req_minor, real req_major, realrpolar);4 real Iqxy(real qx, real qy, real sld, realsolvent_sld,5 real req_minor, real req_major, real rpolar, real theta, real phi, realpsi);1 double form_volume(double req_minor, double req_major, double rpolar); 2 double Iq(double q, double sld, double solvent_sld, 3 double req_minor, double req_major, double rpolar); 4 double Iqxy(double qx, double qy, double sld, double solvent_sld, 5 double req_minor, double req_major, double rpolar, double theta, double phi, double psi); 6 6 7 real form_volume(real req_minor, real req_major, realrpolar)7 double form_volume(double req_minor, double req_major, double rpolar) 8 8 { 9 return REAL(1.333333333333333)*M_PI*req_minor*req_major*rpolar;9 return 1.333333333333333*M_PI*req_minor*req_major*rpolar; 10 10 } 11 11 12 real Iq(realq,13 realsld,14 realsolvent_sld,15 realreq_minor,16 realreq_major,17 realrpolar)12 double Iq(double q, 13 double sld, 14 double solvent_sld, 15 double req_minor, 16 double req_major, 17 double rpolar) 18 18 { 19 // if (req_minor > req_major || req_major > rpolar) return REAL(-1.0); // Exclude invalid region19 // if (req_minor > req_major || req_major > rpolar) return -1.0; // Exclude invalid region 20 20 21 realsn, cn;22 realst, ct;23 //const real lower = REAL(0.0);24 //const real upper = REAL(1.0);25 real outer = REAL(0.0);21 double sn, cn; 22 double st, ct; 23 //const double lower = 0.0; 24 //const double upper = 1.0; 25 double outer = 0.0; 26 26 for (int i=0;i<76;i++) { 27 //const realcos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2;28 const real x = REAL(0.5)*(Gauss76Z[i] + REAL(1.0));27 //const double cos_alpha = (Gauss76Z[i]*(upper-lower) + upper + lower)/2; 28 const double x = 0.5*(Gauss76Z[i] + 1.0); 29 29 SINCOS(M_PI_2*x, sn, cn); 30 const realacosx2 = req_minor*req_minor*cn*cn;31 const realbsinx2 = req_major*req_major*sn*sn;32 const realc2 = rpolar*rpolar;30 const double acosx2 = req_minor*req_minor*cn*cn; 31 const double bsinx2 = req_major*req_major*sn*sn; 32 const double c2 = rpolar*rpolar; 33 33 34 real inner = REAL(0.0);34 double inner = 0.0; 35 35 for (int j=0;j<76;j++) { 36 const real y = REAL(0.5)*(Gauss76Z[j] + REAL(1.0));37 const real t = q*sqrt(acosx2 + bsinx2*(REAL(1.0)-y*y) + c2*y*y);36 const double y = 0.5*(Gauss76Z[j] + 1.0); 37 const double t = q*sqrt(acosx2 + bsinx2*(1.0-y*y) + c2*y*y); 38 38 SINCOS(t, st, ct); 39 const real fq = ( t==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(st-t*ct)/(t*t*t) );39 const double fq = ( t==0.0 ? 1.0 : 3.0*(st-t*ct)/(t*t*t) ); 40 40 inner += Gauss76Wt[j] * fq * fq ; 41 41 } 42 outer += Gauss76Wt[i] * REAL(0.5)* inner;42 outer += Gauss76Wt[i] * 0.5 * inner; 43 43 } 44 //const realfq2 = (upper-lower)/2*outer;45 const real fq2 = REAL(0.5)*outer;46 const reals = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar);47 return REAL(1.0e-4)* fq2 * s * s;44 //const double fq2 = (upper-lower)/2*outer; 45 const double fq2 = 0.5*outer; 46 const double s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar); 47 return 1.0e-4 * fq2 * s * s; 48 48 } 49 49 50 real Iqxy(real qx, realqy,51 realsld,52 realsolvent_sld,53 realreq_minor,54 realreq_major,55 realrpolar,56 realtheta,57 realphi,58 realpsi)50 double Iqxy(double qx, double qy, 51 double sld, 52 double solvent_sld, 53 double req_minor, 54 double req_major, 55 double rpolar, 56 double theta, 57 double phi, 58 double psi) 59 59 { 60 // if (req_minor > req_major || req_major > rpolar) return REAL(-1.0); // Exclude invalid region60 // if (req_minor > req_major || req_major > rpolar) return -1.0; // Exclude invalid region 61 61 62 realstheta, ctheta;63 realsphi, cphi;64 realspsi, cpsi;65 realst, ct;62 double stheta, ctheta; 63 double sphi, cphi; 64 double spsi, cpsi; 65 double st, ct; 66 66 67 const realq = sqrt(qx*qx + qy*qy);68 const realqxhat = qx/q;69 const realqyhat = qy/q;67 const double q = sqrt(qx*qx + qy*qy); 68 const double qxhat = qx/q; 69 const double qyhat = qy/q; 70 70 SINCOS(theta*M_PI_180, stheta, ctheta); 71 71 SINCOS(phi*M_PI_180, sphi, cphi); 72 72 SINCOS(psi*M_PI_180, spsi, cpsi); 73 const realcalpha = ctheta*cphi*qxhat + stheta*qyhat;74 const realcnu = (-cphi*spsi*stheta + sphi*cpsi)*qxhat + spsi*ctheta*qyhat;75 const realcmu = (-stheta*cpsi*cphi - spsi*sphi)*qxhat + ctheta*cpsi*qyhat;76 const realt = q*sqrt(req_minor*req_minor*cnu*cnu73 const double calpha = ctheta*cphi*qxhat + stheta*qyhat; 74 const double cnu = (-cphi*spsi*stheta + sphi*cpsi)*qxhat + spsi*ctheta*qyhat; 75 const double cmu = (-stheta*cpsi*cphi - spsi*sphi)*qxhat + ctheta*cpsi*qyhat; 76 const double t = q*sqrt(req_minor*req_minor*cnu*cnu 77 77 + req_major*req_major*cmu*cmu 78 78 + rpolar*rpolar*calpha*calpha); 79 79 SINCOS(t, st, ct); 80 const real fq = ( t==REAL(0.0) ? REAL(1.0) : REAL(3.0)*(st-t*ct)/(t*t*t) );81 const reals = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar);80 const double fq = ( t==0.0 ? 1.0 : 3.0*(st-t*ct)/(t*t*t) ); 81 const double s = (sld - solvent_sld) * form_volume(req_minor, req_major, rpolar); 82 82 83 return REAL(1.0e-4)* fq * fq * s * s;83 return 1.0e-4 * fq * fq * s * s; 84 84 } 85 85
Note: See TracChangeset
for help on using the changeset viewer.