Changeset e4d8726 in sasmodels
- Timestamp:
- Feb 26, 2016 12:27:24 PM (9 years ago)
- 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. - Files:
-
- 29 added
- 29 edited
- 2 moved
Legend:
- Unmodified
- Added
- Removed
-
.gitignore
r66ebdd6 r4a82d4d 1 1 /build/ 2 2 /dist/ 3 /logs/ 3 4 /*.csv 4 5 *.pyc … … 16 17 /.project 17 18 /.pydevproject 19 /.idea 20 /sasmodels.egg-info/ -
doc/_extensions/dollarmath.py
r19dcb933 r103ea45 12 12 import re 13 13 14 _dollar = re.compile(r"(?:^|(?<=\s ))[$]([^\n]*?)(?<![\\])[$](?:$|(?=\s|[.,;\\]))")14 _dollar = re.compile(r"(?:^|(?<=\s|[(]))[$]([^\n]*?)(?<![\\])[$](?:$|(?=\s|[.,;)\\]))") 15 15 _notdollar = re.compile(r"\\[$]") 16 16 … … 51 51 assert replace_dollar(u"dollar\$ escape")==u"dollar$ escape" 52 52 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" 53 54 assert replace_dollar(u"emb\ $ed$\ ed")==u"emb\ :math:`ed`\ ed" 54 55 assert replace_dollar(u"$first$a")==u"$first$a" 55 56 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)" 56 61 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" 57 64 58 65 if __name__ == "__main__": -
doc/rst_prolog
r19dcb933 r0a4628d 60 60 .. |g/cm^3| replace:: g\ |cdot|\ cm\ :sup:`-3` 61 61 .. |fm^2| replace:: fm\ :sup:`2` 62 .. |Ang*cm^-1| replace:: |Ang|\ |cdot|\ cm\ :sup:`-1` -
multi_compare.sh
r5753e4e r0a4628d 1 #!/bin/ sh1 #!/bin/bash 2 2 3 3 sasview=( ../sasview/build/lib.* ) -
sasmodels/convert.py
r3964f92 r34d6cab 10 10 'broad_peak', 11 11 'two_lorentzian', 12 "two_power_law", 12 13 'gel_fit', 13 14 'gauss_lorentz_gel', … … 151 152 pars['string_thickness_pd_n'] = 0 152 153 pars['number_of_pearls_pd_n'] = 0 154 elif name == 'line': 155 pars['scale'] = 1 156 pars['background'] = 0 153 157 elif name == 'rpa': 154 158 pars['case_num'] = int(pars['case_num']) -
sasmodels/generate.py
reafc9fa rfa8011eb 191 191 The function :func:`make` loads the metadata from the module and returns 192 192 the kernel source. The function :func:`doc` extracts the doc string 193 and adds the parameter table to the top. The function :func:` sources`193 and adds the parameter table to the top. The function :func:`model_sources` 194 194 returns a list of files required by the model. 195 195 """ … … 206 206 import numpy as np 207 207 208 __all__ = ["make", "doc", "sources", "convert_type"]208 #__all__ = ["make", "doc", "model_sources", "convert_type"] 209 209 210 210 C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c') … … 233 233 "degrees": "degree", 234 234 "1/cm": "|cm^-1|", 235 "Ang/cm": "|Ang*cm^-1|", 235 236 "": "None", 236 237 } -
sasmodels/kernel_template.c
r840b859 r960cd80 12 12 #ifndef USE_OPENCL 13 13 # 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) 18 18 #include <limits> 19 19 #include <float.h> 20 20 #define kernel extern "C" __declspec( dllexport ) 21 21 inline double trunc(double x) { return x>=0?floor(x):-floor(-x); } 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 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 50 50 #else 51 51 #define kernel extern "C" … … 260 260 const double vol_weight = VOLUME_WEIGHT_PRODUCT; 261 261 vol += vol_weight*form_volume(VOLUME_PARAMETERS); 262 #endif263 262 norm_vol += vol_weight; 263 #endif 264 264 } 265 265 //else { printf("exclude qx,qy,I:%%g,%%g,%%g\n",qi,scattering); } -
sasmodels/models/be_polyelectrolyte.py
r168052c r0e86967 127 127 :return: 2D-Intensity 128 128 """ 129 i q= Iq(sqrt(qx**2 + qy**2), *args)130 return i q129 intensity = Iq(sqrt(qx**2 + qy**2), *args) 130 return intensity 131 131 132 132 Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values -
sasmodels/models/core_shell_bicelle.py
r8007311 rfa8011eb 7 7 The form factor is normalized by the particle volume. 8 8 9 .. _core-shell- cylinder-geometry:9 .. _core-shell-bicelle-geometry: 10 10 11 11 .. figure:: img/core_shell_bicelle_geometry.png -
sasmodels/models/core_shell_cylinder.py
reb69cce rf0aa7f8 155 155 156 156 def ER(radius, thickness, length): 157 """ 158 Returns the effective radius used in the S*P calculation 159 """ 157 160 radius = radius + thickness 158 161 length = length + 2 * thickness … … 161 164 162 165 def VR(radius, thickness, length): 166 """ 167 Returns volume ratio 168 """ 163 169 whole = pi * (radius + thickness) ** 2 * (length + 2 * thickness) 164 170 core = pi * radius ** 2 * length -
sasmodels/models/correlation_length.py
r5054e80 r3e6c5c1 78 78 # names and the target sasview model name. 79 79 oldname = 'CorrLengthModel' 80 # pylint: disable=bad-continuation 81 oldpars = dict( 82 lorentz_scale='scale_l', porod_scale='scale_p', 80 81 oldpars = dict(lorentz_scale='scale_l', porod_scale='scale_p', 83 82 cor_length='length_l', exponent_p='exponent_p', 84 exponent_l='exponent_l' 85 ) 83 exponent_l='exponent_l') 86 84 87 tests = [ 88 [{}, 0.001, 1009.98], 85 tests = [[{}, 0.001, 1009.98], 89 86 [{}, 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 144 144 145 145 def ER(radius, length): 146 """ 147 Return equivalent radius (ER) 148 """ 146 149 ddd = 0.75 * radius * (2 * radius * length + (length + radius) * (length + pi * radius)) 147 150 return 0.5 * (ddd) ** (1. / 3.) -
sasmodels/models/dab.py
reb69cce r94bd809 1 1 r""" 2 3 2 Calculates the scattering from a randomly distributed, two-phase system based on 4 3 the Debye-Anderson-Brumberger (DAB) model for such systems. The two-phase system -
sasmodels/models/ellipsoid.py
r9c461c7 r431caae 111 111 ---------- 112 112 113 L A Feigin and D I Svergun. *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum, 114 New York, 1987. 113 L A Feigin and D I Svergun. 114 *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, 115 Plenum Press, New York, 1987. 115 116 """ 116 117 -
sasmodels/models/fcc.py
r9aac25d rc0ccea8 114 114 single = False 115 115 116 # pylint: disable=bad-whitespace, line-too-long 116 117 # ["name", "units", default, [lower, upper], "type","description"], 117 118 parameters = [["dnn", "Ang", 220, [-inf, inf], "", "Nearest neighbour distance"], … … 124 125 ["psi", "degrees", 60, [-inf, inf], "orientation", "Out of plane angle"] 125 126 ] 127 # pylint: enable=bad-whitespace, line-too-long 126 128 127 129 source = ["lib/J1.c", "lib/gauss150.c", "lib/sphere_form.c", "fcc.c"] -
sasmodels/models/guinier.py
reb69cce r723cebe 36 36 37 37 # ["name", "units", default, [lower, upper], "type","description"], 38 parameters = [ 39 ["rg", "Ang", 60.0, [0, inf], "", "Radius of Gyration"], 40 ] 38 parameters = [["rg", "Ang", 60.0, [0, inf], "", "Radius of Gyration"]] 41 39 42 40 Iq = """ … … 51 49 52 50 # parameters for demo 53 demo = dict(scale=1.0, rg=60.0)51 demo = dict(scale=1.0, rg=60.0) 54 52 55 53 # For testing against the old sasview models, include the converted parameter … … 59 57 60 58 # parameters for unit tests 61 tests = [ 62 [{'rg' : 31.5}, 0.005, 0.991756] 63 ] 59 tests = [[{'rg' : 31.5}, 0.005, 0.991756]] -
sasmodels/models/hollow_cylinder.py
rec2ca99 re0fd913 93 93 source = ["lib/J1.c", "lib/gauss76.c", "hollow_cylinder.c"] 94 94 95 # pylint: disable=W0613 95 96 def ER(radius, core_radius, length): 96 97 """ -
sasmodels/models/lamellarCailleHG.py
r13ed84c r652a78a 146 146 solvent_sld='sld_solvent') 147 147 # 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 ] 148 tests = [[{'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 68 68 category = "shape:lamellae" 69 69 70 # pylint: disable=bad-whitespace, line-too-long 70 71 # ["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 ] 72 parameters = [["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 82 78 83 79 # No volume normalization despite having a volume parameter … … 113 109 114 110 demo = dict(scale=1, background=0, 115 tail_length=15, head_length=10,111 tail_length=15, head_length=10, 116 112 sld=0.4, head_sld=3.0, solvent_sld=6.0, 117 tail_length_pd= 118 head_length_pd= 113 tail_length_pd=0.2, tail_length_pd_n=40, 114 head_length_pd=0.01, head_length_pd_n=40) 119 115 120 116 oldname = 'LamellarFFHGModel' … … 122 118 sld='sld_tail', head_sld='sld_head', solvent_sld='sld_solvent') 123 119 # 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 120 tests = [[{'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 7 7 * 8 8 */ 9 #define N_POINTS_76 76 9 10 10 11 // Gaussians 11 constant double Gauss76Wt[ 76]={12 constant double Gauss76Wt[N_POINTS_76]={ 12 13 .00126779163408536, //0 13 14 .00294910295364247, … … 88 89 }; 89 90 90 constant double Gauss76Z[ 76]={91 constant double Gauss76Z[N_POINTS_76]={ 91 92 -.999505948362153, //0 92 93 -.997397786355355, -
sasmodels/models/lorentz.py
reb69cce reb5901b 56 56 57 57 # parameters for demo 58 demo = dict(scale=1.0, background=0.0,cor_length=50.0)58 demo = dict(scale=1.0, background=0.0, cor_length=50.0) 59 59 60 60 # For testing against the old sasview models, include the converted parameter … … 64 64 65 65 # parameters for unit tests 66 tests = [ 67 [{'cor_length' : 250},0.01,0.137931] 68 ] 66 tests = [[{'cor_length': 250}, 0.01, 0.137931]] -
sasmodels/models/peak_lorentz.py
r14ba6f6 r04b0b30 11 11 I(q) = \frac{scale}{\bigl(1+\bigl(\frac{q-q_0}{B}\bigr)^2\bigr)} + background 12 12 13 with the peak having height of $I_0$ centered at $q_0$ and having a HWHM (half-width half-maximum) of B. 13 with the peak having height of $I_0$ centered at $q_0$ and having 14 a HWHM (half-width half-maximum) of B. 14 15 15 16 For 2D data the scattering intensity is calculated in the same way as 1D, … … 55 56 56 57 def Iq(q, peak_pos, peak_hwhm): 58 """ 59 Return I(q) 60 """ 57 61 inten = (1/(1+((q-peak_pos)/peak_hwhm)**2)) 58 62 return inten … … 60 64 61 65 def Iqxy(qx, qy, *args): 66 """ 67 Return I(qx, qy) 68 """ 62 69 return Iq(sqrt(qx ** 2 + qy ** 2), *args) 63 70 Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values -
sasmodels/models/rpa.py
r82c299f rfa8011eb 43 43 component. 44 44 45 .. figure:: img/ image215.jpg45 .. figure:: img/rpa_1d.jpg 46 46 47 47 1D plot using the default values (w/500 data points). … … 89 89 # ["name", "units", default, [lower, upper], "type","description"], 90 90 parameters = [ 91 ["case_num", CASES, 0, [0, 10], "", "Component organization" 91 ["case_num", CASES, 0, [0, 10], "", "Component organization"], 92 92 93 93 ["Na", "", 1000.0, [1, inf], "", "Degree of polymerization"], -
sasmodels/models/sphere.py
r9c461c7 rc691551 100 100 101 101 def ER(radius): 102 """ 103 Return equivalent radius (ER) 104 """ 102 105 return radius 103 106 -
sasmodels/models/spherepy.py
reb69cce rd2950f4 118 118 g[low] = sqrt(1 - dlow2 / 4.) * (1 + dlow2 / 8.) + dlow2 / 2.*(1 - dlow2 / 16.) * log(dlow / (2. + sqrt(4. - dlow2))) 119 119 return g 120 sesans.vectorized = True # sesans accepts an darray of z values120 sesans.vectorized = True # sesans accepts an array of z values 121 121 122 122 def ER(radius): -
sasmodels/models/triaxial_ellipsoid.py
r9c461c7 r469e763 82 82 ---------- 83 83 84 L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and Neutron Scattering*, Plenum,85 New York, 1987.84 L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray 85 and Neutron Scattering*, Plenum, New York, 1987. 86 86 """ 87 87 … … 120 120 121 121 def ER(req_minor, req_major, rpolar): 122 """ 123 Returns the effective radius used in the S*P calculation 124 """ 122 125 import numpy as np 123 126 from .ellipsoid import ER as ellipsoid_ER -
sasmodels/models/vesicle.py
r216fa6d rfa8011eb 24 24 is a flat background level (due for example to incoherent scattering in the 25 25 case of neutrons), and $j_1$ is the spherical bessel function 26 $j_1 = ( sin(x) - xcos(x))/ x^2$.26 $j_1 = (\sin(x) - x \cos(x))/ x^2$. 27 27 28 28 The functional form is identical to a "typical" core-shell structure, except … … 35 35 thickness = $R_{\text{tot}} - R_{\text{core}}$. 36 36 37 .. figure: img/vesicle_geometry.jpg 37 .. figure:: img/vesicle_geometry.jpg 38 39 Vesicle geometry. 38 40 39 41 The 2D scattering intensity is the same as *P(q)* above, regardless of the … … 48 50 radius for *S(Q)* when *P(Q)* \* *S(Q)* is applied. 49 51 50 .. image:: img/vesicle_1d.jpg52 .. figure:: img/vesicle_1d.jpg 51 53 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. 55 57 56 58 REFERENCE … … 60 62 """ 61 63 62 import numpy as np63 64 from numpy import pi, inf 64 65 … … 131 132 # NOTE: test results taken from values returned by SasView 3.1.2 132 133 tests = [[{}, 0.0010005303255, 17139.8268799], 133 [{}, 0.200027832249, 0.130387268704 134 [{}, 0.200027832249, 0.130387268704], 134 135 [{}, 'ER', 130.], 135 136 [{}, 'VR', 0.54483386436], -
setup.py
r040575f r3eb3312 11 11 version = "1.0.0a", 12 12 description = "sasmodels package", 13 long_description=open('README. rst').read(),13 long_description=open('README.md').read(), 14 14 author = "SasView Collaboration", 15 15 author_email = "management@sasview.org", -
sasmodels/models/hardsphere.py
r093f754 r97e6d3c 69 69 return(HARDSPH); 70 70 } 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; 73 76 X=fabs(q*effect_radius*2.0); 74 77 … … 77 80 return(HARDSPH); 78 81 } 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; 82 86 G=0.5*volfraction*A; 83 87 84 88 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. 87 91 // 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 90 99 // combining the terms makes things worse at smallest Q in single precision 91 100 //FF = (8-0.8*X2)*A +(3.0-X2/3.)*2*B + (4+0.5*X2)*G +(A/35. +B/40. +G/50.)*X4; 92 101 // 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; 94 103 HARDSPH= 1./(1. + volfraction*FF ); 95 104 return(HARDSPH); 96 105 } 106 X4=X2*X2; 97 107 SINCOS(X,S,C); 98 108 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 ; 101 114 HARDSPH= 1./(1. + 24.*volfraction*FF/X2 ); 102 115 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 104 123 // 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 ); 105 124 // HARDSPH= 1./(1. + 24.*volfraction*FF/X2 );
Note: See TracChangeset
for help on using the changeset viewer.