Changeset ee60aa7 in sasmodels for sasmodels/models/mono_gauss_coil.c
- Timestamp:
- Sep 10, 2018 2:16:46 PM (6 years ago)
- Branches:
- master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- d299327
- Parents:
- 3f818b2
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/models/mono_gauss_coil.c
rd277229 ree60aa7 7 7 effective_radius(int mode, double rg) 8 8 { 9 if (mode == 1) { 9 switch (mode) { 10 case 1: // R_g 10 11 return rg; 11 } else if (mode == 2) {12 case 2: // 2R_g 12 13 return 2.0*rg; 13 } else if (mode == 3) {14 case 3: // 3R_g 14 15 return 3.0*rg; 15 } else {16 case 4: // (5/3)^0.5*R_g 16 17 return sqrt(5.0/3.0)*rg; 17 18 } 18 19 } 19 20 21 static double 22 gauss_coil(double qr) 23 { 24 const double x = qr*qr; 25 26 // Use series expansion at low q for higher accuracy. We could use 27 // smaller polynomials if we sacrifice some digits of precision or 28 // introduce an additional series expansion around x == 1. 29 // See explore/precision.py, gauss_coil function. 30 #if FLOAT_SIZE>4 // DOUBLE_PRECISION 31 // For double precision: use O(5) Pade with 0.5 cutoff (10 mad + 1 divide) 32 if (x < 0.5) { 33 // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 34 const double A1=1./12., A2=2./99., A3=1./2640., A4=1./23760., A5=-1./1995840.; 35 const double B1=5./12., B2=5./66., B3=1./132., B4=1./2376., B5=1./95040.; 36 return (((((A5*x + A4)*x + A3)*x + A2)*x + A1)*x + 1.) 37 / (((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.); 38 } 39 #else 40 // For single precision: use O(7) Taylor with 0.8 cutoff (7 mad) 41 if (x < 0.8) { 42 const double C0 = +1.; 43 const double C1 = -1./3.; 44 const double C2 = +1./12.; 45 const double C3 = -1./60.; 46 const double C4 = +1./360.; 47 const double C5 = -1./2520.; 48 const double C6 = +1./20160.; 49 const double C7 = -1./181440.; 50 return ((((((C7*x + C6)*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 51 } 52 #endif 53 54 return 2.0 * (expm1(-x) + x)/(x*x); 55 } 56 20 57 double Iq(double q, double i_zero, double rg) 21 58 { 22 const double uarg = square(q*rg); 23 const double inten; 24 if (q == 0) { 25 inten = i_zero; 26 } else { 27 inten = 2.0*i_zero * (exp(-uarg) + uarg - 1.0)/square(uarg); 28 } 29 30 return inten; 59 return i_zero * gauss_coil(q*rg); 31 60 }
Note: See TracChangeset
for help on using the changeset viewer.