Changeset 95ce773 in sasmodels


Ignore:
Timestamp:
Mar 18, 2016 12:45:38 PM (9 years ago)
Author:
wojciech
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:
9f7a852
Parents:
a629d8e
Message:

Bessel functions clean-up

Files:
6 edited

Legend:

Unmodified
Added
Removed

  • explore/J1.py

    raceac39 r95ce773  
    5050    model = core.load_model('bessel', dtype=dtype) 
    5151    calculator = direct_model.DirectModel(data.empty_data1D(x), model) 
    52     #ret = calculator(background=0) 
    53     #print ret 
     52 
    5453    return calculator(background=0) 
    5554 

  • explore/J1c.py

    rcbd37a7 r95ce773  
    22Show numerical precision of $2 J_1(x)/x$. 
    33""" 
     4import sys; sys.path.insert(0, '..') 
    45 
    56import numpy as np 
    6 from sympy.mpmath import mp 
     7try: 
     8    from mpmath import mp 
     9except: 
     10    # CRUFT: mpmath split out into its own package 
     11    from sympy.mpmath import mp 
    712#import matplotlib; matplotlib.use('TkAgg') 
    813import pylab 

  • sasmodels/models/bessel.py

    ra5af4e1 r95ce773  
    7878 
    7979Iq = """ 
    80     return J1(q); 
     80    return J1c(q); 
    8181    """ 
    8282 

  • sasmodels/models/lib/j0_cephes.c

    r094e320 r95ce773  
    5353 * except YP, YQ which are designed for absolute error. */ 
    5454 
    55 double j0( double ); 
    56 double j0(double x) { 
     55double J0(double x) { 
    5756 
    5857//Cephes single precission 
    … …  
    6059    double w, z, p, q, xn; 
    6160 
    62     const double TWOOPI = 6.36619772367581343075535E-1; 
     61    //const double TWOOPI = 6.36619772367581343075535E-1; 
    6362    const double SQ2OPI = 7.9788456080286535587989E-1; 
    6463    const double PIO4 = 7.85398163397448309616E-1; 
    … …  
    6766    const double DR2 = 3.04712623436620863991E1; 
    6867 
    69     const double PP[8] = { 
    70         7.96936729297347051624E-4, 
    71         8.28352392107440799803E-2, 
    72         1.23953371646414299388E0, 
    73         5.44725003058768775090E0, 
    74         8.74716500199817011941E0, 
    75         5.30324038235394892183E0, 
    76         9.99999999999999997821E-1, 
    77         0.0 
    78     }; 
    79  
    80     const double PQ[8] = { 
    81         9.24408810558863637013E-4, 
    82         8.56288474354474431428E-2, 
    83         1.25352743901058953537E0, 
    84         5.47097740330417105182E0, 
    85         8.76190883237069594232E0, 
    86         5.30605288235394617618E0, 
    87         1.00000000000000000218E0, 
    88         0.0 
    89     }; 
    90  
    91     const double QP[8] = { 
    92         -1.13663838898469149931E-2, 
    93         -1.28252718670509318512E0, 
    94         -1.95539544257735972385E1, 
    95         -9.32060152123768231369E1, 
    96         -1.77681167980488050595E2, 
    97         -1.47077505154951170175E2, 
    98         -5.14105326766599330220E1, 
    99         -6.05014350600728481186E0, 
    100     }; 
    101  
    102     const double QQ[8] = { 
    103         /*  1.00000000000000000000E0,*/ 
    104         6.43178256118178023184E1, 
    105         8.56430025976980587198E2, 
    106         3.88240183605401609683E3, 
    107         7.24046774195652478189E3, 
    108         5.93072701187316984827E3, 
    109         2.06209331660327847417E3, 
    110         2.42005740240291393179E2, 
    111     }; 
    112  
    113     const double YP[8] = { 
    114         1.55924367855235737965E4, 
    115         -1.46639295903971606143E7, 
    116         5.43526477051876500413E9, 
    117         -9.82136065717911466409E11, 
    118         8.75906394395366999549E13, 
    119         -3.46628303384729719441E15, 
    120         4.42733268572569800351E16, 
    121         -1.84950800436986690637E16, 
    122     }; 
    123  
    124     const double YQ[7] = { 
    125         /* 1.00000000000000000000E0,*/ 
    126         1.04128353664259848412E3, 
    127         6.26107330137134956842E5, 
    128         2.68919633393814121987E8, 
    129         8.64002487103935000337E10, 
    130         2.02979612750105546709E13, 
    131         3.17157752842975028269E15, 
    132         2.50596256172653059228E17, 
    133     }; 
    134  
    135     const double RP[8] = { 
    136         -4.79443220978201773821E9, 
    137         1.95617491946556577543E12, 
    138         -2.49248344360967716204E14, 
    139         9.70862251047306323952E15, 
    140         0.0, 
    141         0.0, 
    142         0.0, 
    143         0.0 
    144     }; 
    145  
    146     const double RQ[8] = { 
    147         /* 1.00000000000000000000E0,*/ 
    148         4.99563147152651017219E2, 
    149         1.73785401676374683123E5, 
    150         4.84409658339962045305E7, 
    151         1.11855537045356834862E10, 
    152         2.11277520115489217587E12, 
    153         3.10518229857422583814E14, 
    154         3.18121955943204943306E16, 
    155         1.71086294081043136091E18, 
    156     }; 
    15768 
    15869    if( x < 0 ) 
    … …  
    16576 
    16677            p = (z - DR1) * (z - DR2); 
    167             p = p * polevl( z, RP, 3)/p1evl( z, RQ, 8 ); 
     78            p = p * polevl( z, RPJ0, 3)/p1evl( z, RQJ0, 8 ); 
    16879            return( p ); 
    16980        } 
    … …  
    17182    w = 5.0/x; 
    17283    q = 25.0/(x*x); 
    173     p = polevl( q, PP, 6)/polevl( q, PQ, 6 ); 
    174     q = polevl( q, QP, 7)/p1evl( q, QQ, 7 ); 
     84    p = polevl( q, PPJ0, 6)/polevl( q, PQJ0, 6 ); 
     85    q = polevl( q, QPJ0, 7)/p1evl( q, QQJ0, 7 ); 
    17586    xn = x - PIO4; 
    17687 
    … …  
    18495    double xx, w, z, p, q, xn; 
    18596 
    186     const double YZ1 =  0.43221455686510834878; 
    187     const double YZ2 = 22.401876406482861405; 
    188     const double YZ3 = 64.130620282338755553; 
     97    //const double YZ1 =  0.43221455686510834878; 
     98    //const double YZ2 = 22.401876406482861405; 
     99    //const double YZ3 = 64.130620282338755553; 
    189100    const double DR1 =  5.78318596294678452118; 
    190101    const double PIO4F = 0.7853981633974483096; 
    191  
    192     double MO[8] = { 
    193         -6.838999669318810E-002, 
    194         1.864949361379502E-001, 
    195         -2.145007480346739E-001, 
    196         1.197549369473540E-001, 
    197         -3.560281861530129E-003, 
    198         -4.969382655296620E-002, 
    199         -3.355424622293709E-006, 
    200         7.978845717621440E-001 
    201     }; 
    202  
    203     double PH[8] = { 
    204         3.242077816988247E+001, 
    205         -3.630592630518434E+001, 
    206         1.756221482109099E+001, 
    207         -4.974978466280903E+000, 
    208         1.001973420681837E+000, 
    209         -1.939906941791308E-001, 
    210         6.490598792654666E-002, 
    211         -1.249992184872738E-001 
    212     }; 
    213  
    214     double YP[8] = { 
    215         9.454583683980369E-008, 
    216         -9.413212653797057E-006, 
    217         5.344486707214273E-004, 
    218         -1.584289289821316E-002, 
    219         1.707584643733568E-001, 
    220         0.0, 
    221         0.0, 
    222         0.0 
    223     }; 
    224  
    225     double JP[8] = { 
    226         -6.068350350393235E-008, 
    227         6.388945720783375E-006, 
    228         -3.969646342510940E-004, 
    229         1.332913422519003E-002, 
    230         -1.729150680240724E-001, 
    231         0.0, 
    232         0.0, 
    233         0.0 
    234     }; 
    235102 
    236103    if( x < 0 ) 
    … …  
    244111                    return( 1.0 - 0.25*z ); 
    245112 
    246             p = (z-DR1) * polevl( z, JP, 4); 
     113            p = (z-DR1) * polevl( z, JPJ0, 4); 
    247114            return( p ); 
    248115        } 
    … …  
    251118    w = sqrt(q); 
    252119 
    253     p = w * polevl( q, MO, 7); 
     120    p = w * polevl( q, MOJ0, 7); 
    254121    w = q*q; 
    255     xn = q * polevl( w, PH, 7) - PIO4F; 
     122    xn = q * polevl( w, PHJ0, 7) - PIO4F; 
    256123    p = p * cos(xn + xx); 
    257124    return(p); 

  • sasmodels/models/lib/j1_cephes.c

    rfad5dc1 r95ce773  
    4040*/ 
    4141 
    42 constant double RP[8] = { 
    43     -8.99971225705559398224E8, 
    44     4.52228297998194034323E11, 
    45     -7.27494245221818276015E13, 
    46     3.68295732863852883286E15, 
    47     0.0, 
    48     0.0, 
    49     0.0, 
    50     0.0 }; 
    51  
    52 constant double RQ[8] = { 
    53         6.20836478118054335476E2, 
    54         2.56987256757748830383E5, 
    55         8.35146791431949253037E7, 
    56         2.21511595479792499675E10, 
    57         4.74914122079991414898E12, 
    58         7.84369607876235854894E14, 
    59         8.95222336184627338078E16, 
    60         5.32278620332680085395E18 
    61     }; 
    62  
    63 constant double PP[8] = { 
    64     7.62125616208173112003E-4, 
    65     7.31397056940917570436E-2, 
    66     1.12719608129684925192E0, 
    67     5.11207951146807644818E0, 
    68     8.42404590141772420927E0, 
    69     5.21451598682361504063E0, 
    70     1.00000000000000000254E0, 
    71     0.0} ; 
    72  
    73  
    74 constant double PQ[8] = { 
    75     5.71323128072548699714E-4, 
    76     6.88455908754495404082E-2, 
    77     1.10514232634061696926E0, 
    78     5.07386386128601488557E0, 
    79     8.39985554327604159757E0, 
    80     5.20982848682361821619E0, 
    81     9.99999999999999997461E-1, 
    82     0.0 }; 
    83  
    84 constant double QP[8] = { 
    85     5.10862594750176621635E-2, 
    86     4.98213872951233449420E0, 
    87     7.58238284132545283818E1, 
    88     3.66779609360150777800E2, 
    89     7.10856304998926107277E2, 
    90     5.97489612400613639965E2, 
    91     2.11688757100572135698E2, 
    92     2.52070205858023719784E1 }; 
    93  
    94 constant double QQ[8] = { 
    95     7.42373277035675149943E1, 
    96     1.05644886038262816351E3, 
    97     4.98641058337653607651E3, 
    98     9.56231892404756170795E3, 
    99     7.99704160447350683650E3, 
    100     2.82619278517639096600E3, 
    101     3.36093607810698293419E2, 
    102     0.0 }; 
    103  
    104 constant double JP[8] = { 
    105         -4.878788132172128E-009, 
    106         6.009061827883699E-007, 
    107         -4.541343896997497E-005, 
    108         1.937383947804541E-003, 
    109         -3.405537384615824E-002, 
    110         0.0, 
    111         0.0, 
    112         0.0 
    113     }; 
    114  
    115 constant double MO1[8] = { 
    116         6.913942741265801E-002, 
    117         -2.284801500053359E-001, 
    118         3.138238455499697E-001, 
    119         -2.102302420403875E-001, 
    120         5.435364690523026E-003, 
    121         1.493389585089498E-001, 
    122         4.976029650847191E-006, 
    123         7.978845453073848E-001 
    124     }; 
    125  
    126 constant double PH1[8] = { 
    127         -4.497014141919556E+001, 
    128         5.073465654089319E+001, 
    129         -2.485774108720340E+001, 
    130         7.222973196770240E+000, 
    131         -1.544842782180211E+000, 
    132         3.503787691653334E-001, 
    133         -1.637986776941202E-001, 
    134         3.749989509080821E-001 
    135     }; 
    136  
    13742double J1(double x) { 
    13843 
    … …  
    15459        { 
    15560            z = x * x; 
    156             w = polevl( z, RP, 3 ) / p1evl( z, RQ, 8 ); 
     61            w = polevl( z, RPJ1, 3 ) / p1evl( z, RQJ1, 8 ); 
    15762            w = w * x * (z - Z1) * (z - Z2); 
    15863            return( w ); 
    … …  
    16267    z = w * w; 
    16368 
    164     p = polevl( z, PP, 6)/polevl( z, PQ, 6 ); 
    165     q = polevl( z, QP, 7)/p1evl( z, QQ, 7 ); 
     69    p = polevl( z, PPJ1, 6)/polevl( z, PQJ1, 6 ); 
     70    q = polevl( z, QPJ1, 7)/p1evl( z, QQJ1, 7 ); 
    16671 
    16772    xn = x - THPIO4; 
    … …  
    18994        { 
    19095            z = xx * xx; 
    191             p = (z-Z1) * xx * polevl( z, JP, 4 ); 
     96            p = (z-Z1) * xx * polevl( z, JPJ1, 4 ); 
    19297            return( p ); 
    19398        } 
    … …  
    196101    w = sqrt(q); 
    197102 
    198     p = w * polevl( q, MO1, 7); 
     103    p = w * polevl( q, MO1J1, 7); 
    199104    w = q*q; 
    200     xn = q * polevl( w, PH1, 7) - THPIO4F; 
     105    xn = q * polevl( w, PH1J1, 7) - THPIO4F; 
    201106    p = p * cos(xn + xx); 
    202107 
    … …  
    205110} 
    206111 
    207  
    208 /*double J1c(double x) { 
    209     return 
    210 }*/ 
     112//Finally J1c function that equals 2*J1(x)/x 
     113double J1c(double x) { 
     114    return (x != 0.0 ) ? 2.0*J1(x)/x : 1.0; 
     115} 

  • sasmodels/models/lib/polevl.c

    rfad5dc1 r95ce773  
    5151*/ 
    5252 
     53constant double RPJ1[8] = { 
     54    -8.99971225705559398224E8, 
     55    4.52228297998194034323E11, 
     56    -7.27494245221818276015E13, 
     57    3.68295732863852883286E15, 
     58    0.0, 
     59    0.0, 
     60    0.0, 
     61    0.0 }; 
     62 
     63constant double RQJ1[8] = { 
     64    6.20836478118054335476E2, 
     65    2.56987256757748830383E5, 
     66    8.35146791431949253037E7, 
     67    2.21511595479792499675E10, 
     68    4.74914122079991414898E12, 
     69    7.84369607876235854894E14, 
     70    8.95222336184627338078E16, 
     71    5.32278620332680085395E18 
     72    }; 
     73 
     74constant double PPJ1[8] = { 
     75    7.62125616208173112003E-4, 
     76    7.31397056940917570436E-2, 
     77    1.12719608129684925192E0, 
     78    5.11207951146807644818E0, 
     79    8.42404590141772420927E0, 
     80    5.21451598682361504063E0, 
     81    1.00000000000000000254E0, 
     82    0.0} ; 
     83 
     84 
     85constant double PQJ1[8] = { 
     86    5.71323128072548699714E-4, 
     87    6.88455908754495404082E-2, 
     88    1.10514232634061696926E0, 
     89    5.07386386128601488557E0, 
     90    8.39985554327604159757E0, 
     91    5.20982848682361821619E0, 
     92    9.99999999999999997461E-1, 
     93    0.0 }; 
     94 
     95constant double QPJ1[8] = { 
     96    5.10862594750176621635E-2, 
     97    4.98213872951233449420E0, 
     98    7.58238284132545283818E1, 
     99    3.66779609360150777800E2, 
     100    7.10856304998926107277E2, 
     101    5.97489612400613639965E2, 
     102    2.11688757100572135698E2, 
     103    2.52070205858023719784E1 }; 
     104 
     105constant double QQJ1[8] = { 
     106    7.42373277035675149943E1, 
     107    1.05644886038262816351E3, 
     108    4.98641058337653607651E3, 
     109    9.56231892404756170795E3, 
     110    7.99704160447350683650E3, 
     111    2.82619278517639096600E3, 
     112    3.36093607810698293419E2, 
     113    0.0 }; 
     114 
     115constant double JPJ1[8] = { 
     116    -4.878788132172128E-009, 
     117    6.009061827883699E-007, 
     118    -4.541343896997497E-005, 
     119    1.937383947804541E-003, 
     120    -3.405537384615824E-002, 
     121    0.0, 
     122    0.0, 
     123    0.0 
     124    }; 
     125 
     126constant double MO1J1[8] = { 
     127    6.913942741265801E-002, 
     128    -2.284801500053359E-001, 
     129    3.138238455499697E-001, 
     130    -2.102302420403875E-001, 
     131    5.435364690523026E-003, 
     132    1.493389585089498E-001, 
     133    4.976029650847191E-006, 
     134    7.978845453073848E-001 
     135    }; 
     136 
     137constant double PH1J1[8] = { 
     138    -4.497014141919556E+001, 
     139    5.073465654089319E+001, 
     140    -2.485774108720340E+001, 
     141    7.222973196770240E+000, 
     142    -1.544842782180211E+000, 
     143    3.503787691653334E-001, 
     144    -1.637986776941202E-001, 
     145    3.749989509080821E-001 
     146    }; 
     147 
     148 constant double PPJ0[8] = { 
     149        7.96936729297347051624E-4, 
     150        8.28352392107440799803E-2, 
     151        1.23953371646414299388E0, 
     152        5.44725003058768775090E0, 
     153        8.74716500199817011941E0, 
     154        5.30324038235394892183E0, 
     155        9.99999999999999997821E-1, 
     156        0.0 
     157    }; 
     158 
     159 constant double PQJ0[8] = { 
     160        9.24408810558863637013E-4, 
     161        8.56288474354474431428E-2, 
     162        1.25352743901058953537E0, 
     163        5.47097740330417105182E0, 
     164        8.76190883237069594232E0, 
     165        5.30605288235394617618E0, 
     166        1.00000000000000000218E0, 
     167        0.0 
     168    }; 
     169 
     170 constant double QPJ0[8] = { 
     171        -1.13663838898469149931E-2, 
     172        -1.28252718670509318512E0, 
     173        -1.95539544257735972385E1, 
     174        -9.32060152123768231369E1, 
     175        -1.77681167980488050595E2, 
     176        -1.47077505154951170175E2, 
     177        -5.14105326766599330220E1, 
     178        -6.05014350600728481186E0, 
     179    }; 
     180 
     181 constant double QQJ0[8] = { 
     182        /*  1.00000000000000000000E0,*/ 
     183        6.43178256118178023184E1, 
     184        8.56430025976980587198E2, 
     185        3.88240183605401609683E3, 
     186        7.24046774195652478189E3, 
     187        5.93072701187316984827E3, 
     188        2.06209331660327847417E3, 
     189        2.42005740240291393179E2, 
     190    }; 
     191 
     192 constant double YPJ0[8] = { 
     193        1.55924367855235737965E4, 
     194        -1.46639295903971606143E7, 
     195        5.43526477051876500413E9, 
     196        -9.82136065717911466409E11, 
     197        8.75906394395366999549E13, 
     198        -3.46628303384729719441E15, 
     199        4.42733268572569800351E16, 
     200        -1.84950800436986690637E16, 
     201 }; 
     202 
     203 
     204 constant double YQJ0[7] = { 
     205        /* 1.00000000000000000000E0,*/ 
     206        1.04128353664259848412E3, 
     207        6.26107330137134956842E5, 
     208        2.68919633393814121987E8, 
     209        8.64002487103935000337E10, 
     210        2.02979612750105546709E13, 
     211        3.17157752842975028269E15, 
     212        2.50596256172653059228E17, 
     213  }; 
     214 
     215 constant double RPJ0[8] = { 
     216        -4.79443220978201773821E9, 
     217        1.95617491946556577543E12, 
     218        -2.49248344360967716204E14, 
     219        9.70862251047306323952E15, 
     220        0.0, 
     221        0.0, 
     222        0.0, 
     223        0.0 
     224  }; 
     225 
     226 constant double RQJ0[8] = { 
     227        /* 1.00000000000000000000E0,*/ 
     228        4.99563147152651017219E2, 
     229        1.73785401676374683123E5, 
     230        4.84409658339962045305E7, 
     231        1.11855537045356834862E10, 
     232        2.11277520115489217587E12, 
     233        3.10518229857422583814E14, 
     234        3.18121955943204943306E16, 
     235        1.71086294081043136091E18, 
     236  }; 
     237 
     238 constant double MOJ0[8] = { 
     239        -6.838999669318810E-002, 
     240        1.864949361379502E-001, 
     241        -2.145007480346739E-001, 
     242        1.197549369473540E-001, 
     243        -3.560281861530129E-003, 
     244        -4.969382655296620E-002, 
     245        -3.355424622293709E-006, 
     246        7.978845717621440E-001 
     247  }; 
     248 
     249 constant double PHJ0[8] = { 
     250        3.242077816988247E+001, 
     251        -3.630592630518434E+001, 
     252        1.756221482109099E+001, 
     253        -4.974978466280903E+000, 
     254        1.001973420681837E+000, 
     255        -1.939906941791308E-001, 
     256        6.490598792654666E-002, 
     257        -1.249992184872738E-001 
     258  }; 
     259 
     260 constant double JPJ0[8] = { 
     261        -6.068350350393235E-008, 
     262        6.388945720783375E-006, 
     263        -3.969646342510940E-004, 
     264        1.332913422519003E-002, 
     265        -1.729150680240724E-001, 
     266        0.0, 
     267        0.0, 
     268        0.0 
     269 }; 
    53270 
    54271double polevl( double x, constant double *coef, int N ) { 
Note: See TracChangeset for help on using the changeset viewer.