Changeset 7dde87f in sasmodels
- Timestamp:
- Dec 4, 2017 8:19:00 AM (7 years ago)
- Branches:
- master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- d0dc9a3
- Parents:
- a261a83 (diff), 0f113fb (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:
-
- 1 added
- 46 edited
Legend:
- Unmodified
- Added
- Removed
-
extra/pylint.rc
rb669b49 r6e68289 21 21 # List of plugins (as comma separated values of python modules names) to load, 22 22 # usually to register additional checkers. 23 load-plugins=pylint_numpy,pylint_pyopencl,pylint_sas23 #load-plugins=pylint_numpy,pylint_scipy,pylint_pyopencl,pylint_sas 24 24 25 25 # Use multiple processes to speed up Pylint. … … 280 280 # (useful for modules/projects where namespaces are manipulated during runtime 281 281 # and thus existing member attributes cannot be deduced by static analysis 282 ignored-modules=numpy,np,numpy.random, 283 bumps,sas, 282 ignored-modules=bumps,sas,numpy,numpy.random,scipy,scipy.special 284 283 285 284 # List of classes names for which member attributes should not be checked -
sasmodels/models/_spherepy.py
r2d81cfe ref07e95 40 40 John Wiley and Sons, New York, (1955) 41 41 42 * 2013/09/09 and 2014/01/06 - Description reviewed by S King and P Parker.*42 * **Last Reviewed by:** S King and P Parker **Date:** 2013/09/09 and 2014/01/06 43 43 """ 44 44 -
sasmodels/models/be_polyelectrolyte.py
r2d81cfe ref07e95 67 67 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 68 68 * **Last Modified by:** Paul Kienzle **Date:** July 24, 2016 69 * **Last Reviewed by:** Paul Butler and Richard Heenan **Date:** 70 October 07, 2016 69 * **Last Reviewed by:** Paul Butler and Richard Heenan **Date:** October 07, 2016 71 70 """ 72 71 -
sasmodels/models/fractal_core_shell.py
r2d81cfe ref07e95 55 55 56 56 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 57 * **Last Modified by:** Paul Butler and Paul Kienzle ** on:** November 27, 201658 * **Last Reviewed by:** Paul Butler and Paul Kienzle ** on:** November 27, 201657 * **Last Modified by:** Paul Butler and Paul Kienzle **Date:** November 27, 2016 58 * **Last Reviewed by:** Paul Butler and Paul Kienzle **Date:** November 27, 2016 59 59 """ 60 60 -
sasmodels/models/multilayer_vesicle.py
r2d81cfe ref07e95 107 107 * **Last Modified by:** Paul Kienzle **Date:** Feb 7, 2017 108 108 * **Last Reviewed by:** Paul Butler **Date:** March 12, 2017 109 110 109 """ 111 110 -
sasmodels/models/parallelepiped.py
r2d81cfe ref07e95 167 167 ---------------------------- 168 168 169 * **Author:** This model is based on form factor calculations implemented 170 in a c-library provided by the NIST Center for Neutron Research (Kline, 2006). 169 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 171 170 * **Last Modified by:** Paul Kienzle **Date:** April 05, 2017 172 171 * **Last Reviewed by:** Richard Heenan **Date:** April 06, 2017 -
sasmodels/models/polymer_micelle.py
r2d81cfe r791281c 28 28 .. math:: 29 29 P(q) &= N^2\beta^2_s\Phi(qr)^2 + N\beta^2_cP_c(q) 30 + 2N^2\beta_s\beta_cS_{sc} s_c(q) + N(N-1)\beta_c^2S_{cc}(q) \\30 + 2N^2\beta_s\beta_cS_{sc}(q) + N(N-1)\beta_c^2S_{cc}(q) \\ 31 31 \beta_s &= V_\text{core}(\rho_\text{core} - \rho_\text{solvent}) \\ 32 32 \beta_c &= V_\text{corona}(\rho_\text{corona} - \rho_\text{solvent}) … … 69 69 70 70 J Pedersen, *J. Appl. Cryst.*, 33 (2000) 637-640 71 72 * **Modified by:** Richard Heenan **Date:** March 20, 2016 73 * **Verified by:** Paul Kienzle **Date:** November 29, 2017 74 * **Description modified by:** Paul Kienzle **Date:** November 29, 2017 75 * **Description reviewed by:** Steve King **Date:** November 30, 2017 71 76 """ 72 77 -
sasmodels/models/pringle.py
r2d81cfe ref07e95 43 43 Derivation by Stefan Alexandru Rautu. 44 44 45 **Author:** Andrew Jackson **on:** 2008 46 47 **Last Modified by:** Wojciech Wpotrzebowski **on:** March 20, 2016 48 49 **Last Reviewed by:** Andrew Jackson **on:** September 26, 2016 45 * **Author:** Andrew Jackson **Date:** 2008 46 * **Last Modified by:** Wojciech Wpotrzebowski **Date:** March 20, 2016 47 * **Last Reviewed by:** Andrew Jackson **Date:** September 26, 2016 50 48 """ 51 49 -
sasmodels/models/raspberry.py
r2d81cfe ref07e95 102 102 Science*, 343(1) (2010) 36-41 103 103 104 **Author:** Andrew Jackson **on:** 2008 105 106 **Modified by:** Andrew Jackson **on:** March 20, 2016 107 108 **Reviewed by:** Andrew Jackson **on:** March 20, 2016 104 * **Author:** Andrew Jackson **Date:** 2008 105 * **Modified by:** Andrew Jackson **Date:** March 20, 2016 106 * **Reviewed by:** Andrew Jackson **Date:** March 20, 2016 109 107 """ 110 108 -
sasmodels/models/sphere.py
r2d81cfe ref07e95 40 40 John Wiley and Sons, New York, (1955) 41 41 42 * 2013/09/09 and 2014/01/06 - Description reviewed by S King and P Parker.*42 * **Last Reviewed by:** S King and P Parker **Date:** 2013/09/09 and 2014/01/06 43 43 """ 44 44 -
sasmodels/models/spinodal.py
r2d81cfe ref07e95 26 26 27 27 * **Author:** Dirk Honecker **Date:** Oct 7, 2016 28 * **Last Modified by:**29 * **Last Reviewed by:**30 28 """ 31 29 -
sasmodels/models/stacked_disks.py
r2d81cfe ref07e95 106 106 107 107 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 108 * **Last Modified by:** Paul Butler and Paul Kienzle ** on:** November 26, 2016109 * **Last Reviewed by:** Paul Butler and Paul Kienzle ** on:** November 26, 2016108 * **Last Modified by:** Paul Butler and Paul Kienzle **Date:** November 26, 2016 109 * **Last Reviewed by:** Paul Butler and Paul Kienzle **Date:** November 26, 2016 110 110 """ 111 111 -
sasmodels/models/triaxial_ellipsoid.py
r2d81cfe r37f08d2 162 162 Returns the effective radius used in the S*P calculation 163 163 """ 164 import numpy as np165 164 from .ellipsoid import ER as ellipsoid_ER 166 165 -
sasmodels/models/two_lorentzian.py
r2d81cfe ref07e95 27 27 None. 28 28 29 **Author:** NIST IGOR/DANSE **on:** pre 2010 30 31 **Last Modified by:** Piotr rozyczko **on:** January 29, 2016 32 33 **Last Reviewed by:** Paul Butler **on:** March 21, 2016 29 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 30 * **Last Modified by:** Piotr rozyczko **Date:** January 29, 2016 31 * **Last Reviewed by:** Paul Butler **Date:** March 21, 2016 34 32 """ 35 33 -
sasmodels/models/two_power_law.py
r2d81cfe ref07e95 37 37 None. 38 38 39 **Author:** NIST IGOR/DANSE **on:** pre 2010 40 41 **Last Modified by:** Wojciech Wpotrzebowski **on:** February 18, 2016 42 43 **Last Reviewed by:** Paul Butler **on:** March 21, 2016 39 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 40 * **Last Modified by:** Wojciech Wpotrzebowski **Date:** February 18, 2016 41 * **Last Reviewed by:** Paul Butler **Date:** March 21, 2016 44 42 """ 45 43 -
sasmodels/models/vesicle.py
r2d81cfe ref07e95 56 56 Sons, New York, (1955) 57 57 58 **Author:** NIST IGOR/DANSE **on:** pre 2010 59 60 **Last Modified by:** Paul Butler **on:** March 20, 2016 61 62 **Last Reviewed by:** Paul Butler **on:** March 20, 2016 58 * **Author:** NIST IGOR/DANSE **Date:** pre 2010 59 * **Last Modified by:** Paul Butler **Date:** March 20, 2016 60 * **Last Reviewed by:** Paul Butler **Date:** March 20, 2016 63 61 """ 64 62 -
.gitignore
rc26897a r9248bf7 22 22 /example/Fit_*/ 23 23 /example/batch_fit.csv 24 /sasmodels/models/lib/gauss*.c -
sasmodels/compare.py
r0e55afe ra261a83 42 42 from .data import plot_theory, empty_data1D, empty_data2D, load_data 43 43 from .direct_model import DirectModel, get_mesh 44 from .generate import FLOAT_RE 44 from .generate import FLOAT_RE, set_integration_size 45 45 from .weights import plot_weights 46 46 … … 706 706 return data, index 707 707 708 def make_engine(model_info, data, dtype, cutoff ):708 def make_engine(model_info, data, dtype, cutoff, ngauss=0): 709 709 # type: (ModelInfo, Data, str, float) -> Calculator 710 710 """ … … 714 714 than OpenCL. 715 715 """ 716 if ngauss: 717 set_integration_size(model_info, ngauss) 718 716 719 if dtype is None or not dtype.endswith('!'): 717 720 return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff) … … 954 957 'poly', 'mono', 'cutoff=', 955 958 'magnetic', 'nonmagnetic', 956 'accuracy=', 959 'accuracy=', 'ngauss=', 957 960 'neval=', # for timing... 958 961 … … 1089 1092 'show_weights' : False, 1090 1093 'sphere' : 0, 1094 'ngauss' : '0', 1091 1095 } 1092 1096 for arg in flags: … … 1115 1119 elif arg.startswith('-engine='): opts['engine'] = arg[8:] 1116 1120 elif arg.startswith('-neval='): opts['count'] = arg[7:] 1121 elif arg.startswith('-ngauss='): opts['ngauss'] = arg[8:] 1117 1122 elif arg.startswith('-random='): 1118 1123 opts['seed'] = int(arg[8:]) … … 1169 1174 1170 1175 comparison = any(PAR_SPLIT in v for v in values) 1176 1171 1177 if PAR_SPLIT in name: 1172 1178 names = name.split(PAR_SPLIT, 2) … … 1181 1187 return None 1182 1188 1189 if PAR_SPLIT in opts['ngauss']: 1190 opts['ngauss'] = [int(k) for k in opts['ngauss'].split(PAR_SPLIT, 2)] 1191 comparison = True 1192 else: 1193 opts['ngauss'] = [int(opts['ngauss'])]*2 1194 1183 1195 if PAR_SPLIT in opts['engine']: 1184 1196 opts['engine'] = opts['engine'].split(PAR_SPLIT, 2) … … 1199 1211 opts['cutoff'] = [float(opts['cutoff'])]*2 1200 1212 1201 base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0]) 1213 base = make_engine(model_info[0], data, opts['engine'][0], 1214 opts['cutoff'][0], opts['ngauss'][0]) 1202 1215 if comparison: 1203 comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1]) 1216 comp = make_engine(model_info[1], data, opts['engine'][1], 1217 opts['cutoff'][1], opts['ngauss'][1]) 1204 1218 else: 1205 1219 comp = None -
sasmodels/generate.py
r2d81cfe ra261a83 270 270 """ 271 271 272 273 def set_integration_size(info, n): 274 # type: (ModelInfo, int) -> None 275 """ 276 Update the model definition, replacing the gaussian integration with 277 a gaussian integration of a different size. 278 279 Note: this really ought to be a method in modelinfo, but that leads to 280 import loops. 281 """ 282 if (info.source and any(lib.startswith('lib/gauss') for lib in info.source)): 283 import os.path 284 from .gengauss import gengauss 285 path = os.path.join(MODEL_PATH, "lib", "gauss%d.c"%n) 286 if not os.path.exists(path): 287 gengauss(n, path) 288 info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss') 289 else lib for lib in info.source] 272 290 273 291 def format_units(units): -
sasmodels/models/barbell.c
rbecded3 r74768cb 23 23 const double qab_r = radius_bell*qab; // Q*R*sin(theta) 24 24 double total = 0.0; 25 for (int i = 0; i < 76; i++){26 const double t = G auss76Z[i]*zm + zb;25 for (int i = 0; i < GAUSS_N; i++){ 26 const double t = GAUSS_Z[i]*zm + zb; 27 27 const double radical = 1.0 - t*t; 28 28 const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 29 29 const double Fq = cos(m*t + b) * radical * bj; 30 total += G auss76Wt[i] * Fq;30 total += GAUSS_W[i] * Fq; 31 31 } 32 32 // translate dx in [-1,1] to dx in [lower,upper] … … 73 73 const double zb = M_PI_4; 74 74 double total = 0.0; 75 for (int i = 0; i < 76; i++){76 const double alpha= G auss76Z[i]*zm + zb;75 for (int i = 0; i < GAUSS_N; i++){ 76 const double alpha= GAUSS_Z[i]*zm + zb; 77 77 double sin_alpha, cos_alpha; // slots to hold sincos function output 78 78 SINCOS(alpha, sin_alpha, cos_alpha); 79 79 const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 80 total += G auss76Wt[i] * Aq * Aq * sin_alpha;80 total += GAUSS_W[i] * Aq * Aq * sin_alpha; 81 81 } 82 82 // translate dx in [-1,1] to dx in [lower,upper] -
sasmodels/models/bcc_paracrystal.c
rea60e08 r74768cb 81 81 82 82 double outer_sum = 0.0; 83 for(int i=0; i< 150; i++) {83 for(int i=0; i<GAUSS_N; i++) { 84 84 double inner_sum = 0.0; 85 const double theta = G auss150Z[i]*theta_m + theta_b;85 const double theta = GAUSS_Z[i]*theta_m + theta_b; 86 86 double sin_theta, cos_theta; 87 87 SINCOS(theta, sin_theta, cos_theta); 88 88 const double qc = q*cos_theta; 89 89 const double qab = q*sin_theta; 90 for(int j=0;j< 150;j++) {91 const double phi = G auss150Z[j]*phi_m + phi_b;90 for(int j=0;j<GAUSS_N;j++) { 91 const double phi = GAUSS_Z[j]*phi_m + phi_b; 92 92 double sin_phi, cos_phi; 93 93 SINCOS(phi, sin_phi, cos_phi); … … 95 95 const double qb = qab*sin_phi; 96 96 const double form = bcc_Zq(qa, qb, qc, dnn, d_factor); 97 inner_sum += G auss150Wt[j] * form;97 inner_sum += GAUSS_W[j] * form; 98 98 } 99 99 inner_sum *= phi_m; // sum(f(x)dx) = sum(f(x)) dx 100 outer_sum += G auss150Wt[i] * inner_sum * sin_theta;100 outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 101 101 } 102 102 outer_sum *= theta_m; -
sasmodels/models/capped_cylinder.c
rbecded3 r74768cb 30 30 const double qab_r = radius_cap*qab; // Q*R*sin(theta) 31 31 double total = 0.0; 32 for (int i=0; i< 76 ;i++) {33 const double t = G auss76Z[i]*zm + zb;32 for (int i=0; i<GAUSS_N; i++) { 33 const double t = GAUSS_Z[i]*zm + zb; 34 34 const double radical = 1.0 - t*t; 35 35 const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 36 36 const double Fq = cos(m*t + b) * radical * bj; 37 total += G auss76Wt[i] * Fq;37 total += GAUSS_W[i] * Fq; 38 38 } 39 39 // translate dx in [-1,1] to dx in [lower,upper] … … 95 95 const double zb = M_PI_4; 96 96 double total = 0.0; 97 for (int i=0; i< 76;i++) {98 const double theta = G auss76Z[i]*zm + zb;97 for (int i=0; i<GAUSS_N ;i++) { 98 const double theta = GAUSS_Z[i]*zm + zb; 99 99 double sin_theta, cos_theta; // slots to hold sincos function output 100 100 SINCOS(theta, sin_theta, cos_theta); … … 103 103 const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 104 104 // scale by sin_theta for spherical coord integration 105 total += G auss76Wt[i] * Aq * Aq * sin_theta;105 total += GAUSS_W[i] * Aq * Aq * sin_theta; 106 106 } 107 107 // translate dx in [-1,1] to dx in [lower,upper] -
sasmodels/models/core_shell_bicelle.c
rbecded3 r74768cb 52 52 53 53 double total = 0.0; 54 for(int i=0;i< N_POINTS_76;i++) {55 double theta = (G auss76Z[i] + 1.0)*uplim;54 for(int i=0;i<GAUSS_N;i++) { 55 double theta = (GAUSS_Z[i] + 1.0)*uplim; 56 56 double sin_theta, cos_theta; // slots to hold sincos function output 57 57 SINCOS(theta, sin_theta, cos_theta); 58 58 double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 59 59 halflength, sld_core, sld_face, sld_rim, sld_solvent); 60 total += G auss76Wt[i]*fq*fq*sin_theta;60 total += GAUSS_W[i]*fq*fq*sin_theta; 61 61 } 62 62 -
sasmodels/models/core_shell_bicelle_elliptical.c
r82592da rd4db147 37 37 //initialize integral 38 38 double outer_sum = 0.0; 39 for(int i=0;i< 76;i++) {39 for(int i=0;i<GAUSS_N;i++) { 40 40 //setup inner integral over the ellipsoidal cross-section 41 41 //const double va = 0.0; 42 42 //const double vb = 1.0; 43 //const double cos_theta = ( G auss76Z[i]*(vb-va) + va + vb )/2.0;44 const double cos_theta = ( G auss76Z[i] + 1.0 )/2.0;43 //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 44 const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0; 45 45 const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 46 46 const double qab = q*sin_theta; … … 49 49 const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 50 50 double inner_sum=0.0; 51 for(int j=0;j< 76;j++) {51 for(int j=0;j<GAUSS_N;j++) { 52 52 //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 53 53 // inner integral limits 54 54 //const double vaj=0.0; 55 55 //const double vbj=M_PI; 56 //const double phi = ( G auss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;57 const double phi = ( G auss76Z[j] +1.0)*M_PI_2;56 //const double phi = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 57 const double phi = ( GAUSS_Z[j] +1.0)*M_PI_2; 58 58 const double rr = sqrt(r2A - r2B*cos(phi)); 59 59 const double be1 = sas_2J1x_x(rr*qab); … … 61 61 const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 62 62 63 inner_sum += G auss76Wt[j] * fq * fq;63 inner_sum += GAUSS_W[j] * fq * fq; 64 64 } 65 65 //now calculate outer integral 66 outer_sum += G auss76Wt[i] * inner_sum;66 outer_sum += GAUSS_W[i] * inner_sum; 67 67 } 68 68 -
sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c
r82592da rd4db147 7 7 double length) 8 8 { 9 return M_PI*( (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length + 9 return M_PI*( (r_minor + thick_rim)*(r_minor*x_core + thick_rim)* length + 10 10 square(r_minor)*x_core*2.0*thick_face ); 11 11 } … … 47 47 //initialize integral 48 48 double outer_sum = 0.0; 49 for(int i=0;i< 76;i++) {49 for(int i=0;i<GAUSS_N;i++) { 50 50 //setup inner integral over the ellipsoidal cross-section 51 51 // since we generate these lots of times, why not store them somewhere? 52 //const double cos_alpha = ( G auss76Z[i]*(vb-va) + va + vb )/2.0;53 const double cos_alpha = ( G auss76Z[i] + 1.0 )/2.0;52 //const double cos_alpha = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 53 const double cos_alpha = ( GAUSS_Z[i] + 1.0 )/2.0; 54 54 const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 55 55 double inner_sum=0; … … 58 58 si1 = sas_sinx_x(sinarg1); 59 59 si2 = sas_sinx_x(sinarg2); 60 for(int j=0;j< 76;j++) {60 for(int j=0;j<GAUSS_N;j++) { 61 61 //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 62 //const double beta = ( G auss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0;63 const double beta = ( G auss76Z[j] +1.0)*M_PI_2;62 //const double beta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 63 const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2; 64 64 const double rr = sqrt(r2A - r2B*cos(beta)); 65 65 double besarg1 = q*rr*sin_alpha; … … 67 67 be1 = sas_2J1x_x(besarg1); 68 68 be2 = sas_2J1x_x(besarg2); 69 inner_sum += G auss76Wt[j] *square(dr1*si1*be1 +69 inner_sum += GAUSS_W[j] *square(dr1*si1*be1 + 70 70 dr2*si1*be2 + 71 71 dr3*si2*be1); 72 72 } 73 73 //now calculate outer integral 74 outer_sum += G auss76Wt[i] * inner_sum;74 outer_sum += GAUSS_W[i] * inner_sum; 75 75 } 76 76 -
sasmodels/models/core_shell_cylinder.c
rbecded3 r74768cb 30 30 const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 31 31 double total = 0.0; 32 for (int i=0; i< 76;i++) {32 for (int i=0; i<GAUSS_N ;i++) { 33 33 // translate a point in [-1,1] to a point in [0, pi/2] 34 //const double theta = ( G auss76Z[i]*(upper-lower) + upper + lower )/2.0;34 //const double theta = ( GAUSS_Z[i]*(upper-lower) + upper + lower )/2.0; 35 35 double sin_theta, cos_theta; 36 const double theta = G auss76Z[i]*M_PI_4 + M_PI_4;36 const double theta = GAUSS_Z[i]*M_PI_4 + M_PI_4; 37 37 SINCOS(theta, sin_theta, cos_theta); 38 38 const double qab = q*sin_theta; … … 40 40 const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 41 41 + _cyl(shell_vd, shell_r*qab, shell_h*qc); 42 total += G auss76Wt[i] * fq * fq * sin_theta;42 total += GAUSS_W[i] * fq * fq * sin_theta; 43 43 } 44 44 // translate dx in [-1,1] to dx in [lower,upper] -
sasmodels/models/core_shell_ellipsoid.c
rbecded3 r74768cb 59 59 const double b = 0.5; 60 60 double total = 0.0; //initialize intergral 61 for(int i=0;i< 76;i++) {62 const double cos_theta = G auss76Z[i]*m + b;61 for(int i=0;i<GAUSS_N;i++) { 62 const double cos_theta = GAUSS_Z[i]*m + b; 63 63 const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 64 64 double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta, … … 66 66 equat_shell, polar_shell, 67 67 sld_core_shell, sld_shell_solvent); 68 total += G auss76Wt[i] * fq * fq;68 total += GAUSS_W[i] * fq * fq; 69 69 } 70 70 total *= m; -
sasmodels/models/core_shell_parallelepiped.c
r4493288 ra261a83 60 60 // outer integral (with gauss points), integration limits = 0, 1 61 61 double outer_sum = 0; //initialize integral 62 for( int i=0; i< 76; i++) {63 const double cos_alpha = 0.5 * ( G auss76Z[i] + 1.0 );62 for( int i=0; i<GAUSS_N; i++) { 63 const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 64 64 const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 65 65 … … 68 68 const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 69 69 double inner_sum = 0.0; 70 for(int j=0; j< 76; j++) {71 const double beta = 0.5 * ( G auss76Z[j] + 1.0 );70 for(int j=0; j<GAUSS_N; j++) { 71 const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 ); 72 72 double sin_beta, cos_beta; 73 73 SINCOS(M_PI_2*beta, sin_beta, cos_beta); … … 89 89 #endif 90 90 91 inner_sum += G auss76Wt[j] * f * f;91 inner_sum += GAUSS_W[j] * f * f; 92 92 } 93 93 inner_sum *= 0.5; 94 94 // now sum up the outer integral 95 outer_sum += G auss76Wt[i] * inner_sum;95 outer_sum += GAUSS_W[i] * inner_sum; 96 96 } 97 97 outer_sum *= 0.5; -
sasmodels/models/cylinder.c
rbecded3 r74768cb 21 21 22 22 double total = 0.0; 23 for (int i=0; i< 76;i++) {24 const double theta = G auss76Z[i]*zm + zb;23 for (int i=0; i<GAUSS_N ;i++) { 24 const double theta = GAUSS_Z[i]*zm + zb; 25 25 double sin_theta, cos_theta; // slots to hold sincos function output 26 26 // theta (theta,phi) the projection of the cylinder on the detector plane 27 27 SINCOS(theta , sin_theta, cos_theta); 28 28 const double form = fq(q*sin_theta, q*cos_theta, radius, length); 29 total += G auss76Wt[i] * form * form * sin_theta;29 total += GAUSS_W[i] * form * form * sin_theta; 30 30 } 31 31 // translate dx in [-1,1] to dx in [lower,upper] -
sasmodels/models/ellipsoid.c
rbecded3 r74768cb 22 22 23 23 // translate a point in [-1,1] to a point in [0, 1] 24 // const double u = G auss76Z[i]*(upper-lower)/2 + (upper+lower)/2;24 // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 25 25 const double zm = 0.5; 26 26 const double zb = 0.5; 27 27 double total = 0.0; 28 for (int i=0;i< 76;i++) {29 const double u = G auss76Z[i]*zm + zb;28 for (int i=0;i<GAUSS_N;i++) { 29 const double u = GAUSS_Z[i]*zm + zb; 30 30 const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 31 31 const double f = sas_3j1x_x(q*r); 32 total += G auss76Wt[i] * f * f;32 total += GAUSS_W[i] * f * f; 33 33 } 34 34 // translate dx in [-1,1] to dx in [lower,upper] -
sasmodels/models/elliptical_cylinder.c
r82592da rd4db147 22 22 //initialize integral 23 23 double outer_sum = 0.0; 24 for(int i=0;i< 76;i++) {24 for(int i=0;i<GAUSS_N;i++) { 25 25 //setup inner integral over the ellipsoidal cross-section 26 const double cos_val = ( G auss76Z[i]*(vb-va) + va + vb )/2.0;26 const double cos_val = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 27 27 const double sin_val = sqrt(1.0 - cos_val*cos_val); 28 28 //const double arg = radius_minor*sin_val; 29 29 double inner_sum=0; 30 for(int j=0;j<76;j++) { 31 //20 gauss points for the inner integral, increase to 76, RKH 6Nov2017 32 const double theta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 30 for(int j=0;j<GAUSS_N;j++) { 31 const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 33 32 const double r = sin_val*sqrt(rA - rB*cos(theta)); 34 33 const double be = sas_2J1x_x(q*r); 35 inner_sum += G auss76Wt[j] * be * be;34 inner_sum += GAUSS_W[j] * be * be; 36 35 } 37 36 //now calculate the value of the inner integral … … 40 39 //now calculate outer integral 41 40 const double si = sas_sinx_x(q*0.5*length*cos_val); 42 outer_sum += G auss76Wt[i] * inner_sum * si * si;41 outer_sum += GAUSS_W[i] * inner_sum * si * si; 43 42 } 44 43 outer_sum *= 0.5*(vb-va); -
sasmodels/models/elliptical_cylinder.py
r2d81cfe ra261a83 121 121 # pylint: enable=bad-whitespace, line-too-long 122 122 123 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "lib/gauss20.c", 124 "elliptical_cylinder.c"] 123 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 125 124 126 125 demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, -
sasmodels/models/fcc_paracrystal.c
rf728001 r74768cb 53 53 54 54 double outer_sum = 0.0; 55 for(int i=0; i< 150; i++) {55 for(int i=0; i<GAUSS_N; i++) { 56 56 double inner_sum = 0.0; 57 const double theta = G auss150Z[i]*theta_m + theta_b;57 const double theta = GAUSS_Z[i]*theta_m + theta_b; 58 58 double sin_theta, cos_theta; 59 59 SINCOS(theta, sin_theta, cos_theta); 60 60 const double qc = q*cos_theta; 61 61 const double qab = q*sin_theta; 62 for(int j=0;j< 150;j++) {63 const double phi = G auss150Z[j]*phi_m + phi_b;62 for(int j=0;j<GAUSS_N;j++) { 63 const double phi = GAUSS_Z[j]*phi_m + phi_b; 64 64 double sin_phi, cos_phi; 65 65 SINCOS(phi, sin_phi, cos_phi); … … 67 67 const double qb = qab*sin_phi; 68 68 const double form = fcc_Zq(qa, qb, qc, dnn, d_factor); 69 inner_sum += G auss150Wt[j] * form;69 inner_sum += GAUSS_W[j] * form; 70 70 } 71 71 inner_sum *= phi_m; // sum(f(x)dx) = sum(f(x)) dx 72 outer_sum += G auss150Wt[i] * inner_sum * sin_theta;72 outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 73 73 } 74 74 outer_sum *= theta_m; -
sasmodels/models/flexible_cylinder_elliptical.c
r592343f r74768cb 17 17 double sum=0.0; 18 18 19 for(int i=0;i< N_POINTS_76;i++) {20 const double zi = ( G auss76Z[i] + 1.0 )*M_PI_4;19 for(int i=0;i<GAUSS_N;i++) { 20 const double zi = ( GAUSS_Z[i] + 1.0 )*M_PI_4; 21 21 double sn, cn; 22 22 SINCOS(zi, sn, cn); 23 23 const double arg = q*sqrt(a*a*sn*sn + b*b*cn*cn); 24 24 const double yyy = sas_2J1x_x(arg); 25 sum += G auss76Wt[i] * yyy * yyy;25 sum += GAUSS_W[i] * yyy * yyy; 26 26 } 27 27 sum *= 0.5; -
sasmodels/models/hollow_cylinder.c
rbecded3 r74768cb 38 38 39 39 double summ = 0.0; //initialize intergral 40 for (int i=0;i< 76;i++) {41 const double cos_theta = 0.5*( G auss76Z[i] * (upper-lower) + lower + upper );40 for (int i=0;i<GAUSS_N;i++) { 41 const double cos_theta = 0.5*( GAUSS_Z[i] * (upper-lower) + lower + upper ); 42 42 const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 43 43 const double form = _fq(q*sin_theta, q*cos_theta, 44 44 radius, thickness, length); 45 summ += G auss76Wt[i] * form * form;45 summ += GAUSS_W[i] * form * form; 46 46 } 47 47 -
sasmodels/models/hollow_rectangular_prism.c
r8de1477 r74768cb 39 39 40 40 double outer_sum = 0.0; 41 for(int i=0; i< 76; i++) {41 for(int i=0; i<GAUSS_N; i++) { 42 42 43 const double theta = 0.5 * ( G auss76Z[i]*(v1b-v1a) + v1a + v1b );43 const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 44 44 double sin_theta, cos_theta; 45 45 SINCOS(theta, sin_theta, cos_theta); … … 49 49 50 50 double inner_sum = 0.0; 51 for(int j=0; j< 76; j++) {51 for(int j=0; j<GAUSS_N; j++) { 52 52 53 const double phi = 0.5 * ( G auss76Z[j]*(v2b-v2a) + v2a + v2b );53 const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 54 54 double sin_phi, cos_phi; 55 55 SINCOS(phi, sin_phi, cos_phi); … … 66 66 const double AP2 = vol_core * termA2 * termB2 * termC2; 67 67 68 inner_sum += G auss76Wt[j] * square(AP1-AP2);68 inner_sum += GAUSS_W[j] * square(AP1-AP2); 69 69 } 70 70 inner_sum *= 0.5 * (v2b-v2a); 71 71 72 outer_sum += G auss76Wt[i] * inner_sum * sin(theta);72 outer_sum += GAUSS_W[i] * inner_sum * sin(theta); 73 73 } 74 74 outer_sum *= 0.5*(v1b-v1a); -
sasmodels/models/hollow_rectangular_prism_thin_walls.c
rab2aea8 r74768cb 1 1 double form_volume(double length_a, double b2a_ratio, double c2a_ratio); 2 double Iq(double q, double sld, double solvent_sld, double length_a, 2 double Iq(double q, double sld, double solvent_sld, double length_a, 3 3 double b2a_ratio, double c2a_ratio); 4 4 … … 29 29 const double v2a = 0.0; 30 30 const double v2b = M_PI_2; //phi integration limits 31 31 32 32 double outer_sum = 0.0; 33 for(int i=0; i< 76; i++) {34 const double theta = 0.5 * ( G auss76Z[i]*(v1b-v1a) + v1a + v1b );33 for(int i=0; i<GAUSS_N; i++) { 34 const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 35 35 36 36 double sin_theta, cos_theta; … … 44 44 45 45 double inner_sum = 0.0; 46 for(int j=0; j< 76; j++) {47 const double phi = 0.5 * ( G auss76Z[j]*(v2b-v2a) + v2a + v2b );46 for(int j=0; j<GAUSS_N; j++) { 47 const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 48 48 49 49 double sin_phi, cos_phi; … … 62 62 * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi ); 63 63 64 inner_sum += G auss76Wt[j] * square(AL+AT);64 inner_sum += GAUSS_W[j] * square(AL+AT); 65 65 } 66 66 67 67 inner_sum *= 0.5 * (v2b-v2a); 68 outer_sum += G auss76Wt[i] * inner_sum * sin_theta;68 outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 69 69 } 70 70 -
sasmodels/models/lib/gauss150.c
r994d77f r74768cb 7 7 * 8 8 */ 9 #ifdef GAUSS_N 10 # undef GAUSS_N 11 # undef GAUSS_Z 12 # undef GAUSS_W 13 #endif 14 #define GAUSS_N 150 15 #define GAUSS_Z Gauss150Z 16 #define GAUSS_W Gauss150Wt 17 18 // Note: using array size 152 so that it is a multiple of 4 9 19 10 20 // Gaussians 11 constant double Gauss150Z[15 0]={21 constant double Gauss150Z[152]={ 12 22 -0.9998723404457334, 13 23 -0.9993274305065947, … … 159 169 0.9983473449340834, 160 170 0.9993274305065947, 161 0.9998723404457334 171 0.9998723404457334, 172 0., 173 0. 162 174 }; 163 175 164 constant double Gauss150Wt[15 0]={176 constant double Gauss150Wt[152]={ 165 177 0.0003276086705538, 166 178 0.0007624720924706, … … 312 324 0.0011976474864367, 313 325 0.0007624720924706, 314 0.0003276086705538 326 0.0003276086705538, 327 0., 328 0. 315 329 }; -
sasmodels/models/lib/gauss20.c
r994d77f r74768cb 7 7 * 8 8 */ 9 #ifdef GAUSS_N 10 # undef GAUSS_N 11 # undef GAUSS_Z 12 # undef GAUSS_W 13 #endif 14 #define GAUSS_N 20 15 #define GAUSS_Z Gauss20Z 16 #define GAUSS_W Gauss20Wt 9 17 10 18 // Gaussians -
sasmodels/models/lib/gauss76.c
r66d119f r74768cb 7 7 * 8 8 */ 9 #define N_POINTS_76 76 9 #ifdef GAUSS_N 10 # undef GAUSS_N 11 # undef GAUSS_Z 12 # undef GAUSS_W 13 #endif 14 #define GAUSS_N 76 15 #define GAUSS_Z Gauss76Z 16 #define GAUSS_W Gauss76Wt 10 17 11 18 // Gaussians 12 constant double Gauss76Wt[ N_POINTS_76]={19 constant double Gauss76Wt[76]={ 13 20 .00126779163408536, //0 14 21 .00294910295364247, … … 89 96 }; 90 97 91 constant double Gauss76Z[ N_POINTS_76]={98 constant double Gauss76Z[76]={ 92 99 -.999505948362153, //0 93 100 -.997397786355355, -
sasmodels/models/parallelepiped.c
r9b7b23f r74768cb 23 23 double outer_total = 0; //initialize integral 24 24 25 for( int i=0; i< 76; i++) {26 const double sigma = 0.5 * ( G auss76Z[i] + 1.0 );25 for( int i=0; i<GAUSS_N; i++) { 26 const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 ); 27 27 const double mu_proj = mu * sqrt(1.0-sigma*sigma); 28 28 … … 30 30 // corresponding to angles from 0 to pi/2. 31 31 double inner_total = 0.0; 32 for(int j=0; j< 76; j++) {33 const double uu = 0.5 * ( G auss76Z[j] + 1.0 );32 for(int j=0; j<GAUSS_N; j++) { 33 const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 ); 34 34 double sin_uu, cos_uu; 35 35 SINCOS(M_PI_2*uu, sin_uu, cos_uu); 36 36 const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 37 37 const double si2 = sas_sinx_x(mu_proj * cos_uu); 38 inner_total += G auss76Wt[j] * square(si1 * si2);38 inner_total += GAUSS_W[j] * square(si1 * si2); 39 39 } 40 40 inner_total *= 0.5; 41 41 42 42 const double si = sas_sinx_x(mu * c_scaled * sigma); 43 outer_total += G auss76Wt[i] * inner_total * si * si;43 outer_total += GAUSS_W[i] * inner_total * si * si; 44 44 } 45 45 outer_total *= 0.5; -
sasmodels/models/pringle.c
r1e7b0db0 r74768cb 29 29 double sumC = 0.0; // initialize integral 30 30 double r; 31 for (int i=0; i < 76; i++) {32 r = G auss76Z[i]*zm + zb;31 for (int i=0; i < GAUSS_N; i++) { 32 r = GAUSS_Z[i]*zm + zb; 33 33 34 34 const double qrs = r*q_sin_psi; 35 35 const double qrrc = r*r*q_cos_psi; 36 36 37 double y = G auss76Wt[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs);37 double y = GAUSS_W[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs); 38 38 double S, C; 39 39 SINCOS(alpha*qrrc, S, C); … … 86 86 87 87 double sum = 0.0; 88 for (int i = 0; i < 76; i++) {89 double psi = G auss76Z[i]*zm + zb;88 for (int i = 0; i < GAUSS_N; i++) { 89 double psi = GAUSS_Z[i]*zm + zb; 90 90 double sin_psi, cos_psi; 91 91 SINCOS(psi, sin_psi, cos_psi); … … 93 93 double sinc_term = square(sas_sinx_x(q * thickness * cos_psi / 2.0)); 94 94 double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term; 95 sum += G auss76Wt[i] * pringle_kernel;95 sum += GAUSS_W[i] * pringle_kernel; 96 96 } 97 97 -
sasmodels/models/rectangular_prism.c
r8de1477 r74768cb 28 28 29 29 double outer_sum = 0.0; 30 for(int i=0; i< 76; i++) {31 const double theta = 0.5 * ( G auss76Z[i]*(v1b-v1a) + v1a + v1b );30 for(int i=0; i<GAUSS_N; i++) { 31 const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 32 32 double sin_theta, cos_theta; 33 33 SINCOS(theta, sin_theta, cos_theta); … … 36 36 37 37 double inner_sum = 0.0; 38 for(int j=0; j< 76; j++) {39 double phi = 0.5 * ( G auss76Z[j]*(v2b-v2a) + v2a + v2b );38 for(int j=0; j<GAUSS_N; j++) { 39 double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 40 40 double sin_phi, cos_phi; 41 41 SINCOS(phi, sin_phi, cos_phi); … … 45 45 const double termB = sas_sinx_x(q * b_half * sin_theta * cos_phi); 46 46 const double AP = termA * termB * termC; 47 inner_sum += G auss76Wt[j] * AP * AP;47 inner_sum += GAUSS_W[j] * AP * AP; 48 48 } 49 49 inner_sum = 0.5 * (v2b-v2a) * inner_sum; 50 outer_sum += G auss76Wt[i] * inner_sum * sin_theta;50 outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 51 51 } 52 52 -
sasmodels/models/sc_paracrystal.c
rf728001 r74768cb 54 54 55 55 double outer_sum = 0.0; 56 for(int i=0; i< 150; i++) {56 for(int i=0; i<GAUSS_N; i++) { 57 57 double inner_sum = 0.0; 58 const double theta = G auss150Z[i]*theta_m + theta_b;58 const double theta = GAUSS_Z[i]*theta_m + theta_b; 59 59 double sin_theta, cos_theta; 60 60 SINCOS(theta, sin_theta, cos_theta); 61 61 const double qc = q*cos_theta; 62 62 const double qab = q*sin_theta; 63 for(int j=0;j< 150;j++) {64 const double phi = G auss150Z[j]*phi_m + phi_b;63 for(int j=0;j<GAUSS_N;j++) { 64 const double phi = GAUSS_Z[j]*phi_m + phi_b; 65 65 double sin_phi, cos_phi; 66 66 SINCOS(phi, sin_phi, cos_phi); … … 68 68 const double qb = qab*sin_phi; 69 69 const double form = sc_Zq(qa, qb, qc, dnn, d_factor); 70 inner_sum += G auss150Wt[j] * form;70 inner_sum += GAUSS_W[j] * form; 71 71 } 72 72 inner_sum *= phi_m; // sum(f(x)dx) = sum(f(x)) dx 73 outer_sum += G auss150Wt[i] * inner_sum * sin_theta;73 outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 74 74 } 75 75 outer_sum *= theta_m; -
sasmodels/models/stacked_disks.c
rbecded3 r74768cb 81 81 double halfheight = 0.5*thick_core; 82 82 83 for(int i=0; i< N_POINTS_76; i++) {84 double zi = (G auss76Z[i] + 1.0)*M_PI_4;83 for(int i=0; i<GAUSS_N; i++) { 84 double zi = (GAUSS_Z[i] + 1.0)*M_PI_4; 85 85 double sin_alpha, cos_alpha; // slots to hold sincos function output 86 86 SINCOS(zi, sin_alpha, cos_alpha); … … 95 95 solvent_sld, 96 96 d); 97 summ += G auss76Wt[i] * yyy * sin_alpha;97 summ += GAUSS_W[i] * yyy * sin_alpha; 98 98 } 99 99 -
sasmodels/models/triaxial_ellipsoid.c
rbecded3 r74768cb 21 21 const double zb = M_PI_4; 22 22 double outer = 0.0; 23 for (int i=0;i< 76;i++) {24 //const double u = G auss76Z[i]*(upper-lower)/2 + (upper + lower)/2;25 const double phi = G auss76Z[i]*zm + zb;23 for (int i=0;i<GAUSS_N;i++) { 24 //const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper + lower)/2; 25 const double phi = GAUSS_Z[i]*zm + zb; 26 26 const double pa_sinsq_phi = pa*square(sin(phi)); 27 27 … … 29 29 const double um = 0.5; 30 30 const double ub = 0.5; 31 for (int j=0;j< 76;j++) {31 for (int j=0;j<GAUSS_N;j++) { 32 32 // translate a point in [-1,1] to a point in [0, 1] 33 const double usq = square(G auss76Z[j]*um + ub);33 const double usq = square(GAUSS_Z[j]*um + ub); 34 34 const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0-usq) + 1.0 + pc*usq); 35 35 const double fq = sas_3j1x_x(q*r); 36 inner += G auss76Wt[j] * fq * fq;36 inner += GAUSS_W[j] * fq * fq; 37 37 } 38 outer += G auss76Wt[i] * inner; // correcting for dx later38 outer += GAUSS_W[i] * inner; // correcting for dx later 39 39 } 40 40 // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi
Note: See TracChangeset
for help on using the changeset viewer.