Changes in / [d4db147:d8ac2ad] in sasmodels
- Location:
- sasmodels
- Files:
-
- 1 deleted
- 29 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/compare.py
rff31782 r376b0ee 44 44 from .direct_model import DirectModel, get_mesh 45 45 from .convert import revert_name, revert_pars, constrain_new_to_old 46 from .generate import FLOAT_RE , set_integration_size46 from .generate import FLOAT_RE 47 47 from .weights import plot_weights 48 48 … … 801 801 return data, index 802 802 803 def make_engine(model_info, data, dtype, cutoff , ngauss=0):803 def make_engine(model_info, data, dtype, cutoff): 804 804 # type: (ModelInfo, Data, str, float) -> Calculator 805 805 """ … … 809 809 than OpenCL. 810 810 """ 811 if ngauss:812 set_integration_size(model_info, ngauss)813 814 811 if dtype == 'sasview': 815 812 return eval_sasview(model_info, data) … … 1048 1045 'poly', 'mono', 'cutoff=', 1049 1046 'magnetic', 'nonmagnetic', 1050 'accuracy=', 'ngauss=',1047 'accuracy=', 1051 1048 'neval=', # for timing... 1052 1049 … … 1182 1179 'show_weights' : False, 1183 1180 'sphere' : 0, 1184 'ngauss' : '0',1185 1181 } 1186 1182 for arg in flags: … … 1209 1205 elif arg.startswith('-engine='): opts['engine'] = arg[8:] 1210 1206 elif arg.startswith('-neval='): opts['count'] = arg[7:] 1211 elif arg.startswith('-ngauss='): opts['ngauss'] = arg[8:]1212 1207 elif arg.startswith('-random='): 1213 1208 opts['seed'] = int(arg[8:]) … … 1265 1260 1266 1261 comparison = any(PAR_SPLIT in v for v in values) 1267 1268 1262 if PAR_SPLIT in name: 1269 1263 names = name.split(PAR_SPLIT, 2) … … 1278 1272 return None 1279 1273 1280 if PAR_SPLIT in opts['ngauss']:1281 opts['ngauss'] = [int(k) for k in opts['ngauss'].split(PAR_SPLIT, 2)]1282 comparison = True1283 else:1284 opts['ngauss'] = [int(opts['ngauss'])]*21285 1286 1274 if PAR_SPLIT in opts['engine']: 1287 1275 opts['engine'] = opts['engine'].split(PAR_SPLIT, 2) … … 1302 1290 opts['cutoff'] = [float(opts['cutoff'])]*2 1303 1291 1304 base = make_engine(model_info[0], data, opts['engine'][0], 1305 opts['cutoff'][0], opts['ngauss'][0]) 1292 base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0]) 1306 1293 if comparison: 1307 comp = make_engine(model_info[1], data, opts['engine'][1], 1308 opts['cutoff'][1], opts['ngauss'][1]) 1294 comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1]) 1309 1295 else: 1310 1296 comp = None -
sasmodels/generate.py
rff31782 rff10479 268 268 """ 269 269 270 271 def set_integration_size(info, n):272 # type: (ModelInfo, int) -> None273 """274 Update the model definition, replacing the gaussian integration with275 a gaussian integration of a different size.276 277 Note: this really ought to be a method in modelinfo, but that leads to278 import loops.279 """280 if (info.source and any(lib.startswith('lib/gauss') for lib in info.source)):281 import os.path282 from .gengauss import gengauss283 path = os.path.join(MODEL_PATH, "lib", "gauss%d.c"%n)284 if not os.path.exists(path):285 gengauss(n, path)286 info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss')287 else lib for lib in info.source]288 270 289 271 def format_units(units): -
sasmodels/models/barbell.c
r74768cb rbecded3 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 < GAUSS_N; i++){26 const double t = G AUSS_Z[i]*zm + zb;25 for (int i = 0; i < 76; i++){ 26 const double t = Gauss76Z[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 AUSS_W[i] * Fq;30 total += Gauss76Wt[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 < GAUSS_N; i++){76 const double alpha= G AUSS_Z[i]*zm + zb;75 for (int i = 0; i < 76; i++){ 76 const double alpha= Gauss76Z[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 AUSS_W[i] * Aq * Aq * sin_alpha;80 total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 81 81 } 82 82 // translate dx in [-1,1] to dx in [lower,upper] -
sasmodels/models/bcc_paracrystal.c
r74768cb rea60e08 81 81 82 82 double outer_sum = 0.0; 83 for(int i=0; i< GAUSS_N; i++) {83 for(int i=0; i<150; i++) { 84 84 double inner_sum = 0.0; 85 const double theta = G AUSS_Z[i]*theta_m + theta_b;85 const double theta = Gauss150Z[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< GAUSS_N;j++) {91 const double phi = G AUSS_Z[j]*phi_m + phi_b;90 for(int j=0;j<150;j++) { 91 const double phi = Gauss150Z[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 AUSS_W[j] * form;97 inner_sum += Gauss150Wt[j] * form; 98 98 } 99 99 inner_sum *= phi_m; // sum(f(x)dx) = sum(f(x)) dx 100 outer_sum += G AUSS_W[i] * inner_sum * sin_theta;100 outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 101 101 } 102 102 outer_sum *= theta_m; -
sasmodels/models/capped_cylinder.c
r74768cb rbecded3 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< GAUSS_N;i++) {33 const double t = G AUSS_Z[i]*zm + zb;32 for (int i=0; i<76 ;i++) { 33 const double t = Gauss76Z[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 AUSS_W[i] * Fq;37 total += Gauss76Wt[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< GAUSS_N;i++) {98 const double theta = G AUSS_Z[i]*zm + zb;97 for (int i=0; i<76 ;i++) { 98 const double theta = Gauss76Z[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 AUSS_W[i] * Aq * Aq * sin_theta;105 total += Gauss76Wt[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
r74768cb rbecded3 52 52 53 53 double total = 0.0; 54 for(int i=0;i< GAUSS_N;i++) {55 double theta = (G AUSS_Z[i] + 1.0)*uplim;54 for(int i=0;i<N_POINTS_76;i++) { 55 double theta = (Gauss76Z[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 AUSS_W[i]*fq*fq*sin_theta;60 total += Gauss76Wt[i]*fq*fq*sin_theta; 61 61 } 62 62 -
sasmodels/models/core_shell_bicelle_elliptical.c
r74768cb r82592da 37 37 //initialize integral 38 38 double outer_sum = 0.0; 39 for(int i=0;i< GAUSS_N;i++) {39 for(int i=0;i<76;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 AUSS_Z[i]*(vb-va) + va + vb )/2.0;44 const double cos_theta = ( G AUSS_Z[i] + 1.0 )/2.0;43 //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 44 const double cos_theta = ( Gauss76Z[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< GAUSS_N;j++) {51 for(int j=0;j<76;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 AUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;57 const double phi = ( G AUSS_Z[j] +1.0)*M_PI_2;56 //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 57 const double phi = ( Gauss76Z[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 AUSS_W[j] * fq * fq;63 inner_sum += Gauss76Wt[j] * fq * fq; 64 64 } 65 65 //now calculate outer integral 66 outer_sum += G AUSS_W[i] * inner_sum;66 outer_sum += Gauss76Wt[i] * inner_sum; 67 67 } 68 68 -
sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c
r74768cb r82592da 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< GAUSS_N;i++) {49 for(int i=0;i<76;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 AUSS_Z[i]*(vb-va) + va + vb )/2.0;53 const double cos_alpha = ( G AUSS_Z[i] + 1.0 )/2.0;52 //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 53 const double cos_alpha = ( Gauss76Z[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< GAUSS_N;j++) {60 for(int j=0;j<76;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 AUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0;63 const double beta = ( G AUSS_Z[j] +1.0)*M_PI_2;62 //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 63 const double beta = ( Gauss76Z[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 AUSS_W[j] *square(dr1*si1*be1 +69 inner_sum += Gauss76Wt[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 AUSS_W[i] * inner_sum;74 outer_sum += Gauss76Wt[i] * inner_sum; 75 75 } 76 76 -
sasmodels/models/core_shell_cylinder.c
r74768cb rbecded3 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< GAUSS_N;i++) {32 for (int i=0; i<76 ;i++) { 33 33 // translate a point in [-1,1] to a point in [0, pi/2] 34 //const double theta = ( G AUSS_Z[i]*(upper-lower) + upper + lower )/2.0;34 //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 35 35 double sin_theta, cos_theta; 36 const double theta = G AUSS_Z[i]*M_PI_4 + M_PI_4;36 const double theta = Gauss76Z[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 AUSS_W[i] * fq * fq * sin_theta;42 total += Gauss76Wt[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
r74768cb rbecded3 59 59 const double b = 0.5; 60 60 double total = 0.0; //initialize intergral 61 for(int i=0;i< GAUSS_N;i++) {62 const double cos_theta = G AUSS_Z[i]*m + b;61 for(int i=0;i<76;i++) { 62 const double cos_theta = Gauss76Z[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 AUSS_W[i] * fq * fq;68 total += Gauss76Wt[i] * fq * fq; 69 69 } 70 70 total *= m; -
sasmodels/models/core_shell_parallelepiped.c
r74768cb rfc0b7aa 80 80 // outer integral (with gauss points), integration limits = 0, 1 81 81 double outer_sum = 0; //initialize integral 82 for( int i=0; i< GAUSS_N; i++) {83 double sigma = 0.5 * ( G AUSS_Z[i] + 1.0 );82 for( int i=0; i<76; i++) { 83 double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 84 84 double mu_proj = mu * sqrt(1.0-sigma*sigma); 85 85 … … 88 88 const double siCt = sas_sinx_x(mu * sigma * tC_over_b); 89 89 double inner_sum = 0.0; 90 for(int j=0; j< GAUSS_N; j++) {91 const double uu = 0.5 * ( G AUSS_Z[j] + 1.0 );90 for(int j=0; j<76; j++) { 91 const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 92 92 double sin_uu, cos_uu; 93 93 SINCOS(M_PI_2*uu, sin_uu, cos_uu); … … 109 109 #endif 110 110 111 inner_sum += G AUSS_W[j] * f * f;111 inner_sum += Gauss76Wt[j] * f * f; 112 112 } 113 113 inner_sum *= 0.5; 114 114 // now sum up the outer integral 115 outer_sum += G AUSS_W[i] * inner_sum;115 outer_sum += Gauss76Wt[i] * inner_sum; 116 116 } 117 117 outer_sum *= 0.5; -
sasmodels/models/cylinder.c
r74768cb rbecded3 21 21 22 22 double total = 0.0; 23 for (int i=0; i< GAUSS_N;i++) {24 const double theta = G AUSS_Z[i]*zm + zb;23 for (int i=0; i<76 ;i++) { 24 const double theta = Gauss76Z[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 AUSS_W[i] * form * form * sin_theta;29 total += Gauss76Wt[i] * form * form * sin_theta; 30 30 } 31 31 // translate dx in [-1,1] to dx in [lower,upper] -
sasmodels/models/ellipsoid.c
r74768cb rbecded3 22 22 23 23 // translate a point in [-1,1] to a point in [0, 1] 24 // const double u = G AUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2;24 // const double u = Gauss76Z[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< GAUSS_N;i++) {29 const double u = G AUSS_Z[i]*zm + zb;28 for (int i=0;i<76;i++) { 29 const double u = Gauss76Z[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 AUSS_W[i] * f * f;32 total += Gauss76Wt[i] * f * f; 33 33 } 34 34 // translate dx in [-1,1] to dx in [lower,upper] -
sasmodels/models/elliptical_cylinder.c
r74768cb r82592da 22 22 //initialize integral 23 23 double outer_sum = 0.0; 24 for(int i=0;i< GAUSS_N;i++) {24 for(int i=0;i<76;i++) { 25 25 //setup inner integral over the ellipsoidal cross-section 26 const double cos_val = ( G AUSS_Z[i]*(vb-va) + va + vb )/2.0;26 const double cos_val = ( Gauss76Z[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<GAUSS_N;j++) { 31 const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.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; 32 33 const double r = sin_val*sqrt(rA - rB*cos(theta)); 33 34 const double be = sas_2J1x_x(q*r); 34 inner_sum += G AUSS_W[j] * be * be;35 inner_sum += Gauss76Wt[j] * be * be; 35 36 } 36 37 //now calculate the value of the inner integral … … 39 40 //now calculate outer integral 40 41 const double si = sas_sinx_x(q*0.5*length*cos_val); 41 outer_sum += G AUSS_W[i] * inner_sum * si * si;42 outer_sum += Gauss76Wt[i] * inner_sum * si * si; 42 43 } 43 44 outer_sum *= 0.5*(vb-va); -
sasmodels/models/elliptical_cylinder.py
r74768cb reda8b30 43 43 P(q) = \text{scale} <F^2> / V 44 44 45 For 2d data the orientation of the particle is required, described using a different set 46 of angles as in the diagrams below, for further details of the calculation and angular 45 For 2d data the orientation of the particle is required, described using a different set 46 of angles as in the diagrams below, for further details of the calculation and angular 47 47 dispersions see :ref:`orientation` . 48 48 … … 121 121 # pylint: enable=bad-whitespace, line-too-long 122 122 123 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 123 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "lib/gauss20.c", 124 "elliptical_cylinder.c"] 124 125 125 126 demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, -
sasmodels/models/fcc_paracrystal.c
r74768cb rf728001 53 53 54 54 double outer_sum = 0.0; 55 for(int i=0; i< GAUSS_N; i++) {55 for(int i=0; i<150; i++) { 56 56 double inner_sum = 0.0; 57 const double theta = G AUSS_Z[i]*theta_m + theta_b;57 const double theta = Gauss150Z[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< GAUSS_N;j++) {63 const double phi = G AUSS_Z[j]*phi_m + phi_b;62 for(int j=0;j<150;j++) { 63 const double phi = Gauss150Z[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 AUSS_W[j] * form;69 inner_sum += Gauss150Wt[j] * form; 70 70 } 71 71 inner_sum *= phi_m; // sum(f(x)dx) = sum(f(x)) dx 72 outer_sum += G AUSS_W[i] * inner_sum * sin_theta;72 outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 73 73 } 74 74 outer_sum *= theta_m; -
sasmodels/models/flexible_cylinder_elliptical.c
r74768cb r592343f 17 17 double sum=0.0; 18 18 19 for(int i=0;i< GAUSS_N;i++) {20 const double zi = ( G AUSS_Z[i] + 1.0 )*M_PI_4;19 for(int i=0;i<N_POINTS_76;i++) { 20 const double zi = ( Gauss76Z[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 AUSS_W[i] * yyy * yyy;25 sum += Gauss76Wt[i] * yyy * yyy; 26 26 } 27 27 sum *= 0.5; -
sasmodels/models/hollow_cylinder.c
r74768cb rbecded3 38 38 39 39 double summ = 0.0; //initialize intergral 40 for (int i=0;i< GAUSS_N;i++) {41 const double cos_theta = 0.5*( G AUSS_Z[i] * (upper-lower) + lower + upper );40 for (int i=0;i<76;i++) { 41 const double cos_theta = 0.5*( Gauss76Z[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 AUSS_W[i] * form * form;45 summ += Gauss76Wt[i] * form * form; 46 46 } 47 47 -
sasmodels/models/hollow_rectangular_prism.c
r74768cb r8de1477 39 39 40 40 double outer_sum = 0.0; 41 for(int i=0; i< GAUSS_N; i++) {41 for(int i=0; i<76; i++) { 42 42 43 const double theta = 0.5 * ( G AUSS_Z[i]*(v1b-v1a) + v1a + v1b );43 const double theta = 0.5 * ( Gauss76Z[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< GAUSS_N; j++) {51 for(int j=0; j<76; j++) { 52 52 53 const double phi = 0.5 * ( G AUSS_Z[j]*(v2b-v2a) + v2a + v2b );53 const double phi = 0.5 * ( Gauss76Z[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 AUSS_W[j] * square(AP1-AP2);68 inner_sum += Gauss76Wt[j] * square(AP1-AP2); 69 69 } 70 70 inner_sum *= 0.5 * (v2b-v2a); 71 71 72 outer_sum += G AUSS_W[i] * inner_sum * sin(theta);72 outer_sum += Gauss76Wt[i] * inner_sum * sin(theta); 73 73 } 74 74 outer_sum *= 0.5*(v1b-v1a); -
sasmodels/models/hollow_rectangular_prism_thin_walls.c
r74768cb rab2aea8 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< GAUSS_N; i++) {34 const double theta = 0.5 * ( G AUSS_Z[i]*(v1b-v1a) + v1a + v1b );33 for(int i=0; i<76; i++) { 34 const double theta = 0.5 * ( Gauss76Z[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< GAUSS_N; j++) {47 const double phi = 0.5 * ( G AUSS_Z[j]*(v2b-v2a) + v2a + v2b );46 for(int j=0; j<76; j++) { 47 const double phi = 0.5 * ( Gauss76Z[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 AUSS_W[j] * square(AL+AT);64 inner_sum += Gauss76Wt[j] * square(AL+AT); 65 65 } 66 66 67 67 inner_sum *= 0.5 * (v2b-v2a); 68 outer_sum += G AUSS_W[i] * inner_sum * sin_theta;68 outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 69 69 } 70 70 -
sasmodels/models/lib/gauss150.c
r74768cb r994d77f 7 7 * 8 8 */ 9 #ifdef GAUSS_N10 # undef GAUSS_N11 # undef GAUSS_Z12 # undef GAUSS_W13 #endif14 #define GAUSS_N 15015 #define GAUSS_Z Gauss150Z16 #define GAUSS_W Gauss150Wt17 18 // Note: using array size 152 so that it is a multiple of 419 9 20 10 // Gaussians 21 constant double Gauss150Z[15 2]={11 constant double Gauss150Z[150]={ 22 12 -0.9998723404457334, 23 13 -0.9993274305065947, … … 169 159 0.9983473449340834, 170 160 0.9993274305065947, 171 0.9998723404457334, 172 0., 173 0. 161 0.9998723404457334 174 162 }; 175 163 176 constant double Gauss150Wt[15 2]={164 constant double Gauss150Wt[150]={ 177 165 0.0003276086705538, 178 166 0.0007624720924706, … … 324 312 0.0011976474864367, 325 313 0.0007624720924706, 326 0.0003276086705538, 327 0., 328 0. 314 0.0003276086705538 329 315 }; -
sasmodels/models/lib/gauss20.c
r74768cb r994d77f 7 7 * 8 8 */ 9 #ifdef GAUSS_N10 # undef GAUSS_N11 # undef GAUSS_Z12 # undef GAUSS_W13 #endif14 #define GAUSS_N 2015 #define GAUSS_Z Gauss20Z16 #define GAUSS_W Gauss20Wt17 9 18 10 // Gaussians -
sasmodels/models/lib/gauss76.c
r74768cb r66d119f 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 76 15 #define GAUSS_Z Gauss76Z 16 #define GAUSS_W Gauss76Wt 9 #define N_POINTS_76 76 17 10 18 11 // Gaussians 19 constant double Gauss76Wt[ 76]={12 constant double Gauss76Wt[N_POINTS_76]={ 20 13 .00126779163408536, //0 21 14 .00294910295364247, … … 96 89 }; 97 90 98 constant double Gauss76Z[ 76]={91 constant double Gauss76Z[N_POINTS_76]={ 99 92 -.999505948362153, //0 100 93 -.997397786355355, -
sasmodels/models/parallelepiped.c
r74768cb r9b7b23f 23 23 double outer_total = 0; //initialize integral 24 24 25 for( int i=0; i< GAUSS_N; i++) {26 const double sigma = 0.5 * ( G AUSS_Z[i] + 1.0 );25 for( int i=0; i<76; i++) { 26 const double sigma = 0.5 * ( Gauss76Z[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< GAUSS_N; j++) {33 const double uu = 0.5 * ( G AUSS_Z[j] + 1.0 );32 for(int j=0; j<76; j++) { 33 const double uu = 0.5 * ( Gauss76Z[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 AUSS_W[j] * square(si1 * si2);38 inner_total += Gauss76Wt[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 AUSS_W[i] * inner_total * si * si;43 outer_total += Gauss76Wt[i] * inner_total * si * si; 44 44 } 45 45 outer_total *= 0.5; -
sasmodels/models/pringle.c
r74768cb r1e7b0db0 29 29 double sumC = 0.0; // initialize integral 30 30 double r; 31 for (int i=0; i < GAUSS_N; i++) {32 r = G AUSS_Z[i]*zm + zb;31 for (int i=0; i < 76; i++) { 32 r = Gauss76Z[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 AUSS_W[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs);37 double y = Gauss76Wt[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 < GAUSS_N; i++) {89 double psi = G AUSS_Z[i]*zm + zb;88 for (int i = 0; i < 76; i++) { 89 double psi = Gauss76Z[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 AUSS_W[i] * pringle_kernel;95 sum += Gauss76Wt[i] * pringle_kernel; 96 96 } 97 97 -
sasmodels/models/rectangular_prism.c
r74768cb r8de1477 28 28 29 29 double outer_sum = 0.0; 30 for(int i=0; i< GAUSS_N; i++) {31 const double theta = 0.5 * ( G AUSS_Z[i]*(v1b-v1a) + v1a + v1b );30 for(int i=0; i<76; i++) { 31 const double theta = 0.5 * ( Gauss76Z[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< GAUSS_N; j++) {39 double phi = 0.5 * ( G AUSS_Z[j]*(v2b-v2a) + v2a + v2b );38 for(int j=0; j<76; j++) { 39 double phi = 0.5 * ( Gauss76Z[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 AUSS_W[j] * AP * AP;47 inner_sum += Gauss76Wt[j] * AP * AP; 48 48 } 49 49 inner_sum = 0.5 * (v2b-v2a) * inner_sum; 50 outer_sum += G AUSS_W[i] * inner_sum * sin_theta;50 outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 51 51 } 52 52 -
sasmodels/models/sc_paracrystal.c
r74768cb rf728001 54 54 55 55 double outer_sum = 0.0; 56 for(int i=0; i< GAUSS_N; i++) {56 for(int i=0; i<150; i++) { 57 57 double inner_sum = 0.0; 58 const double theta = G AUSS_Z[i]*theta_m + theta_b;58 const double theta = Gauss150Z[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< GAUSS_N;j++) {64 const double phi = G AUSS_Z[j]*phi_m + phi_b;63 for(int j=0;j<150;j++) { 64 const double phi = Gauss150Z[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 AUSS_W[j] * form;70 inner_sum += Gauss150Wt[j] * form; 71 71 } 72 72 inner_sum *= phi_m; // sum(f(x)dx) = sum(f(x)) dx 73 outer_sum += G AUSS_W[i] * inner_sum * sin_theta;73 outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 74 74 } 75 75 outer_sum *= theta_m; -
sasmodels/models/stacked_disks.c
r74768cb rbecded3 81 81 double halfheight = 0.5*thick_core; 82 82 83 for(int i=0; i< GAUSS_N; i++) {84 double zi = (G AUSS_Z[i] + 1.0)*M_PI_4;83 for(int i=0; i<N_POINTS_76; i++) { 84 double zi = (Gauss76Z[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 AUSS_W[i] * yyy * sin_alpha;97 summ += Gauss76Wt[i] * yyy * sin_alpha; 98 98 } 99 99 -
sasmodels/models/triaxial_ellipsoid.c
r74768cb rbecded3 21 21 const double zb = M_PI_4; 22 22 double outer = 0.0; 23 for (int i=0;i< GAUSS_N;i++) {24 //const double u = G AUSS_Z[i]*(upper-lower)/2 + (upper + lower)/2;25 const double phi = G AUSS_Z[i]*zm + zb;23 for (int i=0;i<76;i++) { 24 //const double u = Gauss76Z[i]*(upper-lower)/2 + (upper + lower)/2; 25 const double phi = Gauss76Z[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< GAUSS_N;j++) {31 for (int j=0;j<76;j++) { 32 32 // translate a point in [-1,1] to a point in [0, 1] 33 const double usq = square(G AUSS_Z[j]*um + ub);33 const double usq = square(Gauss76Z[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 AUSS_W[j] * fq * fq;36 inner += Gauss76Wt[j] * fq * fq; 37 37 } 38 outer += G AUSS_W[i] * inner; // correcting for dx later38 outer += Gauss76Wt[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.