Changeset bea7acb in sasmodels


Ignore:
Timestamp:
Mar 10, 2016 8:56:46 AM (4 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:
fc0fcd0, 3a45c2c
Parents:
0d0aee1 (diff), b15849c (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:
3 added
5 edited

Legend:

Unmodified
Added
Removed
  • explore/J1c.py

    re6f1410 r0a6da3c  
    88import pylab 
    99 
    10 mp.dps = 150 # number of digits to use in estimating true value 
    1110 
    1211SHOW_DIFF = True  # True if show diff rather than function value 
    1312LINEAR_X = False  # True if q is linearly spaced instead of log spaced 
    1413 
    15 def mp_J1c(vec): 
     14def mp_J1c(vec, bits=500): 
    1615    """ 
    1716    Direct calculation using sympy multiprecision library. 
    1817    """ 
    19     return [_mp_J1c(mp.mpf(x)) for x in vec] 
     18    with mp.workprec(bits): 
     19        return [_mp_J1c(mp.mpf(x)) for x in vec] 
    2020 
    2121def _mp_J1c(x): 
     
    2525    return mp.mpf(2)*mp.j1(x)/x 
    2626 
    27 def np_j1c(x, dtype): 
     27def np_J1c(x, dtype): 
    2828    """ 
    2929    Direct calculation using scipy. 
    3030    """ 
    31     from scipy.special import j1 
     31    from scipy.special import j1 as J1 
    3232    x = np.asarray(x, dtype) 
    33     return np.asarray(2, dtype)*j1(x)/x 
     33    return np.asarray(2, dtype)*J1(x)/x 
    3434 
    35 def cephes_j1c(x, dtype, n): 
     35def cephes_J1c(x, dtype, n): 
    3636    """ 
    3737    Calculation using pade approximant. 
    3838    """ 
     39    f = np.float64 if np.dtype(dtype) == np.float64 else np.float32 
    3940    x = np.asarray(x, dtype) 
    4041    ans = np.empty_like(x) 
     
    4243 
    4344    # Branch a 
    44     a_idx = ax < 8.0 
     45    a_idx = ax < f(8.0) 
    4546    a_xsq = x[a_idx]**2 
    4647    a_coeff1 = list(reversed((72362614232.0, -7895059235.0, 242396853.1, -2972611.439, 15704.48260, -30.16036606))) 
    4748    a_coeff2 = list(reversed((144725228442.0, 2300535178.0, 18583304.74, 99447.43394, 376.9991397, 1.0))) 
    48     a_ans1 = np.polyval(a_coeff1[n:], a_xsq) 
    49     a_ans2 = np.polyval(a_coeff2[n:], a_xsq) 
    50     ans[a_idx] = 2*a_ans1/a_ans2 
     49    a_ans1 = np.polyval(np.asarray(a_coeff1[n:], dtype), a_xsq) 
     50    a_ans2 = np.polyval(np.asarray(a_coeff2[n:], dtype), a_xsq) 
     51    ans[a_idx] = f(2.0)*a_ans1/a_ans2 
    5152 
    5253    # Branch b 
     
    5556    b_x = x[b_idx] 
    5657 
    57     b_y = 64.0/(b_ax**2) 
    58     b_xx = b_ax - 2.356194491 
     58    b_y = f(64.0)/(b_ax**2) 
     59    b_xx = b_ax - f(2.356194491) 
    5960    b_coeff1 = list(reversed((1.0, 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, -0.240337019e-6))) 
    6061    b_coeff2 = list(reversed((0.04687499995, -0.2002690873e-3, 0.8449199096e-5, -0.88228987e-6, 0.105787412e-6))) 
    61     b_ans1 = np.polyval(b_coeff1[n:], b_y) 
    62     b_ans2 = np.polyval(b_coeff2[n:], b_y) 
     62    b_ans1 = np.polyval(np.asarray(b_coeff1[n:], dtype),b_y) 
     63    b_ans2 = np.polyval(np.asarray(b_coeff2[n:], dtype), b_y) 
    6364    b_sn, b_cn = np.sin(b_xx), np.cos(b_xx) 
    64     ans[b_idx] = np.sign(b_x)*np.sqrt(0.636619772/b_ax) * (b_cn*b_ans1 - (8.0/b_ax)*b_sn*b_ans2)*2.0/b_x 
     65    ans[b_idx] = np.sign(b_x)*np.sqrt(f(0.636619772)/b_ax) * (b_cn*b_ans1 - (f(8.0)/b_ax)*b_sn*b_ans2)*f(2.0)/b_x 
    6566 
    6667    return ans 
     68 
     69def div_J1c(x, dtype): 
     70    f = np.float64 if np.dtype(dtype) == np.float64 else np.float32 
     71    x = np.asarray(x, dtype) 
     72    return f(2.0)*np.asarray([_J1(xi, f)/xi for xi in x], dtype) 
     73 
     74def _J1(x, f): 
     75    ax = abs(x) 
     76    if ax < f(8.0): 
     77        y = x*x 
     78        ans1 = x*(f(72362614232.0) 
     79                  + y*(f(-7895059235.0) 
     80                  + y*(f(242396853.1) 
     81                  + y*(f(-2972611.439) 
     82                  + y*(f(15704.48260) 
     83                  + y*(f(-30.16036606))))))) 
     84        ans2 = (f(144725228442.0) 
     85                  + y*(f(2300535178.0) 
     86                  + y*(f(18583304.74) 
     87                  + y*(f(99447.43394) 
     88                  + y*(f(376.9991397) 
     89                  + y))))) 
     90        return ans1/ans2 
     91    else: 
     92        y = f(64.0)/(ax*ax) 
     93        xx = ax - f(2.356194491) 
     94        ans1 = (f(1.0) 
     95                  + y*(f(0.183105e-2) 
     96                  + y*(f(-0.3516396496e-4) 
     97                  + y*(f(0.2457520174e-5) 
     98                  + y*f(-0.240337019e-6))))) 
     99        ans2 = (f(0.04687499995) 
     100                  + y*(f(-0.2002690873e-3) 
     101                  + y*(f(0.8449199096e-5) 
     102                  + y*(f(-0.88228987e-6) 
     103                  + y*f(0.105787412e-6))))) 
     104        sn, cn = np.sin(xx), np.cos(xx) 
     105        ans = np.sqrt(f(0.636619772)/ax) * (cn*ans1 - (f(8.0)/ax)*sn*ans2) 
     106        return -ans if (x < f(0.0)) else ans 
    67107 
    68108def plotdiff(x, target, actual, label): 
     
    84124    """ 
    85125    target = np.asarray(mp_J1c(x), 'double') 
    86     direct = np_j1c(x, precision) 
    87     approx0 = cephes_j1c(x, precision, 0) 
    88     approx1 = cephes_j1c(x, precision, 1) 
    89     plotdiff(x, target, direct, 'direct '+precision) 
    90     plotdiff(x, target, approx0, 'cephes '+precision) 
    91     #plotdiff(x, target, approx1, 'reduced '+precision) 
     126    #plotdiff(x, target, mp_J1c(x, 11), 'mp 11 bits') 
     127    plotdiff(x, target, np_J1c(x, precision), 'direct '+precision) 
     128    plotdiff(x, target, cephes_J1c(x, precision, 0), 'cephes '+precision) 
     129    #plotdiff(x, target, cephes_J1c(x, precision, 1), 'cephes '+precision) 
     130    #plotdiff(x, target, div_J1c(x, precision), 'cephes 2 J1(x)/x '+precision) 
    92131    pylab.xlabel("qr (1/Ang)") 
    93132    if SHOW_DIFF: 
     
    117156 
    118157if __name__ == "__main__": 
    119     print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3])) 
     158    #print "\n".join(str(x) for x in mp_J1c([1e-6,1e-5,1e-4,1e-3])) 
    120159    main() 
    121160    pylab.show() 
  • sasmodels/models/power_law.py

    r841753c rb15849c  
    11#power_law model 
    22#conversion of PowerLawAbsModel.py 
     3#converted by Steve King, Dec 2015 
    34 
    45r""" 
  • doc/rst_prolog

    r0a4628d r0d0aee1  
    5050.. |cm^-2| replace:: cm\ :sup:`-2` 
    5151.. |cm^3| replace:: cm\ :sup:`3` 
     52.. |1e15cm^3| replace:: 10\ :sup:`15`\ cm\ :sup:`3` 
    5253.. |cm^-3| replace:: cm\ :sup:`-3` 
    5354.. |sr^-1| replace:: sr\ :sup:`-1` 
  • sasmodels/convert.py

    r7d4b2ae r0d0aee1  
    5252    new model definition end with sld. 
    5353    """ 
    54     return dict((p, (v*1e6 if p.endswith('sld') else v)) 
     54    return dict((p, (v*1e6 if p.endswith('sld') else v*1e-15 if 'ndensity' in p else v)) 
    5555                for p, v in pars.items()) 
    5656 
     
    7070    new model definition end with sld. 
    7171    """ 
    72     return dict((p, (v*1e-6 if p.endswith('sld') else v)) 
     72    return dict((p, (v*1e-6 if p.endswith('sld') else v*1e15 if 'ndensity' in p else v)) 
    7373                for p, v in pars.items()) 
    7474 
  • sasmodels/generate.py

    rc437dbb r0d0aee1  
    242242    "1/Ang^2": "|Ang^-2|", 
    243243    "1e-6/Ang^2": "|1e-6Ang^-2|", 
     244    "1e15/cm^3": "|1e15cm^3|", 
    244245    "degrees": "degree", 
    245246    "1/cm": "|cm^-1|", 
     247    "1/cm^3": "|cm^-3|", 
    246248    "Ang/cm": "|Ang*cm^-1|", 
     249    "Ang^3": "|Ang^3|", 
    247250    "": "None", 
    248251    } 
Note: See TracChangeset for help on using the changeset viewer.