Changeset 3a220e6 in sasmodels


Ignore:
Timestamp:
Oct 3, 2017 3:56:55 PM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
237c9cf
Parents:
6e72989
Message:

tweak cutoff point for debye taylor series expansion

Files:
2 edited

Legend:

Unmodified
Added
Removed
  • explore/precision.py

    r6e72989 r3a220e6  
    8181        elif xrange == "linear": 
    8282            lin_min, lin_max, lin_steps = 1, 1000, 2000 
     83            lin_min, lin_max, lin_steps = 0.001, 2, 2000 
    8384        elif xrange == "log": 
    8485            log_min, log_max, log_steps = -3, 5, 400 
     
    328329    ocl_function=make_ocl(""" 
    329330    const double qsq = q*q; 
    330     if (qsq < 0.1) { 
     331    if (qsq < -0.75*0.75) { // Pade approximation around 0 
     332        const double x = qsq; 
     333        if (0) { 
     334            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 4}] 
     335            return (x*x/180. + 1.)/((1./30.*x + 1./3.)*x + 1); 
     336        } else if (0) { 
     337            // padeapproximant[2*exp[-x^2] + x^2-1)/x^4, {x, 0, 6}] 
     338            // works for single precision, with qsq < 1 
     339            const double A1=1./24., A2=1./84, A3=-1./3360; 
     340            const double B1=3./8., B2=3./56., B3=1./336.; 
     341            return (((A3*x + A2)*x + A1)*x + 1.)/(((B3*x + B2)*x + B1)*x + 1.); 
     342        } else if (0) { 
     343            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 
     344            const double A1=1./15., A2=1./60, A3=0., A4=1./75600.; 
     345            const double B1=2./5., B2=1./15., B3=1./180., B4=1./5040.; 
     346            return ((((A4*x + A3)*x + A2)*x + A1)*x + 1.) 
     347                  /((((B4*x + B3)*x + B2)*x + B1)*x + 1.); 
     348 
     349        } else { 
     350            // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 
     351            // works for double precision, with qsq < 9/16 
     352            const double A1=1./12., A2=2./99., A3=1./2640., A4=1./23760., A5=-1./1995840.; 
     353            const double B1=5./12., B2=5./66., B3=1./132., B4=1./2376., B5=1./95040.; 
     354            return (((((A5*x + A4)*x + A3)*x + A2)*x + A1)*x + 1.) 
     355                  /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.); 
     356        } 
     357    } else if (qsq < 0.1) { // Taylor series around 0. 
    331358        const double x = qsq; 
    332359        const double C0 = +1.; 
     
    339366        const double C7 = -1./181440.; 
    340367        //return ((((C5*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 
    341         return (((((C6*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 
    342         //return ((((((C7*x + C6)*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 
    343     } /* else if (qsq < 1.0) { 
    344         // taylor series around q^2 = 0.5 
    345         const double x = qsq - 0.5; 
    346         const double sqrt_e = sqrt(M_E); 
    347         const double C0 = 8./sqrt_e - 4.; 
    348         const double C1 = 24. - 40./sqrt_e; 
    349         const double C2 = 132./sqrt_e - 80.; 
    350         const double C3 = 224. - 1108./(3.*sqrt_e); 
    351         const double C4 = 2849./(3.*sqrt_e) - 576.; 
    352         const double C5 = 1408 - 11607./(5.*sqrt_e); 
    353         return ((((C5*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 
    354     } */ else { 
     368        //return (((((C6*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 
     369        return ((((((C7*x + C6)*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 
     370    } else { 
    355371        return 2.*(exp(-qsq) + qsq - 1.)/(qsq*qsq); 
    356372    } 
  • sasmodels/models/lib/wrc_cyl.c

    r6e72989 r3a220e6  
    121121{ 
    122122#if FLOAT_SIZE>4 
    123 #define DEBYE_CUTOFF 0.01  // 1e-13 error 
     123#define DEBYE_CUTOFF 0.1  // 1e-14 error 
    124124#else 
    125 #define DEBYE_CUTOFF 1.0  // 1e-7 error 
     125#define DEBYE_CUTOFF 0.9  // 4e-7 error 
    126126#endif 
     127 
     128/* For double precision, the following gets 1e-15 error rather than 1e-14 
     129    if (qsq < 9./16.) { 
     130        // PadeApproximant[2*Exp[-x^2] + x^2-1)/x^4, {x, 0, 8}] 
     131        const double A1=1./12., A2=2./99., A3=1./2640., A4=1./23760., A5=-1./1995840.; 
     132        const double B1=5./12., B2=5./66., B3=1./132., B4=1./2376., B5=1./95040.; 
     133        const double x = qsq; 
     134        return (((((A5*x + A4)*x + A3)*x + A2)*x + A1)*x + 1.) 
     135                /(((((B5*x + B4)*x + B3)*x + B2)*x + B1)*x + 1.); 
     136    } 
     137*/ 
    127138 
    128139    if (qsq < DEBYE_CUTOFF) { 
     
    137148        const double C7 = -1./181440.; 
    138149        return ((((((C7*x + C6)*x + C5)*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 
    139     /* failed attempt to improve double precision better than 1e-13 
    140     } else if (qsq < 1.0) { 
    141         // taylor series around q^2 = 0.5 
    142         const double x = qsq - 0.5; 
    143         const double sqrt_e = sqrt(M_E); 
    144         const double C0 = 8./sqrt_e - 4.; 
    145         const double C1 = 24. - 40./sqrt_e; 
    146         const double C2 = 132./sqrt_e - 80.; 
    147         const double C3 = 224. - 1108./(3.*sqrt_e); 
    148         const double C4 = 2849./(3.*sqrt_e) - 576.; 
    149         const double C5 = 1408 - 11607./(5.*sqrt_e); 
    150         return ((((C5*x + C4)*x + C3)*x + C2)*x + C1)*x + C0; 
    151     */ 
    152150    } else { 
    153151        return 2.*(exp(-qsq) + qsq - 1.)/(qsq*qsq); 
Note: See TracChangeset for help on using the changeset viewer.