Changeset 5a483877 in sasmodels


Ignore:
Timestamp:
Mar 8, 2016 8:59:33 AM (9 years ago)
Author:
piotr
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:
0d0aee1, 0a6da3c
Parents:
7d4b2ae (diff), 8dca856 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge remote-tracking branch 'origin/master'

Files:
10 added
5 edited
1 moved

Legend:

Unmodified
Added
Removed

  • sasmodels/models/lib/lanczos_gamma.c

    r5c1d341 re6f1410  
    22double lanczos_gamma(double q) 
    33{ 
    4     // Lanczos approximation to the Gamma function. 
     4    // Lanczos approximation to the Log Gamma function. 
    55 
    66    double x,y,tmp,ser; 

  • sasmodels/models/lib/sph_j1c.c

    rc437dbb re6f1410  
    55* Note that the values differ from sasview ~ 5e-12 rather than 5e-14, but 
    66* in this case it is likely cancellation errors in the original expression 
    7 * using double precision that are the source.  Single precision only 
    8 * requires the first 3 terms.  Double precision requires the 4th term. 
    9 * The fifth term is not needed, and is commented out. 
    10 * Taylor expansion: 
    11 *      1.0 + q2*(-3./30. + q2*(3./840.))+ q2*(-3./45360. + q2*(3./3991680.)))) 
     7* using double precision that are the source. 
    128*/ 
     9 
     10// The choice of the number of terms in the series and the cutoff value for 
     11// switching between series and direct calculation depends on the numeric 
     12// precision. 
     13// 
     14// Point where direct calculation reaches machine precision: 
     15// 
     16//   single machine precision eps 3e-8 at qr=1.1 ** 
     17//   double machine precision eps 4e-16 at qr=1.1 
     18// 
     19// Point where Taylor series reaches machine precision (eps), where taylor 
     20// series matches direct calculation (cross) and the error at that point: 
     21// 
     22//   prec   n eps  cross  error 
     23//   single 3 0.28  0.4   6.2e-7 
     24//   single 4 0.68  0.7   2.3e-7 
     25//   single 5 1.18  1.2   7.5e-8 
     26//   double 3 0.01  0.03  2.3e-13 
     27//   double 4 0.06  0.1   3.1e-14 
     28//   double 5 0.16  0.2   5.0e-15 
     29// 
     30// ** Note: relative error on single precision starts increase on the direct 
     31// method at qr=1.1, rising from 3e-8 to 5e-5 by qr=1e3.  This should be 
     32// safe for the sans range, with objects of 100 nm supported to a q of 0.1 
     33// while maintaining 5 digits of precision.  For usans/sesans, the objects 
     34// are larger but the q is smaller, so again it should be fine. 
     35// 
     36// See explore/sph_j1c.py for code to explore these ranges. 
     37 
     38// Use 4th order series 
     39#if FLOAT_SIZE>4 
     40#define SPH_J1C_CUTOFF 0.1 
     41#else 
     42#define SPH_J1C_CUTOFF 0.7 
     43#endif 
    1344 
    1445double sph_j1c(double q); 
    1546double sph_j1c(double q) 
    1647{ 
    17     const double q2 = q*q; 
    18     double sin_q, cos_q; 
    19  
    20     SINCOS(q, sin_q, cos_q); 
    21  
    22 #if FLOAT_SIZE>4 
    23 // For double precision, choose a cutoff of 0.18, which is the lower limit 
    24 // for the trig expression to 14 digits.  If only 12 digits are needed, then 
    25 // only 4 terms of the Taylor expansion are needed. 
    26 #define CUTOFF 0.18 
    27 #else 
    28 // For single precision, choose a cutoff halfway between the single precision 
    29 // lower limit for the trig expression of 1.03 and the upper limit of 1.3 
    30 // for the Taylor series expansion with 5 terms (or 9 if you count the zeros). 
    31 // For 5 digits of precision, you can drop two terms of the taylor series 
    32 // and choose a cutoff of 0.3. 
    33 #define CUTOFF 1.15 
    34 #endif 
    35  
    36     const double bessel = (q < CUTOFF) 
    37         ? (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.))))) 
    38         : 3.0*(sin_q/q - cos_q)/q2; 
    39     return bessel; 
    40  
    41 /* 
    42     // Code to test various expressions 
    43     // tested using sympy.mpmath.mp with 
    44     //    mp.dps = 50 (which is good for q down to 1e-6) 
    45     //    def j1c(x): return 3*(mp.sin(x)/x - mp.cos(x))/x**2 
    46     double ret; 
    47     if (sizeof(q2) > 8) { 
    48  
    49         ret = 3.0*(sinl(q)/q - cosl(q))/q2; 
    50 printf("%.15Lf %.15Lg\n", q, ret); 
    51     //} else if (q < 0.384038453352533) { 
    52     //} else if (q < 0.18) { 
    53     } else if (q < 18e-10) { 
    54         // NB: all are good to 5 digits below 0.18f 
    55         // good is defined as 14 digits for double and 7 for single 
    56         //ret = 1.0 + q2*q2*(3./840.) - q2*(3./30.); // good below 0.02d 0.34f 
    57         ret = square((1. + 3./5600.*q2*q2) - q2/20.); // good below 0.03d 0.08f 
    58         //ret = 1.0 - q2*(3./30.); // good below 0.02d 0.07f 
    59         //ret = 1.0 + q2*(-3./30. + q2*(3./840.)); // good below 0.02d 0.34f 
    60         //ret = 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))); // good below 0.1d 0.8f, 12 digits below 0.18d 
    61         //ret = 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.)))); // good below 0.18d 1.3f 
    62 printf("%.15g %.15g\n", q, ret); 
     48    if (q < SPH_J1C_CUTOFF) { 
     49        const double q2 = q*q; 
     50        return (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))));// + q2*(3./3991680.))))); 
    6351    } else { 
    64         // NB: can use a cutoff of 0.1f if the goal is 5 digits rather than 7 
    65         ret = 3.0*(sin_q/q - cos_q)/q2; // good above 0.18d, 1.03f 
    66 printf("%.15g %.15g\n", q, ret); 
     52        double sin_q, cos_q; 
     53        SINCOS(q, sin_q, cos_q); 
     54        return 3.0*(sin_q/q - cos_q)/(q*q); 
    6755    } 
    68     return ret; 
    69 */ 
    7056} 

  • sasmodels/convert.py

    r44bd2be r7d4b2ae  
    1515    'be_polyelectrolyte', 
    1616    'correlation_length', 
    17     'binary_hard_sphere' 
     17    'binary_hard_sphere', 
     18    'fractal_core_shell' 
    1819] 
    1920 

  • sasmodels/models/core_shell_sphere.c

    redc9f8d r7d4b2ae  
    55 
    66double Iq(double q, double radius, double thickness, double core_sld, double shell_sld, double solvent_sld) { 
    7     // Core first, then add in shell 
    8     const double core_qr = q * radius; 
    9     const double core_contrast = core_sld - shell_sld; 
    10     const double core_bes = sph_j1c(core_qr); 
    11     const double core_volume = 4.0 * M_PI / 3.0 * radius * radius * radius; 
    12     double f = core_volume * core_bes * core_contrast; 
    137 
    14     // Now the shell 
    15     const double shell_qr = q * (radius + thickness); 
    16     const double shell_contrast = shell_sld - solvent_sld; 
    17     const double shell_bes = sph_j1c(shell_qr); 
    18     const double shell_volume = 4.0 * M_PI / 3.0 * pow((radius + thickness), 3); 
    19     f += shell_volume * shell_bes * shell_contrast; 
    20     return f * f * 1.0e-4; 
     8 
     9    double intensity = core_shell_kernel(q, 
     10                              radius, 
     11                              thickness, 
     12                              core_sld, 
     13                              shell_sld, 
     14                              solvent_sld); 
     15    return intensity; 
    2116} 
    2217 

  • sasmodels/models/core_shell_sphere.py

    rfa8011eb r7d4b2ae  
    7373# pylint: enable=bad-whitespace, line-too-long 
    7474 
    75 source = ["lib/sph_j1c.c", "core_shell_sphere.c"] 
     75source = ["lib/sph_j1c.c", "lib/core_shell.c", "core_shell_sphere.c"] 
    7676 
    7777demo = dict(scale=1, background=0, radius=60, thickness=10, 
Note: See TracChangeset for help on using the changeset viewer.