Changeset e4d8726 in sasmodels


Ignore:
Timestamp:
Feb 26, 2016 12:27:24 PM (7 years ago)
Author:
richardh
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:
139c528
Parents:
97e6d3c (diff), 960cd80 (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 branch 'master' of https://github.com/SasView/sasmodels

Files:
29 added
29 edited
2 moved

Legend:

Unmodified
Added
Removed
  • .gitignore

    r66ebdd6 r4a82d4d  
    11/build/ 
    22/dist/ 
     3/logs/ 
    34/*.csv 
    45*.pyc 
     
    1617/.project 
    1718/.pydevproject 
     19/.idea 
     20/sasmodels.egg-info/ 
  • doc/_extensions/dollarmath.py

    r19dcb933 r103ea45  
    1212import re 
    1313 
    14 _dollar = re.compile(r"(?:^|(?<=\s))[$]([^\n]*?)(?<![\\])[$](?:$|(?=\s|[.,;\\]))") 
     14_dollar = re.compile(r"(?:^|(?<=\s|[(]))[$]([^\n]*?)(?<![\\])[$](?:$|(?=\s|[.,;)\\]))") 
    1515_notdollar = re.compile(r"\\[$]") 
    1616 
     
    5151    assert replace_dollar(u"dollar\$ escape")==u"dollar$ escape" 
    5252    assert replace_dollar(u"dollar \$escape\$ too")==u"dollar $escape$ too" 
     53    assert replace_dollar(u"spaces $in the$ math")==u"spaces :math:`in the` math" 
    5354    assert replace_dollar(u"emb\ $ed$\ ed")==u"emb\ :math:`ed`\ ed" 
    5455    assert replace_dollar(u"$first$a")==u"$first$a" 
    5556    assert replace_dollar(u"a$last$")==u"a$last$" 
     57    assert replace_dollar(u"$37")==u"$37" 
     58    assert replace_dollar(u"($37)")==u"($37)" 
     59    assert replace_dollar(u"$37 - $43")==u"$37 - $43" 
     60    assert replace_dollar(u"($37, $38)")==u"($37, $38)" 
    5661    assert replace_dollar(u"a $mid$dle a")==u"a $mid$dle a" 
     62    assert replace_dollar(u"a ($in parens$) a")==u"a (:math:`in parens`) a" 
     63    assert replace_dollar(u"a (again $in parens$) a")==u"a (again :math:`in parens`) a" 
    5764 
    5865if __name__ == "__main__": 
  • doc/rst_prolog

    r19dcb933 r0a4628d  
    6060.. |g/cm^3| replace:: g\ |cdot|\ cm\ :sup:`-3` 
    6161.. |fm^2| replace:: fm\ :sup:`2` 
     62.. |Ang*cm^-1| replace:: |Ang|\ |cdot|\ cm\ :sup:`-1` 
  • multi_compare.sh

    r5753e4e r0a4628d  
    1 #!/bin/sh 
     1#!/bin/bash 
    22 
    33sasview=( ../sasview/build/lib.* ) 
  • sasmodels/convert.py

    r3964f92 r34d6cab  
    1010    'broad_peak', 
    1111    'two_lorentzian', 
     12    "two_power_law", 
    1213    'gel_fit', 
    1314    'gauss_lorentz_gel', 
     
    151152        pars['string_thickness_pd_n'] = 0 
    152153        pars['number_of_pearls_pd_n'] = 0 
     154    elif name == 'line': 
     155        pars['scale'] = 1 
     156        pars['background'] = 0 
    153157    elif name == 'rpa': 
    154158        pars['case_num'] = int(pars['case_num']) 
  • sasmodels/generate.py

    reafc9fa rfa8011eb  
    191191The function :func:`make` loads the metadata from the module and returns 
    192192the kernel source.  The function :func:`doc` extracts the doc string 
    193 and adds the parameter table to the top.  The function :func:`sources` 
     193and adds the parameter table to the top.  The function :func:`model_sources` 
    194194returns a list of files required by the model. 
    195195""" 
     
    206206import numpy as np 
    207207 
    208 __all__ = ["make", "doc", "sources", "convert_type"] 
     208#__all__ = ["make", "doc", "model_sources", "convert_type"] 
    209209 
    210210C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c') 
     
    233233    "degrees": "degree", 
    234234    "1/cm": "|cm^-1|", 
     235    "Ang/cm": "|Ang*cm^-1|", 
    235236    "": "None", 
    236237    } 
  • sasmodels/kernel_template.c

    r840b859 r960cd80  
    1212#ifndef USE_OPENCL 
    1313#  ifdef __cplusplus 
    14      #include <cstdio> 
    15      #include <cmath> 
    16      using namespace std; 
    17      #if defined(_MSC_VER) 
     14      #include <cstdio> 
     15      #include <cmath> 
     16      using namespace std; 
     17      #if defined(_MSC_VER) 
    1818         #include <limits> 
    1919         #include <float.h> 
    2020         #define kernel extern "C" __declspec( dllexport ) 
    2121         inline double trunc(double x) { return x>=0?floor(x):-floor(-x); } 
    22              inline double fmin(double x, double y) { return x>y ? y : x; } 
    23              inline double fmax(double x, double y) { return x<y ? y : x; } 
    24              inline double isnan(double x) { return _isnan(x); } 
    25              #define NAN (std::numeric_limits<double>::quiet_NaN()) // non-signalling NaN 
    26              static double cephes_expm1(double x) { 
    27                  // Adapted from the cephes math library. 
    28                  // Copyright 1984 - 1992 by Stephen L. Moshier 
    29                  if (x != x || x == 0.0) { 
    30                      return x; // NaN and +/- 0 
    31                  } else if (x < -0.5 || x > 0.5) { 
    32                      return exp(x) - 1.0; 
    33                  } else { 
    34                      const double xsq = x*x; 
    35                      const double p = ((( 
    36                           +1.2617719307481059087798E-4)*xsq 
    37                       +3.0299440770744196129956E-2)*xsq 
    38                       +9.9999999999999999991025E-1); 
    39                  const double q = (((( 
    40                       +3.0019850513866445504159E-6)*xsq 
    41                       +2.5244834034968410419224E-3)*xsq 
    42                       +2.2726554820815502876593E-1)*xsq 
    43                       +2.0000000000000000000897E0); 
    44                  double r = x * p; 
    45                      r =  r / (q - r); 
    46                      return r+r; 
    47                  } 
    48              } 
    49              #define expm1 cephes_expm1 
     22         inline double fmin(double x, double y) { return x>y ? y : x; } 
     23         inline double fmax(double x, double y) { return x<y ? y : x; } 
     24         inline double isnan(double x) { return _isnan(x); } 
     25         #define NAN (std::numeric_limits<double>::quiet_NaN()) // non-signalling NaN 
     26         static double cephes_expm1(double x) { 
     27            // Adapted from the cephes math library. 
     28            // Copyright 1984 - 1992 by Stephen L. Moshier 
     29            if (x != x || x == 0.0) { 
     30               return x; // NaN and +/- 0 
     31            } else if (x < -0.5 || x > 0.5) { 
     32               return exp(x) - 1.0; 
     33            } else { 
     34               const double xsq = x*x; 
     35               const double p = ((( 
     36                  +1.2617719307481059087798E-4)*xsq 
     37                  +3.0299440770744196129956E-2)*xsq 
     38                  +9.9999999999999999991025E-1); 
     39               const double q = (((( 
     40                  +3.0019850513866445504159E-6)*xsq 
     41                  +2.5244834034968410419224E-3)*xsq 
     42                  +2.2726554820815502876593E-1)*xsq 
     43                  +2.0000000000000000000897E0); 
     44               double r = x * p; 
     45               r =  r / (q - r); 
     46               return r+r; 
     47             } 
     48         } 
     49         #define expm1 cephes_expm1 
    5050     #else 
    5151         #define kernel extern "C" 
     
    260260        const double vol_weight = VOLUME_WEIGHT_PRODUCT; 
    261261        vol += vol_weight*form_volume(VOLUME_PARAMETERS); 
    262       #endif 
    263262        norm_vol += vol_weight; 
     263      #endif 
    264264      } 
    265265      //else { printf("exclude qx,qy,I:%%g,%%g,%%g\n",qi,scattering); } 
  • sasmodels/models/be_polyelectrolyte.py

    r168052c r0e86967  
    127127    :return:     2D-Intensity 
    128128    """ 
    129     iq = Iq(sqrt(qx**2 + qy**2), *args) 
    130     return iq 
     129    intensity = Iq(sqrt(qx**2 + qy**2), *args) 
     130    return intensity 
    131131 
    132132Iqxy.vectorized = True  # Iqxy accepts an array of qx, qy values 
  • sasmodels/models/core_shell_bicelle.py

    r8007311 rfa8011eb  
    77The form factor is normalized by the particle volume. 
    88 
    9 .. _core-shell-cylinder-geometry: 
     9.. _core-shell-bicelle-geometry: 
    1010 
    1111.. figure:: img/core_shell_bicelle_geometry.png 
  • sasmodels/models/core_shell_cylinder.py

    reb69cce rf0aa7f8  
    155155 
    156156def ER(radius, thickness, length): 
     157    """ 
     158        Returns the effective radius used in the S*P calculation 
     159    """ 
    157160    radius = radius + thickness 
    158161    length = length + 2 * thickness 
     
    161164 
    162165def VR(radius, thickness, length): 
     166    """ 
     167        Returns volume ratio 
     168    """ 
    163169    whole = pi * (radius + thickness) ** 2 * (length + 2 * thickness) 
    164170    core = pi * radius ** 2 * length 
  • sasmodels/models/correlation_length.py

    r5054e80 r3e6c5c1  
    7878# names and the target sasview model name. 
    7979oldname = 'CorrLengthModel' 
    80 # pylint: disable=bad-continuation 
    81 oldpars = dict( 
    82                lorentz_scale='scale_l', porod_scale='scale_p', 
     80 
     81oldpars = dict(lorentz_scale='scale_l', porod_scale='scale_p', 
    8382               cor_length='length_l', exponent_p='exponent_p', 
    84                exponent_l='exponent_l' 
    85                ) 
     83               exponent_l='exponent_l') 
    8684 
    87 tests = [ 
    88          [{}, 0.001, 1009.98], 
     85tests = [[{}, 0.001, 1009.98], 
    8986         [{}, 0.150141, 0.174645], 
    90          [{}, 0.442528, 0.0203957] 
    91          ] 
    92 # pylint: enable=bad-continuation 
     87         [{}, 0.442528, 0.0203957]] 
  • sasmodels/models/cylinder.py

    reb69cce r0c3a226  
    144144 
    145145def ER(radius, length): 
     146    """ 
     147        Return equivalent radius (ER) 
     148    """ 
    146149    ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 
    147150    return 0.5 * (ddd) ** (1. / 3.) 
  • sasmodels/models/dab.py

    reb69cce r94bd809  
    11r""" 
    2  
    32Calculates the scattering from a randomly distributed, two-phase system based on 
    43the Debye-Anderson-Brumberger (DAB) model for such systems. The two-phase system 
  • sasmodels/models/ellipsoid.py

    r9c461c7 r431caae  
    111111---------- 
    112112 
    113 L A Feigin and D I Svergun. *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum, 
    114 New York, 1987. 
     113L A Feigin and D I Svergun. 
     114*Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, 
     115Plenum Press, New York, 1987. 
    115116""" 
    116117 
  • sasmodels/models/fcc.py

    r9aac25d rc0ccea8  
    114114single = False 
    115115 
     116# pylint: disable=bad-whitespace, line-too-long 
    116117#             ["name", "units", default, [lower, upper], "type","description"], 
    117118parameters = [["dnn", "Ang", 220, [-inf, inf], "", "Nearest neighbour distance"], 
     
    124125              ["psi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"] 
    125126             ] 
     127# pylint: enable=bad-whitespace, line-too-long 
    126128 
    127129source = ["lib/J1.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc.c"] 
  • sasmodels/models/guinier.py

    reb69cce r723cebe  
    3636 
    3737#             ["name", "units", default, [lower, upper], "type","description"], 
    38 parameters = [ 
    39               ["rg", "Ang", 60.0, [0, inf], "", "Radius of Gyration"], 
    40               ] 
     38parameters = [["rg", "Ang", 60.0, [0, inf], "", "Radius of Gyration"]] 
    4139 
    4240Iq = """ 
     
    5149 
    5250# parameters for demo 
    53 demo = dict(scale=1.0,rg=60.0) 
     51demo = dict(scale=1.0, rg=60.0) 
    5452 
    5553# For testing against the old sasview models, include the converted parameter 
     
    5957 
    6058# parameters for unit tests 
    61 tests = [ 
    62          [{'rg' : 31.5}, 0.005, 0.991756] 
    63          ] 
     59tests = [[{'rg' : 31.5}, 0.005, 0.991756]] 
  • sasmodels/models/hollow_cylinder.py

    rec2ca99 re0fd913  
    9393source = ["lib/J1.c", "lib/gauss76.c", "hollow_cylinder.c"] 
    9494 
     95# pylint: disable=W0613 
    9596def ER(radius, core_radius, length): 
    9697    """ 
  • sasmodels/models/lamellarCailleHG.py

    r13ed84c r652a78a  
    146146    solvent_sld='sld_solvent') 
    147147# 
    148 tests = [ 
    149         [ {'scale': 1.0, 'background' : 0.0, 'tail_length' : 10.0, 'head_length' : 2.0,'Nlayers' : 30.0, 'spacing' : 40., 
    150             'Caille_parameter' : 0.001, 'sld' : 0.4, 'head_sld' : 2.0, 'solvent_sld' : 6.0, 
    151             'tail_length_pd' : 0.0, 'head_length_pd' : 0.0, 'spacing_pd' : 0.0 }, [0.001], [6838238.571488]] 
    152         ] 
     148tests = [[{'scale': 1.0, 'background': 0.0, 'tail_length': 10.0, 'head_length': 2.0, 
     149           'Nlayers': 30.0, 'spacing': 40., 'Caille_parameter': 0.001, 'sld': 0.4, 
     150           'head_sld': 2.0, 'solvent_sld': 6.0, 'tail_length_pd': 0.0, 
     151           'head_length_pd': 0.0, 'spacing_pd': 0.0}, [0.001], [6838238.571488]]] 
  • sasmodels/models/lamellarFFHG.py

    r7f47777 r3eb6b90  
    6868category = "shape:lamellae" 
    6969 
     70# pylint: disable=bad-whitespace, line-too-long 
    7071#             ["name", "units", default, [lower, upper], "type","description"], 
    71 parameters = [["tail_length", "Ang",  15, [0, inf], "volume", 
    72                "Tail thickness"], 
    73               ["head_length", "Ang",  10, [0, inf], "volume", 
    74                "head thickness"], 
    75               ["sld", "1e-6/Ang^2", 0.4, [-inf,inf], "", 
    76                "Tail scattering length density"], 
    77               ["head_sld", "1e-6/Ang^2", 3.0, [-inf,inf], "", 
    78                "Head scattering length density"], 
    79               ["solvent_sld", "1e-6/Ang^2", 6, [-inf,inf], "", 
    80                "Solvent scattering length density"], 
    81              ] 
     72parameters = [["tail_length", "Ang",       15,   [0, inf],  "volume",  "Tail thickness"], 
     73              ["head_length", "Ang",       10,   [0, inf],  "volume",  "head thickness"], 
     74              ["sld",         "1e-6/Ang^2", 0.4, [-inf,inf], "",       "Tail scattering length density"], 
     75              ["head_sld",    "1e-6/Ang^2", 3.0, [-inf,inf], "",       "Head scattering length density"], 
     76              ["solvent_sld", "1e-6/Ang^2", 6,   [-inf,inf], "",       "Solvent scattering length density"]] 
     77# pylint: enable=bad-whitespace, line-too-long 
    8278 
    8379# No volume normalization despite having a volume parameter 
     
    113109 
    114110demo = dict(scale=1, background=0, 
    115             tail_length=15,head_length=10, 
     111            tail_length=15, head_length=10, 
    116112            sld=0.4, head_sld=3.0, solvent_sld=6.0, 
    117             tail_length_pd= 0.2, tail_length_pd_n=40, 
    118             head_length_pd= 0.01, head_length_pd_n=40) 
     113            tail_length_pd=0.2, tail_length_pd_n=40, 
     114            head_length_pd=0.01, head_length_pd_n=40) 
    119115 
    120116oldname = 'LamellarFFHGModel' 
     
    122118               sld='sld_tail', head_sld='sld_head', solvent_sld='sld_solvent') 
    123119# 
    124 tests = [ 
    125         [ {'scale': 1.0, 'background' : 0.0, 'tail_length' : 15.0, 'head_length' : 10.0,'sld' : 0.4, 
    126          'head_sld' : 3.0, 'solvent_sld' : 6.0, }, [0.001], [653143.9209]] 
    127         ] 
    128  
    129  
     120tests = [[{'scale': 1.0, 'background': 0.0, 'tail_length': 15.0, 'head_length': 10.0, 
     121           'sld': 0.4, 'head_sld': 3.0, 'solvent_sld': 6.0}, [0.001], [653143.9209]]] 
  • sasmodels/models/lib/gauss76.c

    r994d77f r66d119f  
    77 * 
    88 */ 
     9#define N_POINTS_76 76 
    910 
    1011// Gaussians 
    11 constant double Gauss76Wt[76]={ 
     12constant double Gauss76Wt[N_POINTS_76]={ 
    1213        .00126779163408536,             //0 
    1314        .00294910295364247, 
     
    8889}; 
    8990 
    90 constant double Gauss76Z[76]={ 
     91constant double Gauss76Z[N_POINTS_76]={ 
    9192        -.999505948362153,              //0 
    9293        -.997397786355355, 
  • sasmodels/models/lorentz.py

    reb69cce reb5901b  
    5656 
    5757# parameters for demo 
    58 demo = dict(scale=1.0,background=0.0,cor_length=50.0) 
     58demo = dict(scale=1.0, background=0.0, cor_length=50.0) 
    5959 
    6060# For testing against the old sasview models, include the converted parameter 
     
    6464 
    6565# parameters for unit tests 
    66 tests = [ 
    67          [{'cor_length' : 250},0.01,0.137931] 
    68          ] 
     66tests = [[{'cor_length': 250}, 0.01, 0.137931]] 
  • sasmodels/models/peak_lorentz.py

    r14ba6f6 r04b0b30  
    1111    I(q) = \frac{scale}{\bigl(1+\bigl(\frac{q-q_0}{B}\bigr)^2\bigr)} + background 
    1212 
    13 with the peak having height of $I_0$ centered at $q_0$ and having a HWHM (half-width half-maximum) of B. 
     13with the peak having height of $I_0$ centered at $q_0$ and having 
     14a HWHM (half-width half-maximum) of B. 
    1415 
    1516For 2D data the scattering intensity is calculated in the same way as 1D, 
     
    5556 
    5657def Iq(q, peak_pos, peak_hwhm): 
     58    """ 
     59        Return I(q) 
     60    """ 
    5761    inten = (1/(1+((q-peak_pos)/peak_hwhm)**2)) 
    5862    return inten 
     
    6064 
    6165def Iqxy(qx, qy, *args): 
     66    """ 
     67        Return I(qx, qy) 
     68    """ 
    6269    return Iq(sqrt(qx ** 2 + qy ** 2), *args) 
    6370Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values 
  • sasmodels/models/rpa.py

    r82c299f rfa8011eb  
    4343component. 
    4444 
    45 .. figure:: img/image215.jpg 
     45.. figure:: img/rpa_1d.jpg 
    4646 
    4747    1D plot using the default values (w/500 data points). 
     
    8989#   ["name", "units", default, [lower, upper], "type","description"], 
    9090parameters = [ 
    91     ["case_num", CASES, 0, [0, 10], "", "Component organization" ], 
     91    ["case_num", CASES, 0, [0, 10], "", "Component organization"], 
    9292 
    9393    ["Na", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
  • sasmodels/models/sphere.py

    r9c461c7 rc691551  
    100100 
    101101def ER(radius): 
     102    """ 
     103        Return equivalent radius (ER) 
     104    """ 
    102105    return radius 
    103106 
  • sasmodels/models/spherepy.py

    reb69cce rd2950f4  
    118118    g[low] = sqrt(1 - dlow2 / 4.) * (1 + dlow2 / 8.) + dlow2 / 2.*(1 - dlow2 / 16.) * log(dlow / (2. + sqrt(4. - dlow2))) 
    119119    return g 
    120 sesans.vectorized = True  # sesans accepts and array of z values 
     120sesans.vectorized = True  # sesans accepts an array of z values 
    121121 
    122122def ER(radius): 
  • sasmodels/models/triaxial_ellipsoid.py

    r9c461c7 r469e763  
    8282---------- 
    8383 
    84 L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum, 
    85 New York, 1987. 
     84L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray 
     85and Neutron Scattering*, Plenum, New York, 1987. 
    8686""" 
    8787 
     
    120120 
    121121def ER(req_minor, req_major, rpolar): 
     122    """ 
     123        Returns the effective radius used in the S*P calculation 
     124    """ 
    122125    import numpy as np 
    123126    from .ellipsoid import ER as ellipsoid_ER 
  • sasmodels/models/vesicle.py

    r216fa6d rfa8011eb  
    2424is a flat background level (due for example to incoherent scattering in the 
    2525case of neutrons), and $j_1$ is the spherical bessel function 
    26 $j_1 = (sin(x) - x cos(x))/ x^2$. 
     26$j_1 = (\sin(x) - x \cos(x))/ x^2$. 
    2727 
    2828The functional form is identical to a "typical" core-shell structure, except 
     
    3535thickness = $R_{\text{tot}} - R_{\text{core}}$. 
    3636 
    37 .. figure: img/vesicle_geometry.jpg 
     37.. figure:: img/vesicle_geometry.jpg 
     38 
     39    Vesicle geometry. 
    3840 
    3941The 2D scattering intensity is the same as *P(q)* above, regardless of the 
     
    4850radius for *S(Q)* when *P(Q)* \* *S(Q)* is applied. 
    4951 
    50 .. image:: img/vesicle_1d.jpg 
     52.. figure:: img/vesicle_1d.jpg 
    5153 
    52 *Figure. 1D plot using the default values given in the table 
    53 (w/200 data point). Polydispersity and instrumental resolution normally 
    54 will smear out most of the rapidly oscillating features.* 
     54    1D plot using the default values given in the table (w/200 data point). 
     55    Polydispersity and instrumental resolution normally will smear out most 
     56    of the rapidly oscillating features. 
    5557 
    5658REFERENCE 
     
    6062""" 
    6163 
    62 import numpy as np 
    6364from numpy import pi, inf 
    6465 
     
    131132# NOTE: test results taken from values returned by SasView 3.1.2 
    132133tests = [[{}, 0.0010005303255, 17139.8268799], 
    133          [{}, 0.200027832249, 0.130387268704 ], 
     134         [{}, 0.200027832249, 0.130387268704], 
    134135         [{}, 'ER', 130.], 
    135136         [{}, 'VR', 0.54483386436], 
  • setup.py

    r040575f r3eb3312  
    1111    version = "1.0.0a", 
    1212    description = "sasmodels package", 
    13     long_description=open('README.rst').read(), 
     13    long_description=open('README.md').read(), 
    1414    author = "SasView Collaboration", 
    1515    author_email = "management@sasview.org", 
  • sasmodels/models/hardsphere.py

    r093f754 r97e6d3c  
    6969               return(HARDSPH); 
    7070      } 
    71       D=pow((1.-volfraction),2); 
    72       A=pow((1.+2.*volfraction)/D, 2); 
     71      // removing use of pow(xxx,2) and rearranging the calcs of A, B & G cut ~40% off execution time ( 0.5 to 0.3 msec) 
     72      X = 1.0/( 1.0 -volfraction); 
     73      D= X*X; 
     74      A= (1.+2.*volfraction)*D; 
     75      A *=A; 
    7376      X=fabs(q*effect_radius*2.0); 
    7477 
     
    7780                 return(HARDSPH); 
    7881      } 
    79       X2=pow(X,2); 
    80       X4=pow(X2,2); 
    81       B=-6.*volfraction* pow((1.+0.5*volfraction)/D ,2); 
     82      X2 =X*X; 
     83      B = (1.0 +0.5*volfraction)*D; 
     84      B *= B; 
     85      B *= -6.*volfraction; 
    8286      G=0.5*volfraction*A; 
    8387 
    8488      if(X < 0.2) { 
    85       // use Taylor series expansion for small X, IT IS VERY PICKY ABOUT THE X CUT OFF VALUE, ought to be lower in double.  
    86       // No obvious way to rearrange the equations to avoid needing a very high number of significant figures.  
     89      // RKH Feb 2016, use Taylor series expansion for small X, IT IS VERY PICKY ABOUT THE X CUT OFF VALUE, ought to be lower in double.  
     90      // else no obvious way to rearrange the equations to avoid needing a very high number of significant figures.  
    8791      // Series expansion found using Mathematica software. Numerical test in .xls showed terms to X^2 are sufficient  
    88       // for 5 or 6 significant figures but I put the X^4 one in anyway  
    89             FF = 8*A +6*B + 4*G - (0.8*A +2.0*B/3.0 +0.5*G)*X2 +(A/35. +B/40. +G/50.)*X4; 
     92      // for 5 or 6 significant figures, but I put the X^4 one in anyway  
     93            //FF = 8*A +6*B + 4*G - (0.8*A +2.0*B/3.0 +0.5*G)*X2 +(A/35. +B/40. +G/50.)*X4; 
     94            // refactoring the polynomial makes it very slightly faster (0.5 not 0.6 msec) 
     95            //FF = 8*A +6*B + 4*G + ( -0.8*A -2.0*B/3.0 -0.5*G +(A/35. +B/40. +G/50.)*X2)*X2; 
     96 
     97            FF = 8.0*A +6.0*B + 4.0*G + ( -0.8*A -B/1.5 -0.5*G +(A/35. +0.0125*B +0.02*G)*X2)*X2; 
     98 
    9099            // combining the terms makes things worse at smallest Q in single precision 
    91100            //FF = (8-0.8*X2)*A +(3.0-X2/3.)*2*B + (4+0.5*X2)*G +(A/35. +B/40. +G/50.)*X4; 
    92101            // note that G = -volfraction*A/2, combining this makes no further difference at smallest Q 
    93             //FF = (8 +2.*volfraction + ( volfraction/4. -0.8 +(volfraction/100. -1./35.)*X2 )*X2 )*A  + (3.0 -X2/3. +X4/40)*2*B; 
     102            //FF = (8 +2.*volfraction + ( volfraction/4. -0.8 +(volfraction/100. -1./35.)*X2 )*X2 )*A  + (3.0 -X2/3. +X4/40.)*2.*B; 
    94103            HARDSPH= 1./(1. + volfraction*FF ); 
    95104            return(HARDSPH); 
    96105      } 
     106      X4=X2*X2; 
    97107      SINCOS(X,S,C); 
    98108 
    99 // RKH Feb 2016, use version from FISH code as it is better than original sasview one at small Q in single precision 
    100       FF=A*(S-X*C)/X + B*(2.*X*S -(X2-2.)*C -2.)/X2 + G*( (4.*X2*X -24.*X)*S -(X4 -12.*X2 +24.)*C +24. )/X4; 
     109// RKH Feb 2016, use version FISH code as is better than original sasview one at small Q in single precision, and more than twice as fast in double. 
     110      //FF=A*(S-X*C)/X + B*(2.*X*S -(X2-2.)*C -2.)/X2 + G*( (4.*X2*X -24.*X)*S -(X4 -12.*X2 +24.)*C +24. )/X4; 
     111      // refactoring the polynomial here & above makes it slightly faster 
     112 
     113      FF=  (( G*( (4.*X2 -24.)*X*S -(X4 -12.*X2 +24.)*C +24. )/X2 + B*(2.*X*S -(X2-2.)*C -2.) )/X + A*(S-X*C))/X ; 
    101114      HARDSPH= 1./(1. + 24.*volfraction*FF/X2 ); 
    102115 
    103 // rearrange the terms, is now about same as sasmodels 
     116      // changing /X and /X2 to *MX1 and *MX2, no significantg difference? 
     117      //MX=1.0/X; 
     118      //MX2=MX*MX; 
     119      //FF=  (( G*( (4.*X2 -24.)*X*S -(X4 -12.*X2 +24.)*C +24. )*MX2 + B*(2.*X*S -(X2-2.)*C -2.) )*MX + A*(S-X*C)) ; 
     120      //HARDSPH= 1./(1. + 24.*volfraction*FF*MX2*MX ); 
     121 
     122// grouping the terms, was about same as sasmodels for single precision issues 
    104123//     FF=A*(S/X-C) + B*(2.*S/X - C +2.0*(C-1.0)/X2) + G*( (4./X -24./X3)*S -(1.0 -12./X2 +24./X4)*C +24./X4 ); 
    105124//     HARDSPH= 1./(1. + 24.*volfraction*FF/X2 ); 
Note: See TracChangeset for help on using the changeset viewer.