Changes in / [ddfdb16:2694cb8] in sasmodels
- Files:
-
- 2 added
- 47 edited
Legend:
- Unmodified
- Added
- Removed
-
.gitignore
r2badeca r67cc0ff 21 21 /.idea 22 22 .vscode 23 xlate.c24 23 /sasmodels.egg-info/ 25 24 /example/Fit_*/ 26 25 /example/batch_fit.csv 26 /sasmodels/models/lib/gauss*.c -
doc/guide/plugin.rst
r0a9fcab rc654160 292 292 **Note: The order of the parameters in the definition will be the order of the 293 293 parameters in the user interface and the order of the parameters in Iq(), 294 Iqxy() and form_volume(). And** *scale* **and** *background* **parameters are 295 implicit to all models, so they do not need to be included in the parameter table.** 294 Iqac(), Iqabc() and form_volume(). And** *scale* **and** *background* 295 **parameters are implicit to all models, so they do not need to be included 296 in the parameter table.** 296 297 297 298 - **"name"** is the name of the parameter shown on the FitPage. … … 362 363 scattered intensity. 363 364 364 - "volume" parameters are passed to Iq(), Iqxy(), and form_volume(), and 365 have polydispersity loops generated automatically. 366 367 - "orientation" parameters are only passed to Iqxy(), and have angular 368 dispersion. 365 - "volume" parameters are passed to Iq(), Iqac(), Iqabc() and form_volume(), 366 and have polydispersity loops generated automatically. 367 368 - "orientation" parameters are not passed, but instead are combined with 369 orientation dispersity to translate *qx* and *qy* to *qa*, *qb* and *qc*. 370 These parameters should appear at the end of the table with the specific 371 names *theta*, *phi* and for asymmetric shapes *psi*, in that order. 369 372 370 373 Some models will have integer parameters, such as number of pearls in the … … 419 422 That is, the individual models do not need to include polydispersity 420 423 calculations, but instead rely on numerical integration to compute the 421 appropriately smeared pattern. Angular dispersion values over polar angle 422 $\theta$ requires an additional $\cos \theta$ weighting due to decreased 423 arc length for the equatorial angle $\phi$ with increasing latitude. 424 appropriately smeared pattern. 424 425 425 426 Python Models … … 468 469 barbell). If I(q; pars) is NaN for any $q$, then those parameters will be 469 470 ignored, and not included in the calculation of the weighted polydispersity. 470 471 Similar to *Iq*, you can define *Iqxy(qx, qy, par1, par2, ...)* where the472 parameter list includes any orientation parameters. If *Iqxy* is not defined,473 then it will default to *Iqxy = Iq(sqrt(qx**2+qy**2), par1, par2, ...)*.474 471 475 472 Models should define *form_volume(par1, par2, ...)* where the parameter … … 497 494 } 498 495 499 *Iqxy* is similar to *Iq*, except it uses parameters *qx, qy* instead of *q*,500 and it includes orientation parameters.501 502 496 *form_volume* defines the volume of the shape. As in python models, it 503 497 includes only the volume parameters. 504 498 505 *Iqxy* will default to *Iq(sqrt(qx**2 + qy**2), par1, ...)* and506 *form_volume* will default to 1.0.507 508 499 **source=['fn.c', ...]** includes the listed C source files in the 509 program before *Iq* and *Iqxy* are defined. This allows you to extend the 510 library of C functions available to your model. 500 program before *Iq* and *form_volume* are defined. This allows you to 501 extend the library of C functions available to your model. 502 503 *c_code* includes arbitrary C code into your kernel, which can be 504 handy for defining helper functions for *Iq* and *form_volume*. Note that 505 you can put the full function definition for *Iq* and *form_volume* 506 (include function declaration) into *c_code* as well, or put them into an 507 external C file and add that file to the list of sources. 511 508 512 509 Models are defined using double precision declarations for the … … 532 529 533 530 #define INVALID(v) (v.bell_radius < v.radius) 531 532 The INVALID define can go into *Iq*, or *c_code*, or an external C file 533 listed in *source*. 534 535 Oriented Shapes 536 ............... 537 538 If the scattering is dependent on the orientation of the shape, then you 539 will need to include *orientation* parameters *theta*, *phi* and *psi* 540 at the end of the parameter table. Shape orientation uses *a*, *b* and *c* 541 axes, corresponding to the *x*, *y* and *z* axes in the laboratory coordinate 542 system, with *z* along the beam and *x*-*y* in the detector plane, with *x* 543 horizontal and *y* vertical. The *psi* parameter rotates the shape 544 about its *c* axis, the *theta* parameter then rotates the *c* axis toward 545 the *x* axis of the detector, then *phi* rotates the shape in the detector 546 plane. (Prior to these rotations, orientation dispersity will be applied 547 as roll-pitch-yaw, rotating *c*, then *b* then *a* in the shape coordinate 548 system.) A particular *qx*, *qy* point on the detector, then corresponds 549 to *qa*, *qb*, *qc* with respect to the shape. 550 551 The oriented C model is called as *Iqabc(qa, qb, qc, par1, par2, ...)* where 552 *par1*, etc. are the parameters to the model. If the shape is rotationally 553 symmetric about *c* then *psi* is not needed, and the model is called 554 as *Iqac(qab, qc, par1, par2, ...)*. In either case, the orientation 555 parameters are not included in the function call. 556 557 For 1D oriented shapes, an integral over all angles is usually needed for 558 the *Iq* function. Given symmetry and the substitution $u = \cos(\alpha)$, 559 $du = -\sin(\alpha)\,d\alpha$ this becomes 560 561 .. math:: 562 563 I(q) &= \frac{1}{4\pi} \int_{-\pi/2}^{pi/2} \int_{-pi}^{pi} 564 F(q_a, q_b, q_c)^2 \sin(\alpha)\,d\beta\,d\alpha \\ 565 &= \frac{8}{4\pi} \int_{0}^{pi/2} \int_{0}^{\pi/2} 566 F^2 \sin(\alpha)\,d\beta\,d\alpha \\ 567 &= \frac{8}{4\pi} \int_1^0 \int_{0}^{\pi/2} - F^2 \,d\beta\,du \\ 568 &= \frac{8}{4\pi} \int_0^1 \int_{0}^{\pi/2} F^2 \,d\beta\,du 569 570 for 571 572 .. math:: 573 574 q_a &= q \sin(\alpha)\sin(\beta) = q \sqrt{1-u^2} \sin(\beta) \\ 575 q_b &= q \sin(\alpha)\cos(\beta) = q \sqrt{1-u^2} \cos(\beta) \\ 576 q_c &= q \cos(\alpha) = q u 577 578 Using the $z, w$ values for Gauss-Legendre integration in "lib/gauss76.c", the 579 numerical integration is then:: 580 581 double outer_sum = 0.0; 582 for (int i = 0; i < GAUSS_N; i++) { 583 const double cos_alpha = 0.5*GAUSS_Z[i] + 0.5; 584 const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 585 const double qc = cos_alpha * q; 586 double inner_sum = 0.0; 587 for (int j = 0; j < GAUSS_N; j++) { 588 const double beta = M_PI_4 * GAUSS_Z[j] + M_PI_4; 589 double sin_beta, cos_beta; 590 SINCOS(beta, sin_beta, cos_beta); 591 const double qa = sin_alpha * sin_beta * q; 592 const double qb = sin_alpha * cos_beta * q; 593 const double form = Fq(qa, qb, qc, ...); 594 inner_sum += GAUSS_W[j] * form * form; 595 } 596 outer_sum += GAUSS_W[i] * inner_sum; 597 } 598 outer_sum *= 0.25; // = 8/(4 pi) * outer_sum * (pi/2) / 4 599 600 The *z* values for the Gauss-Legendre integration extends from -1 to 1, so 601 the double sum of *w[i]w[j]* explains the factor of 4. Correcting for the 602 average *dz[i]dz[j]* gives $(1-0) \cdot (\pi/2-0) = \pi/2$. The $8/(4 \pi)$ 603 factor comes from the integral over the quadrant. With less symmetry (eg., 604 in the bcc and fcc paracrystal models), then an integral over the entire 605 sphere may be necessary. 606 607 For simpler models which are rotationally symmetric a single integral 608 suffices: 609 610 .. math:: 611 612 I(q) &= \frac{1}{\pi}\int_{-\pi/2}^{\pi/2} 613 F(q_{ab}, q_c)^2 \sin(\alpha)\,d\alpha/\pi \\ 614 &= \frac{2}{\pi} \int_0^1 F^2\,du 615 616 for 617 618 .. math:: 619 620 q_{ab} &= q \sin(\alpha) = q \sqrt{1 - u^2} \\ 621 q_c &= q \cos(\alpha) = q u 622 623 624 with integration loop:: 625 626 double sum = 0.0; 627 for (int i = 0; i < GAUSS_N; i++) { 628 const double cos_alpha = 0.5*GAUSS_Z[i] + 0.5; 629 const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 630 const double qab = sin_alpha * q; 631 const double qc = cos_alpha * q; 632 const double form = Fq(qab, qc, ...); 633 sum += GAUSS_W[j] * form * form; 634 } 635 sum *= 0.5; // = 2/pi * sum * (pi/2) / 2 636 637 Magnetism 638 ......... 639 640 Magnetism is supported automatically for all shapes by modifying the 641 effective SLD of particle according to the Halpern-Johnson vector 642 describing the interaction between neutron spin and magnetic field. All 643 parameters marked as type *sld* in the parameter table are treated as 644 possibly magnetic particles with magnitude *M0* and direction 645 *mtheta* and *mphi*. Polarization parameters are also provided 646 automatically for magnetic models to set the spin state of the measurement. 647 648 For more complicated systems where magnetism is not uniform throughout 649 the individual particles, you will need to write your own models. 650 You should not mark the nuclear sld as type *sld*, but instead leave 651 them unmarked and provide your own magnetism and polarization parameters. 652 For 2D measurements you will need $(q_x, q_y)$ values for the measurement 653 to compute the proper magnetism and orientation, which you can implement 654 using *Iqxy(qx, qy, par1, par2, ...)*. 534 655 535 656 Special Functions … … 796 917 show a 50x improvement or more over the equivalent pure python model. 797 918 798 External C Models799 .................800 801 External C models are very much like embedded C models, except that802 *Iq*, *Iqxy* and *form_volume* are defined in an external source file803 loaded using the *source=[...]* statement. You need to supply the function804 declarations for each of these that you need instead of building them805 automatically from the parameter table.806 807 919 808 920 .. _Form_Factors: … … 1006 1118 variable name *Rg* for example because $R_g$ is the right name for the model 1007 1119 parameter then ignore the lint errors. Also, ignore *missing-docstring* 1008 for standard model functions *Iq*, *Iq xy*, etc.1120 for standard model functions *Iq*, *Iqac*, etc. 1009 1121 1010 1122 We will have delinting sessions at the SasView Code Camps, where we can -
explore/asymint.py
r1820208 ra1c32c2 86 86 a, b, c = env.mpf(a), env.mpf(b), env.mpf(c) 87 87 def Fq(qa, qb, qc): 88 siA = env.sas_sinx_x( 0.5*a*qa/2)89 siB = env.sas_sinx_x( 0.5*b*qb/2)90 siC = env.sas_sinx_x( 0.5*c*qc/2)88 siA = env.sas_sinx_x(a*qa/2) 89 siB = env.sas_sinx_x(b*qb/2) 90 siC = env.sas_sinx_x(c*qc/2) 91 91 return siA * siB * siC 92 92 Fq.__doc__ = "parallelepiped a=%g, b=%g c=%g"%(a, b, c) 93 93 volume = a*b*c 94 94 norm = CONTRAST**2*volume/10000 95 return norm, Fq 96 97 def make_core_shell_parallelepiped(a, b, c, da, db, dc, slda, sldb, sldc, env=NPenv): 98 overlapping = False 99 a, b, c = env.mpf(a), env.mpf(b), env.mpf(c) 100 da, db, dc = env.mpf(da), env.mpf(db), env.mpf(dc) 101 slda, sldb, sldc = env.mpf(slda), env.mpf(sldb), env.mpf(sldc) 102 dr0 = CONTRAST 103 drA, drB, drC = slda-SLD_SOLVENT, sldb-SLD_SOLVENT, sldc-SLD_SOLVENT 104 tA, tB, tC = a + 2*da, b + 2*db, c + 2*dc 105 def Fq(qa, qb, qc): 106 siA = a*env.sas_sinx_x(a*qa/2) 107 siB = b*env.sas_sinx_x(b*qb/2) 108 siC = c*env.sas_sinx_x(c*qc/2) 109 siAt = tA*env.sas_sinx_x(tA*qa/2) 110 siBt = tB*env.sas_sinx_x(tB*qb/2) 111 siCt = tC*env.sas_sinx_x(tC*qc/2) 112 if overlapping: 113 return (dr0*siA*siB*siC 114 + drA*(siAt-siA)*siB*siC 115 + drB*siAt*(siBt-siB)*siC 116 + drC*siAt*siBt*(siCt-siC)) 117 else: 118 return (dr0*siA*siB*siC 119 + drA*(siAt-siA)*siB*siC 120 + drB*siA*(siBt-siB)*siC 121 + drC*siA*siB*(siCt-siC)) 122 Fq.__doc__ = "core-shell parallelepiped a=%g, b=%g c=%g"%(a, b, c) 123 if overlapping: 124 volume = a*b*c + 2*da*b*c + 2*tA*db*c + 2*tA*tB*dc 125 else: 126 volume = a*b*c + 2*da*b*c + 2*a*db*c + 2*a*b*dc 127 norm = 1/(volume*10000) 95 128 return norm, Fq 96 129 … … 184 217 NORM, KERNEL = make_parallelepiped(A, B, C) 185 218 NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 219 elif shape == 'core_shell_parallelepiped': 220 #A, B, C = 4450, 14000, 47 221 #A, B, C = 445, 140, 47 # integer for the sake of mpf 222 A, B, C = 6800, 114, 1380 223 DA, DB, DC = 2300, 21, 58 224 SLDA, SLDB, SLDC = "5", "-0.3", "11.5" 225 #A,B,C,DA,DB,DC,SLDA,SLDB,SLDC = 10,20,30,100,200,300,1,2,3 226 #SLD_SOLVENT,CONTRAST = 0, 4 227 if 1: # C shortest 228 B, C = C, B 229 DB, DC = DC, DB 230 SLDB, SLDC = SLDC, SLDB 231 elif 0: # C longest 232 A, C = C, A 233 DA, DC = DC, DA 234 SLDA, SLDC = SLDC, SLDA 235 NORM, KERNEL = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC) 236 NORM_MP, KERNEL_MP = make_core_shell_parallelepiped(A, B, C, DA, DB, DC, SLDA, SLDB, SLDC, env=MPenv) 186 237 elif shape == 'paracrystal': 187 238 LATTICE = 'bcc' … … 342 393 print("gauss-150", *gauss_quad_2d(Q, n=150)) 343 394 print("gauss-500", *gauss_quad_2d(Q, n=500)) 395 print("gauss-1025", *gauss_quad_2d(Q, n=1025)) 396 print("gauss-2049", *gauss_quad_2d(Q, n=2049)) 344 397 #gridded_2d(Q, n=2**8+1) 345 398 gridded_2d(Q, n=2**10+1) 346 #gridded_2d(Q, n=2**1 3+1)399 #gridded_2d(Q, n=2**12+1) 347 400 #gridded_2d(Q, n=2**15+1) 348 if shape != 'paracrystal': # adaptive forms are too slow! 401 if shape not in ('paracrystal', 'core_shell_parallelepiped'): 402 # adaptive forms on models for which the calculations are fast enough 349 403 print("dblquad", *scipy_dblquad_2d(Q)) 350 404 print("semi-romberg-100", *semi_romberg_2d(Q, n=100)) -
sasmodels/autoc.py
r2badeca r67cc0ff 49 49 return 50 50 51 public_methods = "Iq", "Iq xy", "form_volume"51 public_methods = "Iq", "Iqac", "Iqabc", "Iqxy", "form_volume" 52 52 53 53 tagged = [] # type: List[str] … … 72 72 function_name, function = translate.pop(0) 73 73 filename = function.__code__.co_filename 74 escaped_filename = filename.replace('\\', '\\\\') 74 75 offset = function.__code__.co_firstlineno 75 76 refs = function.__code__.co_names … … 91 92 # not special: add function to translate stack 92 93 translate.append((name, obj)) 93 elif isinstance(obj, float):94 elif isinstance(obj, (int, float, list, tuple, np.ndarray)): 94 95 constants[name] = obj 95 snippets.append("const double %s = %.15g;"%(name, obj)) 96 elif isinstance(obj, int): 97 constants[name] = obj 98 snippets.append("const int %s = %d;"%(name, obj)) 99 elif isinstance(obj, (list, tuple, np.ndarray)): 100 constants[name] = obj 101 # extend constant arrays to a multiple of 4; not sure if this 102 # is necessary, but some OpenCL targets broke if the number 103 # of parameters in the parameter table was not a multiple of 4, 104 # so do it for all constant arrays to be safe. 105 if len(obj)%4 != 0: 106 obj = list(obj) + [0.]*(4-len(obj)) 107 vals = ", ".join("%.15g"%v for v in obj) 108 snippets.append("const double %s[] = {%s};" %(name, vals)) 96 # Claim all constants are declared on line 1 97 snippets.append('#line 1 "%s"'%escaped_filename) 98 snippets.append(define_constant(name, obj)) 109 99 elif isinstance(obj, special.Gauss): 110 constants["GAUSS_N"] = obj.n 111 constants["GAUSS_Z"] = obj.z 112 constants["GAUSS_W"] = obj.w 113 libs.append('lib/gauss%d.c'%obj.n) 100 for var, value in zip(("N", "Z", "W"), (obj.n, obj.z, obj.w)): 101 var = "GAUSS_"+var 102 constants[var] = value 103 snippets.append('#line 1 "%s"'%escaped_filename) 104 snippets.append(define_constant(var, value)) 105 #libs.append('lib/gauss%d.c'%obj.n) 114 106 source = (source.replace(name+'.n', 'GAUSS_N') 115 107 .replace(name+'.z', 'GAUSS_Z') … … 129 121 130 122 # translate source 131 functions = py2c.translate( 132 [code[name] for name in ordered_dag(depends) if name in code], 133 constants) 134 snippets.clear() 135 snippets.append(functions) 136 #print("source", info.source) 137 print("\n".join(snippets)) 138 try: 139 c_text = "\n".join(snippets) 140 translated = open ("_autoc.c", "a+") 141 translated.write (c_text) 142 translated.close() 143 except Exception as excp: 144 strErr = "Error:\n" + str(excp.args) 145 print(strErr) 146 #return 147 # raise RuntimeError("not yet converted...") 123 ordered_code = [code[name] for name in ordered_dag(depends) if name in code] 124 functions = py2c.translate(ordered_code, constants) 125 snippets.extend(functions) 148 126 149 127 # update model info 150 128 info.source = unique_libs 151 129 info.c_code = "\n".join(snippets) 152 info.Iq = info.Iqxy = info.form_volume = None 130 info.Iq = info.Iqac = info.Iqabc = info.Iqxy = info.form_volume = None 131 132 def define_constant(name, value): 133 if isinstance(value, int): 134 parts = ["int ", name, " = ", "%d"%value, ";"] 135 elif isinstance(value, float): 136 parts = ["double ", name, " = ", "%.15g"%value, ";"] 137 else: 138 # extend constant arrays to a multiple of 4; not sure if this 139 # is necessary, but some OpenCL targets broke if the number 140 # of parameters in the parameter table was not a multiple of 4, 141 # so do it for all constant arrays to be safe. 142 if len(value)%4 != 0: 143 value = list(value) + [0.]*(4 - len(value)%4) 144 elements = ["%.15g"%v for v in value] 145 parts = ["double ", name, "[]", " = ", 146 "{\n ", ", ".join(elements), "\n};"] 147 return "".join(parts) 153 148 154 149 -
sasmodels/compare.py
r2d81cfe r2d81cfe 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 … … 693 693 data = empty_data2D(q, resolution=res) 694 694 data.accuracy = opts['accuracy'] 695 set_beam_stop(data, 0.0004)695 set_beam_stop(data, qmin) 696 696 index = ~data.mask 697 697 else: … … 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 … … 1274 1288 if model_info != model_info2: 1275 1289 pars2 = randomize_pars(model_info2, pars2) 1276 limit_dimensions(model_info , pars2, maxdim)1290 limit_dimensions(model_info2, pars2, maxdim) 1277 1291 # Share values for parameters with the same name 1278 1292 for k, v in pars.items(): -
sasmodels/details.py
r2d81cfe r108e70e 258 258 # type: (...) -> Sequence[np.ndarray] 259 259 """ 260 **Deprecated** Theta weights will be computed in the kernel wrapper if 261 they are needed. 262 260 263 If there is a theta parameter, update the weights of that parameter so that 261 264 the cosine weighting required for polar integration is preserved. … … 272 275 Returns updated weights vectors 273 276 """ 274 # TODO: explain in a comment why scale and background are missing275 277 # Apparently the parameters.theta_offset similarly skips scale and 276 278 # and background, so the indexing works out, but they are still shipped … … 279 281 index = parameters.theta_offset 280 282 theta = dispersity[index] 281 # TODO: modify the dispersity vector to avoid the theta=-90,90,270,...282 283 theta_weight = abs(cos(radians(theta))) 283 284 weights = tuple(theta_weight*w if k == index else w -
sasmodels/generate.py
rdb03406 r108e70e 7 7 particular dimensions averaged over all orientations. 8 8 9 *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy for a form 10 with particular dimensions for a single orientation. 11 12 *Imagnetic(qx, qy, result[], p1, p2, ...)* returns the scattering for the 13 polarized neutron spin states (up-up, up-down, down-up, down-down) for 14 a form with particular dimensions for a single orientation. 9 *Iqac(qab, qc, p1, p2, ...)* returns the scattering at qab, qc 10 for a rotationally symmetric form with particular dimensions. 11 qab, qc are determined from shape orientation and scattering angles. 12 This call is used if the shape has orientation parameters theta and phi. 13 14 *Iqabc(qa, qb, qc, p1, p2, ...)* returns the scattering at qa, qb, qc 15 for a form with particular dimensions. qa, qb, qc are determined from 16 shape orientation and scattering angles. This call is used if the shape 17 has orientation parameters theta, phi and psi. 18 19 *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy. Use this 20 to create an arbitrary 2D theory function, needed for q-dependent 21 background functions and for models with non-uniform magnetism. 15 22 16 23 *form_volume(p1, p2, ...)* returns the volume of the form with particular … … 31 38 scale and background parameters for each model. 32 39 33 *Iq*, *Iqxy*, *Imagnetic* and *form_volume* should be stylized C-99 34 functions written for OpenCL. All functions need prototype declarations 35 even if the are defined before they are used. OpenCL does not support 36 *#include* preprocessor directives, so instead the list of includes needs 37 to be given as part of the metadata in the kernel module definition. 38 The included files should be listed using a path relative to the kernel 39 module, or if using "lib/file.c" if it is one of the standard includes 40 provided with the sasmodels source. The includes need to be listed in 41 order so that functions are defined before they are used. 40 C code should be stylized C-99 functions written for OpenCL. All functions 41 need prototype declarations even if the are defined before they are used. 42 Although OpenCL supports *#include* preprocessor directives, the list of 43 includes should be given as part of the metadata in the kernel module 44 definition. The included files should be listed using a path relative to the 45 kernel module, or if using "lib/file.c" if it is one of the standard includes 46 provided with the sasmodels source. The includes need to be listed in order 47 so that functions are defined before they are used. 42 48 43 49 Floating point values should be declared as *double*. For single precision … … 107 113 present, the volume ratio is 1. 108 114 109 *form_volume*, *Iq*, *Iq xy*, *Imagnetic* are strings containing the110 C source code for the body of the volume, Iq, and Iqxyfunctions115 *form_volume*, *Iq*, *Iqac*, *Iqabc* are strings containing 116 the C source code for the body of the volume, Iq, and Iqac functions 111 117 respectively. These can also be defined in the last source file. 112 118 113 *Iq* and *Iqxy* also be instead be python functions defining the119 *Iq*, *Iqac*, *Iqabc* also be instead be python functions defining the 114 120 kernel. If they are marked as *Iq.vectorized = True* then the 115 121 kernel is passed the entire *q* vector at once, otherwise it is … … 168 174 from zlib import crc32 169 175 from inspect import currentframe, getframeinfo 176 import logging 170 177 171 178 import numpy as np # type: ignore … … 181 188 pass 182 189 # pylint: enable=unused-import 190 191 logger = logging.getLogger(__name__) 183 192 184 193 # jitter projection to use in the kernel code. See explore/jitter.py … … 270 279 """ 271 280 281 282 def set_integration_size(info, n): 283 # type: (ModelInfo, int) -> None 284 """ 285 Update the model definition, replacing the gaussian integration with 286 a gaussian integration of a different size. 287 288 Note: this really ought to be a method in modelinfo, but that leads to 289 import loops. 290 """ 291 if (info.source and any(lib.startswith('lib/gauss') for lib in info.source)): 292 import os.path 293 from .gengauss import gengauss 294 path = os.path.join(MODEL_PATH, "lib", "gauss%d.c"%n) 295 if not os.path.exists(path): 296 gengauss(n, path) 297 info.source = ["lib/gauss%d.c"%n if lib.startswith('lib/gauss') 298 else lib for lib in info.source] 272 299 273 300 def format_units(units): … … 608 635 609 636 """ 610 def _gen_fn( name, pars, body, filename, line):611 # type: ( str, List[Parameter], str, str, int) -> str637 def _gen_fn(model_info, name, pars): 638 # type: (ModelInfo, str, List[Parameter]) -> str 612 639 """ 613 640 Generate a function given pars and body. … … 621 648 """ 622 649 par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void' 650 body = getattr(model_info, name) 651 filename = model_info.filename 652 # Note: if symbol is defined strangely in the module then default it to 1 653 lineno = model_info.lineno.get(name, 1) 623 654 return _FN_TEMPLATE % { 624 655 'name': name, 'pars': par_decl, 'body': body, 625 'filename': filename.replace('\\', '\\\\'), 'line': line ,656 'filename': filename.replace('\\', '\\\\'), 'line': lineno, 626 657 } 627 658 … … 638 669 639 670 # type in IQXY pattern could be single, float, double, long double, ... 640 _IQXY_PATTERN = re.compile( "^((inline|static) )? *([a-z ]+ )? *Iqxy *([(]|$)",671 _IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]", 641 672 flags=re.MULTILINE) 642 def _have_Iqxy(sources):673 def find_xy_mode(source): 643 674 # type: (List[str]) -> bool 644 675 """ 645 Return t rue if any file defines Iqxy.676 Return the xy mode as qa, qac, qabc or qxy. 646 677 647 678 Note this is not a C parser, and so can be easily confused by 648 679 non-standard syntax. Also, it will incorrectly identify the following 649 as having Iqxy::680 as having 2D models:: 650 681 651 682 /* 652 double Iq xy(qx, qy, ...) { ... fill this in later ... }683 double Iqac(qab, qc, ...) { ... fill this in later ... } 653 684 */ 654 685 655 If you want to comment out an Iqxy function, use // on the front of the 656 line instead. 657 """ 658 for _path, code in sources: 659 if _IQXY_PATTERN.search(code): 660 return True 661 return False 662 663 664 def _add_source(source, code, path): 686 If you want to comment out the function, use // on the front of the 687 line:: 688 689 /* 690 // double Iqac(qab, qc, ...) { ... fill this in later ... } 691 */ 692 693 """ 694 for code in source: 695 m = _IQXY_PATTERN.search(code) 696 if m is not None: 697 return m.group('mode') 698 return 'qa' 699 700 701 def _add_source(source, code, path, lineno=1): 665 702 """ 666 703 Add a file to the list of source code chunks, tagged with path and line. 667 704 """ 668 705 path = path.replace('\\', '\\\\') 669 source.append('#line 1 "%s"' % path)706 source.append('#line %d "%s"' % (lineno, path)) 670 707 source.append(code) 671 708 … … 698 735 user_code = [(f, open(f).read()) for f in model_sources(model_info)] 699 736 700 # What kind of 2D model do we need?701 xy_mode = ('qa' if not _have_Iqxy(user_code) and not isinstance(model_info.Iqxy, str)702 else 'qac' if not partable.is_asymmetric703 else 'qabc')704 705 737 # Build initial sources 706 738 source = [] … … 710 742 711 743 if model_info.c_code: 712 source.append(model_info.c_code) 744 _add_source(source, model_info.c_code, model_info.filename, 745 lineno=model_info.lineno.get('c_code', 1)) 713 746 714 747 # Make parameters for q, qx, qy so that we can use them in declarations 715 q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 748 q, qx, qy, qab, qa, qb, qc \ 749 = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()] 716 750 # Generate form_volume function, etc. from body only 717 751 if isinstance(model_info.form_volume, str): 718 752 pars = partable.form_volume_parameters 719 source.append(_gen_fn('form_volume', pars, model_info.form_volume, 720 model_info.filename, model_info._form_volume_line)) 753 source.append(_gen_fn(model_info, 'form_volume', pars)) 721 754 if isinstance(model_info.Iq, str): 722 755 pars = [q] + partable.iq_parameters 723 source.append(_gen_fn('Iq', pars, model_info.Iq, 724 model_info.filename, model_info._Iq_line)) 756 source.append(_gen_fn(model_info, 'Iq', pars)) 725 757 if isinstance(model_info.Iqxy, str): 726 pars = [qx, qy] + partable.iqxy_parameters 727 source.append(_gen_fn('Iqxy', pars, model_info.Iqxy, 728 model_info.filename, model_info._Iqxy_line)) 758 pars = [qx, qy] + partable.iq_parameters + partable.orientation_parameters 759 source.append(_gen_fn(model_info, 'Iqxy', pars)) 760 if isinstance(model_info.Iqac, str): 761 pars = [qab, qc] + partable.iq_parameters 762 source.append(_gen_fn(model_info, 'Iqac', pars)) 763 if isinstance(model_info.Iqabc, str): 764 pars = [qa, qb, qc] + partable.iq_parameters 765 source.append(_gen_fn(model_info, 'Iqabc', pars)) 766 767 # What kind of 2D model do we need? Is it consistent with the parameters? 768 xy_mode = find_xy_mode(source) 769 if xy_mode == 'qabc' and not partable.is_asymmetric: 770 raise ValueError("asymmetric oriented models need to define Iqabc") 771 elif xy_mode == 'qac' and partable.is_asymmetric: 772 raise ValueError("symmetric oriented models need to define Iqac") 773 elif not partable.orientation_parameters and xy_mode in ('qac', 'qabc'): 774 raise ValueError("Unexpected function I%s for unoriented shape"%xy_mode) 775 elif partable.orientation_parameters and xy_mode not in ('qac', 'qabc'): 776 if xy_mode == 'qxy': 777 logger.warn("oriented shapes should define Iqac or Iqabc") 778 else: 779 raise ValueError("Expected function Iqac or Iqabc for oriented shape") 729 780 730 781 # Define the parameter table … … 752 803 if xy_mode == 'qabc': 753 804 pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 754 call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iq xy(%s)" % pars805 call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqabc(%s)" % pars 755 806 clear_iqxy = "#undef CALL_IQ_ABC" 756 807 elif xy_mode == 'qac': 757 808 pars = ",".join(["_qa", "_qc"] + model_refs) 758 call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iq xy(%s)" % pars809 call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars 759 810 clear_iqxy = "#undef CALL_IQ_AC" 760 el se: # xy_mode == 'qa'811 elif xy_mode == 'qa': 761 812 pars = ",".join(["_qa"] + model_refs) 762 813 call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 763 814 clear_iqxy = "#undef CALL_IQ_A" 815 elif xy_mode == 'qxy': 816 orientation_refs = _call_pars("_v.", partable.orientation_parameters) 817 pars = ",".join(["_qx", "_qy"] + model_refs + orientation_refs) 818 call_iqxy = "#define CALL_IQ_XY(_qx,_qy,_v) Iqxy(%s)" % pars 819 clear_iqxy = "#undef CALL_IQ_XY" 820 if partable.orientation_parameters: 821 call_iqxy += "\n#define HAVE_THETA" 822 clear_iqxy += "\n#undef HAVE_THETA" 823 if partable.is_asymmetric: 824 call_iqxy += "\n#define HAVE_PSI" 825 clear_iqxy += "\n#undef HAVE_PSI" 826 764 827 765 828 magpars = [k-2 for k, p in enumerate(partable.call_parameters) -
sasmodels/kernel_header.c
r8698a0d r108e70e 150 150 inline double cube(double x) { return x*x*x; } 151 151 inline double sas_sinx_x(double x) { return x==0 ? 1.0 : sin(x)/x; } 152 153 // CRUFT: support old style models with orientation received qx, qy and angles 154 155 // To rotate from the canonical position to theta, phi, psi, first rotate by 156 // psi about the major axis, oriented along z, which is a rotation in the 157 // detector plane xy. Next rotate by theta about the y axis, aligning the major 158 // axis in the xz plane. Finally, rotate by phi in the detector plane xy. 159 // To compute the scattering, undo these rotations in reverse order: 160 // rotate in xy by -phi, rotate in xz by -theta, rotate in xy by -psi 161 // The returned q is the length of the q vector and (xhat, yhat, zhat) is a unit 162 // vector in the q direction. 163 // To change between counterclockwise and clockwise rotation, change the 164 // sign of phi and psi. 165 166 #if 1 167 //think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 168 #define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 169 SINCOS(phi*M_PI_180, sn, cn); \ 170 q = sqrt(qx*qx + qy*qy); \ 171 cn = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180)); \ 172 sn = sqrt(1 - cn*cn); \ 173 } while (0) 174 #else 175 // SasView 3.x definition of orientation 176 #define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 177 SINCOS(theta*M_PI_180, sn, cn); \ 178 q = sqrt(qx*qx + qy*qy);\ 179 cn = (q==0. ? 1.0 : (cn*cos(phi*M_PI_180)*qx + sn*qy)/q); \ 180 sn = sqrt(1 - cn*cn); \ 181 } while (0) 182 #endif 183 184 #if 1 185 #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, xhat, yhat, zhat) do { \ 186 q = sqrt(qx*qx + qy*qy); \ 187 const double qxhat = qx/q; \ 188 const double qyhat = qy/q; \ 189 double sin_theta, cos_theta; \ 190 double sin_phi, cos_phi; \ 191 double sin_psi, cos_psi; \ 192 SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 193 SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 194 SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 195 xhat = qxhat*(-sin_phi*sin_psi + cos_theta*cos_phi*cos_psi) \ 196 + qyhat*( cos_phi*sin_psi + cos_theta*sin_phi*cos_psi); \ 197 yhat = qxhat*(-sin_phi*cos_psi - cos_theta*cos_phi*sin_psi) \ 198 + qyhat*( cos_phi*cos_psi - cos_theta*sin_phi*sin_psi); \ 199 zhat = qxhat*(-sin_theta*cos_phi) \ 200 + qyhat*(-sin_theta*sin_phi); \ 201 } while (0) 202 #else 203 // SasView 3.x definition of orientation 204 #define ORIENT_ASYMMETRIC(qx, qy, theta, phi, psi, q, cos_alpha, cos_mu, cos_nu) do { \ 205 q = sqrt(qx*qx + qy*qy); \ 206 const double qxhat = qx/q; \ 207 const double qyhat = qy/q; \ 208 double sin_theta, cos_theta; \ 209 double sin_phi, cos_phi; \ 210 double sin_psi, cos_psi; \ 211 SINCOS(theta*M_PI_180, sin_theta, cos_theta); \ 212 SINCOS(phi*M_PI_180, sin_phi, cos_phi); \ 213 SINCOS(psi*M_PI_180, sin_psi, cos_psi); \ 214 cos_alpha = cos_theta*cos_phi*qxhat + sin_theta*qyhat; \ 215 cos_mu = (-sin_theta*cos_psi*cos_phi - sin_psi*sin_phi)*qxhat + cos_theta*cos_psi*qyhat; \ 216 cos_nu = (-cos_phi*sin_psi*sin_theta + sin_phi*cos_psi)*qxhat + sin_psi*cos_theta*qyhat; \ 217 } while (0) 218 #endif -
sasmodels/kernel_iq.c
r6aee3ab r108e70e 31 31 // CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes 32 32 // CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes 33 // CALL_IQ_XY(qx, qy, table) : call the Iqxy function for arbitrary models 33 34 // INVALID(table) : test if the current point is feesible to calculate. This 34 35 // will be defined in the kernel definition file. … … 469 470 #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc) 470 471 #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 472 #elif defined(CALL_IQ_XY) 473 // direct call to qx,qy calculator 474 double qx, qy; 475 #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 476 #define BUILD_ROTATION() do {} while(0) 477 #define APPLY_ROTATION() do {} while(0) 478 #define CALL_KERNEL() CALL_IQ_XY(qx, qy, local_values.table) 471 479 #endif 472 480 … … 477 485 #if defined(CALL_IQ) || defined(CALL_IQ_A) 478 486 #define APPLY_PROJECTION() const double weight=weight0 487 #elif defined(CALL_IQ_XY) 488 // CRUFT: support oriented model which define Iqxy rather than Iqac or Iqabc 489 // Need to plug the values for the orientation angles back into parameter 490 // table in case they were overridden by the orientation offset. This 491 // means that orientation dispersity will not work for these models, but 492 // it was broken anyway, so no matter. Still want to provide Iqxy in case 493 // the user model wants full control of orientation/magnetism. 494 #if defined(HAVE_PSI) 495 const double theta = values[details->theta_par+2]; 496 const double phi = values[details->theta_par+3]; 497 const double psi = values[details->theta_par+4]; 498 double weight; 499 #define APPLY_PROJECTION() do { \ 500 local_values.table.theta = theta; \ 501 local_values.table.phi = phi; \ 502 local_values.table.psi = psi; \ 503 weight=weight0; \ 504 } while (0) 505 #elif defined(HAVE_THETA) 506 const double theta = values[details->theta_par+2]; 507 const double phi = values[details->theta_par+3]; 508 double weight; 509 #define APPLY_PROJECTION() do { \ 510 local_values.table.theta = theta; \ 511 local_values.table.phi = phi; \ 512 weight=weight0; \ 513 } while (0) 514 #else 515 #define APPLY_PROJECTION() const double weight=weight0 516 #endif 479 517 #else // !spherosymmetric projection 480 518 // Grab the "view" angles (theta, phi, psi) from the initial parameter table. -
sasmodels/kerneldll.py
r2d81cfe r1ddb794 185 185 subprocess.check_output(command, shell=shell, stderr=subprocess.STDOUT) 186 186 except subprocess.CalledProcessError as exc: 187 raise RuntimeError("compile failed.\n%s\n%s"%(command_str, exc.output)) 187 raise RuntimeError("compile failed.\n%s\n%s" 188 % (command_str, exc.output.decode())) 188 189 if not os.path.exists(output): 189 190 raise RuntimeError("compile failed. File is in %r"%source) -
sasmodels/kernelpy.py
r2d81cfe r108e70e 26 26 # pylint: enable=unused-import 27 27 28 logger = logging.getLogger(__name__) 29 28 30 class PyModel(KernelModel): 29 31 """ … … 31 33 """ 32 34 def __init__(self, model_info): 33 # Make sure Iq and Iqxy areavailable and vectorized35 # Make sure Iq is available and vectorized 34 36 _create_default_functions(model_info) 35 37 self.info = model_info 36 38 self.dtype = np.dtype('d') 37 logg ing.info("load python model " + self.info.name)39 logger.info("load python model " + self.info.name) 38 40 39 41 def make_kernel(self, q_vectors): 40 42 q_input = PyInput(q_vectors, dtype=F64) 41 kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq 42 return PyKernel(kernel, self.info, q_input) 43 return PyKernel(self.info, q_input) 43 44 44 45 def release(self): … … 89 90 Callable SAS kernel. 90 91 91 *kernel* is the DllKernel object to call.92 *kernel* is the kernel object to call. 92 93 93 94 *model_info* is the module information … … 104 105 Call :meth:`release` when done with the kernel instance. 105 106 """ 106 def __init__(self, kernel,model_info, q_input):107 def __init__(self, model_info, q_input): 107 108 # type: (callable, ModelInfo, List[np.ndarray]) -> None 108 109 self.dtype = np.dtype('d') … … 110 111 self.q_input = q_input 111 112 self.res = np.empty(q_input.nq, q_input.dtype) 112 self.kernel = kernel113 113 self.dim = '2d' if q_input.is_2d else '1d' 114 114 … … 159 159 # Generate a closure which calls the form_volume if it exists. 160 160 form_volume = model_info.form_volume 161 self._volume = ((lambda: form_volume(*volume_args)) if form_volume 162 else(lambda: 1.0))161 self._volume = ((lambda: form_volume(*volume_args)) if form_volume else 162 (lambda: 1.0)) 163 163 164 164 def __call__(self, call_details, values, cutoff, magnetic): … … 261 261 any functions that are not already marked as vectorized. 262 262 """ 263 # Note: must call create_vector_Iq before create_vector_Iqxy 263 264 _create_vector_Iq(model_info) 264 _create_vector_Iqxy(model_info) # call create_vector_Iq() first265 _create_vector_Iqxy(model_info) 265 266 266 267 … … 280 281 model_info.Iq = vector_Iq 281 282 283 282 284 def _create_vector_Iqxy(model_info): 283 285 """ 284 286 Define Iqxy as a vector function if it exists, or default it from Iq(). 285 287 """ 286 Iq , Iqxy = model_info.Iq, model_info.Iqxy288 Iqxy = getattr(model_info, 'Iqxy', None) 287 289 if callable(Iqxy): 288 290 if not getattr(Iqxy, 'vectorized', False): … … 295 297 vector_Iqxy.vectorized = True 296 298 model_info.Iqxy = vector_Iqxy 297 el if callable(Iq):299 else: 298 300 #print("defaulting Iqxy") 299 301 # Iq is vectorized because create_vector_Iq was already called. 302 Iq = model_info.Iq 300 303 def default_Iqxy(qx, qy, *args): 301 304 """ -
sasmodels/modelinfo.py
rdb03406 r67cc0ff 42 42 43 43 # assumptions about common parameters exist throughout the code, such as: 44 # (1) kernel functions Iq, Iqxy, form_volume, ... don't see them44 # (1) kernel functions Iq, Iqxy, Iqac, Iqabc, form_volume, ... don't see them 45 45 # (2) kernel drivers assume scale is par[0] and background is par[1] 46 46 # (3) mixture models drop the background on components and replace the scale … … 261 261 262 262 *type* indicates how the parameter will be used. "volume" parameters 263 will be used in all functions. "orientation" parameters will be used264 in *Iqxy* and *Imagnetic*. "magnetic* parameters will be used in265 *I magnetic* only. If *type* is the empty string, the parameter will263 will be used in all functions. "orientation" parameters are not passed, 264 but will be used to convert from *qx*, *qy* to *qa*, *qb*, *qc* in calls to 265 *Iqxy* and *Imagnetic*. If *type* is the empty string, the parameter will 266 266 be used in all of *Iq*, *Iqxy* and *Imagnetic*. "sld" parameters 267 267 can automatically be promoted to magnetic parameters, each of which … … 391 391 with vector parameter p sent as p[]. 392 392 393 * [removed] *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...)394 function, with vector parameter p sent as p[].395 396 393 * *form_volume_parameters* is the list of parameters to the form_volume(...) 397 394 function, with vector parameter p sent as p[]. … … 448 445 self.iq_parameters = [p for p in self.kernel_parameters 449 446 if p.type not in ('orientation', 'magnetic')] 450 # note: orientation no longer sent to Iqxy, so its the same as 451 #self.iqxy_parameters = [p for p in self.kernel_parameters 452 # if p.type != 'magnetic'] 447 self.orientation_parameters = [p for p in self.kernel_parameters 448 if p.type == 'orientation'] 453 449 self.form_volume_parameters = [p for p in self.kernel_parameters 454 450 if p.type == 'volume'] … … 495 491 if p.type != 'orientation': 496 492 raise TypeError("psi must be an orientation parameter") 493 elif p.type == 'orientation': 494 raise TypeError("only theta, phi and psi can be orientation parameters") 497 495 if theta >= 0 and phi >= 0: 496 last_par = len(self.kernel_parameters) - 1 498 497 if phi != theta+1: 499 498 raise TypeError("phi must follow theta") 500 499 if psi >= 0 and psi != phi+1: 501 500 raise TypeError("psi must follow phi") 501 #if (psi >= 0 and psi != last_par) or (psi < 0 and phi != last_par): 502 # raise TypeError("orientation parameters must appear at the " 503 # "end of the parameter table") 502 504 elif theta >= 0 or phi >= 0 or psi >= 0: 503 505 raise TypeError("oriented shapes must have both theta and phi and maybe psi") … … 720 722 721 723 724 #: Set of variables defined in the model that might contain C code 725 C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code'] 726 722 727 def _find_source_lines(model_info, kernel_module): 723 728 # type: (ModelInfo, ModuleType) -> None … … 725 730 Identify the location of the C source inside the model definition file. 726 731 727 This code runs through the source of the kernel module looking for 728 lines that start with 'Iq', 'Iqxy' or 'form_volume'. Clearly there are 729 all sorts of reasons why this might not work (e.g., code commented out 730 in a triple-quoted line block, code built using string concatenation, 731 or code defined in the branch of an 'if' block), but it should work 732 properly in the 95% case, and getting the incorrect line number will 733 be harmless. 734 """ 735 # Check if we need line numbers at all 736 if callable(model_info.Iq): 737 return None 738 739 if (model_info.Iq is None 740 and model_info.Iqxy is None 741 and model_info.Imagnetic is None 742 and model_info.form_volume is None): 732 This code runs through the source of the kernel module looking for lines 733 that contain C code (because it is a c function definition). Clearly 734 there are all sorts of reasons why this might not work (e.g., code 735 commented out in a triple-quoted line block, code built using string 736 concatenation, code defined in the branch of an 'if' block, code imported 737 from another file), but it should work properly in the 95% case, and for 738 the remainder, getting the incorrect line number will merely be 739 inconvenient. 740 """ 741 # Only need line numbers if we are creating a C module and the C symbols 742 # are defined. 743 if (callable(model_info.Iq) 744 or not any(hasattr(model_info, s) for s in C_SYMBOLS)): 743 745 return 744 746 745 # find the defintion lines for the different code blocks747 # load the module source if we can 746 748 try: 747 749 source = inspect.getsource(kernel_module) 748 750 except IOError: 749 751 return 750 for k, v in enumerate(source.split('\n')): 751 if v.startswith('Imagnetic'): 752 model_info._Imagnetic_line = k+1 753 elif v.startswith('Iqxy'): 754 model_info._Iqxy_line = k+1 755 elif v.startswith('Iq'): 756 model_info._Iq_line = k+1 757 elif v.startswith('form_volume'): 758 model_info._form_volume_line = k+1 759 752 753 # look for symbol at the start of the line 754 for lineno, line in enumerate(source.split('\n')): 755 for name in C_SYMBOLS: 756 if line.startswith(name): 757 # Add 1 since some compilers complain about "#line 0" 758 model_info.lineno[name] = lineno + 1 759 break 760 760 761 761 def make_model_info(kernel_module): … … 766 766 Fill in default values for parts of the module that are not provided. 767 767 768 Note: vectorized Iq and Iq xyfunctions will be created for python768 Note: vectorized Iq and Iqac/Iqabc functions will be created for python 769 769 models when the model is first called, not when the model is loaded. 770 770 """ … … 796 796 info.c_code = getattr(kernel_module, 'c_code', None) 797 797 info.source = getattr(kernel_module, 'source', []) 798 info.c_code = getattr(kernel_module, 'c_code', None) 798 799 # TODO: check the structure of the tests 799 800 info.tests = getattr(kernel_module, 'tests', []) … … 803 804 info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 804 805 info.Iqxy = getattr(kernel_module, 'Iqxy', None) # type: ignore 806 info.Iqac = getattr(kernel_module, 'Iqac', None) # type: ignore 807 info.Iqabc = getattr(kernel_module, 'Iqabc', None) # type: ignore 805 808 info.Imagnetic = getattr(kernel_module, 'Imagnetic', None) # type: ignore 806 809 info.profile = getattr(kernel_module, 'profile', None) # type: ignore … … 817 820 info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore 818 821 822 info.lineno = {} 819 823 _find_source_lines(info, kernel_module) 820 try: 821 autoc.convert(info, kernel_module) 822 except Exception as exc: 823 raise 824 logger.warn(str(exc)) 824 if getattr(kernel_module, 'py2c', False): 825 try: 826 autoc.convert(info, kernel_module) 827 except Exception as exc: 828 logger.warn(str(exc) + " while converting %s from C to python"%name) 829 830 if callable(info.Iq) and parameters.has_2d: 831 raise ValueError("oriented python models not supported") 825 832 826 833 return info … … 835 842 836 843 The structure should be mostly static, other than the delayed definition 837 of *Iq* and *Iqxy* if they need to be defined.844 of *Iq*, *Iqac* and *Iqabc* if they need to be defined. 838 845 """ 839 846 #: Full path to the file defining the kernel, if any. … … 917 924 structure_factor = None # type: bool 918 925 #: List of C source files used to define the model. The source files 919 #: should define the *Iq* function, and possibly *Iq xy*, though a default920 #: *Iqxy = Iq(sqrt(qx**2+qy**2)* will be created if no *Iqxy* is provided.921 #: Files containing the most basic functions must appear first in the list,922 #: followed by the files that use those functions. Form factors are923 #: indicated by providing a:attr:`ER` function.926 #: should define the *Iq* function, and possibly *Iqac* or *Iqabc* if the 927 #: model defines orientation parameters. Files containing the most basic 928 #: functions must appear first in the list, followed by the files that 929 #: use those functions. Form factors are indicated by providing 930 #: an :attr:`ER` function. 924 931 source = None # type: List[str] 925 932 #: The set of tests that must pass. The format of the tests is described … … 970 977 #: include the decimal point. See :mod:`generate` for more details. 971 978 Iq = None # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 972 #: Returns *I(qx, qy, a, b, ...)*. The interface follows :attr:`Iq`. 973 Iqxy = None # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 979 #: Returns *I(qab, qc, a, b, ...)*. The interface follows :attr:`Iq`. 980 Iqac = None # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 981 #: Returns *I(qa, qb, qc, a, b, ...)*. The interface follows :attr:`Iq`. 982 Iqabc = None # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 974 983 #: Returns *I(qx, qy, a, b, ...)*. The interface follows :attr:`Iq`. 975 984 Imagnetic = None # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] … … 987 996 #: Returns a random parameter set for the model 988 997 random = None # type: Optional[Callable[[], Dict[str, float]]] 989 990 # line numbers within the python file for bits of C source, if defined 991 # NB: some compilers fail with a "#line 0" directive, so default to 1. 992 _Imagnetic_line = 1 993 _Iqxy_line = 1 994 _Iq_line = 1 995 _form_volume_line = 1 996 998 #: Line numbers for symbols defining C code 999 lineno = None # type: Dict[str, int] 997 1000 998 1001 def __init__(self): -
sasmodels/models/_cylpy.py
r59ee4db r67cc0ff 138 138 ] 139 139 140 py2c = True 141 140 142 def invalid(v): 141 143 return v.radius < 0 or v.length < 0 … … 169 171 return 1.0e-4 * s * s * orient_avg_1D(q, radius, length) 170 172 171 def Iq xy(qab, qc, sld, solvent_sld, radius, length):173 def Iqac(qab, qc, sld, solvent_sld, radius, length): 172 174 s = (sld-solvent_sld) * form_volume(radius, length) 173 175 form = fq(qab, qc, radius, length) -
sasmodels/models/_spherepy.py
ref07e95 r108e70e 88 88 Iq.vectorized = True # Iq accepts an array of q values 89 89 90 def Iqxy(qx, qy, sld, sld_solvent, radius):91 return Iq(sqrt(qx ** 2 + qy ** 2), sld, sld_solvent, radius)92 Iqxy.vectorized = True # Iqxy accepts arrays of qx, qy values93 94 90 def sesans(z, sld, sld_solvent, radius): 95 91 """ -
sasmodels/models/barbell.c
rbecded3 r108e70e 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] … … 90 90 91 91 static double 92 Iq xy(double qab, double qc,92 Iqac(double qab, double qc, 93 93 double sld, double solvent_sld, 94 94 double radius_bell, double radius, double length) -
sasmodels/models/bcc_paracrystal.c
rea60e08 r108e70e 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; … … 107 107 108 108 109 static double Iq xy(double qa, double qb, double qc,109 static double Iqabc(double qa, double qb, double qc, 110 110 double dnn, double d_factor, double radius, 111 111 double sld, double solvent_sld) -
sasmodels/models/capped_cylinder.c
rbecded3 r108e70e 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] … … 115 115 116 116 static double 117 Iq xy(double qab, double qc,117 Iqac(double qab, double qc, 118 118 double sld, double solvent_sld, double radius, 119 119 double radius_cap, double length) -
sasmodels/models/core_shell_bicelle.c
rbecded3 r108e70e 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 … … 67 67 68 68 static double 69 Iq xy(double qab, double qc,69 Iqac(double qab, double qc, 70 70 double radius, 71 71 double thick_rim, -
sasmodels/models/core_shell_bicelle_elliptical.c
r82592da r108e70e 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 … … 71 71 72 72 static double 73 Iq xy(double qa, double qb, double qc,73 Iqabc(double qa, double qb, double qc, 74 74 double r_minor, 75 75 double x_core, -
sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c
r82592da r108e70e 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 … … 79 79 80 80 static double 81 Iq xy(double qa, double qb, double qc,81 Iqabc(double qa, double qb, double qc, 82 82 double r_minor, 83 83 double x_core, … … 114 114 return 1.0e-4 * Aq*exp(-0.5*(square(qa) + square(qb) + square(qc) )*square(sigma)); 115 115 } 116 -
sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py
r110f69c r108e70e 149 149 ["sld_rim", "1e-6/Ang^2", 1, [-inf, inf], "sld", "Cylinder rim scattering length density"], 150 150 ["sld_solvent", "1e-6/Ang^2", 6, [-inf, inf], "sld", "Solvent scattering length density"], 151 ["sigma", "Ang", 0, [0, inf], "", "interfacial roughness"], 151 152 ["theta", "degrees", 90.0, [-360, 360], "orientation", "cylinder axis to beam angle"], 152 153 ["phi", "degrees", 0, [-360, 360], "orientation", "rotation about beam"], 153 154 ["psi", "degrees", 0, [-360, 360], "orientation", "rotation about cylinder axis"], 154 ["sigma", "Ang", 0, [0, inf], "", "interfacial roughness"]155 155 ] 156 156 -
sasmodels/models/core_shell_cylinder.c
rbecded3 r108e70e 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] … … 48 48 49 49 50 double Iq xy(double qab, double qc,50 double Iqac(double qab, double qc, 51 51 double core_sld, 52 52 double shell_sld, -
sasmodels/models/core_shell_ellipsoid.c
rbecded3 r108e70e 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; … … 75 75 76 76 static double 77 Iq xy(double qab, double qc,77 Iqac(double qab, double qc, 78 78 double radius_equat_core, 79 79 double x_core, -
sasmodels/models/core_shell_parallelepiped.c
rc69d6d6 r108e70e 1 // Set OVERLAPPING to 1 in order to fill in the edges of the box, with 2 // c endcaps and b overlapping a. With the proper choice of parameters, 3 // (setting rim slds to sld, core sld to solvent, rim thickness to thickness 4 // and subtracting 2*thickness from length, this should match the hollow 5 // rectangular prism.) Set it to 0 for the documented behaviour. 6 #define OVERLAPPING 0 1 7 static double 2 8 form_volume(double length_a, double length_b, double length_c, 3 9 double thick_rim_a, double thick_rim_b, double thick_rim_c) 4 10 { 5 //return length_a * length_b * length_c; 6 return length_a * length_b * length_c + 7 2.0 * thick_rim_a * length_b * length_c + 8 2.0 * thick_rim_b * length_a * length_c + 9 2.0 * thick_rim_c * length_a * length_b; 11 return 12 #if OVERLAPPING 13 // Hollow rectangular prism only includes the volume of the shell 14 // so uncomment the next line when comparing. Solid rectangular 15 // prism, or parallelepiped want filled cores, so comment when 16 // comparing. 17 //-length_a * length_b * length_c + 18 (length_a + 2.0*thick_rim_a) * 19 (length_b + 2.0*thick_rim_b) * 20 (length_c + 2.0*thick_rim_c); 21 #else 22 length_a * length_b * length_c + 23 2.0 * thick_rim_a * length_b * length_c + 24 2.0 * length_a * thick_rim_b * length_c + 25 2.0 * length_a * length_b * thick_rim_c; 26 #endif 10 27 } 11 28 … … 24 41 double thick_rim_c) 25 42 { 26 // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c _scaled43 // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c 27 44 // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 28 //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 45 // Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 46 // Code rewritten (PAK) 29 47 30 const double mu = 0.5 * q * length_b;48 const double half_q = 0.5*q; 31 49 32 //calculate volume before rescaling (in original code, but not used) 33 //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 34 //double vol = length_a * length_b * length_c + 35 // 2.0 * thick_rim_a * length_b * length_c + 36 // 2.0 * thick_rim_b * length_a * length_c + 37 // 2.0 * thick_rim_c * length_a * length_b; 50 const double tA = length_a + 2.0*thick_rim_a; 51 const double tB = length_b + 2.0*thick_rim_b; 52 const double tC = length_c + 2.0*thick_rim_c; 38 53 39 // Scale sides by B 40 const double a_scaled = length_a / length_b; 41 const double c_scaled = length_c / length_b; 42 43 double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b; 44 double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b; 45 double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present 46 47 double Vin = length_a * length_b * length_c; 48 //double Vot = (length_a * length_b * length_c + 49 // 2.0 * thick_rim_a * length_b * length_c + 50 // 2.0 * length_a * thick_rim_b * length_c + 51 // 2.0 * length_a * length_b * thick_rim_c); 52 double V1 = (2.0 * thick_rim_a * length_b * length_c); // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 53 double V2 = (2.0 * length_a * thick_rim_b * length_c); // incorrect V2(aa*bb*cc+2*aa*tb*cc) 54 double V3 = (2.0 * length_a * length_b * thick_rim_c); //not present 55 56 // Scale factors (note that drC is not used later) 57 const double drho0 = (core_sld-solvent_sld); 58 const double drhoA = (arim_sld-solvent_sld); 59 const double drhoB = (brim_sld-solvent_sld); 60 const double drhoC = (crim_sld-solvent_sld); // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot; 61 62 63 // Precompute scale factors for combining cross terms from the shape 64 const double scale23 = drhoA*V1; 65 const double scale14 = drhoB*V2; 66 const double scale24 = drhoC*V3; 67 const double scale11 = drho0*Vin; 68 const double scale12 = drho0*Vin - scale23 - scale14 - scale24; 54 // Scale factors 55 const double dr0 = (core_sld-solvent_sld); 56 const double drA = (arim_sld-solvent_sld); 57 const double drB = (brim_sld-solvent_sld); 58 const double drC = (crim_sld-solvent_sld); 69 59 70 60 // outer integral (with gauss points), integration limits = 0, 1 71 double outer_total = 0; //initialize integral 61 double outer_sum = 0; //initialize integral 62 for( int i=0; i<GAUSS_N; i++) { 63 const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 64 const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 72 65 73 for( int i=0; i<76; i++) { 74 double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 75 double mu_proj = mu * sqrt(1.0-sigma*sigma); 66 // inner integral (with gauss points), integration limits = 0, pi/2 67 const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 68 const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 69 double inner_sum = 0.0; 70 for(int j=0; j<GAUSS_N; j++) { 71 const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 ); 72 double sin_beta, cos_beta; 73 SINCOS(M_PI_2*beta, sin_beta, cos_beta); 74 const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 75 const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); 76 const double siAt = tA * sas_sinx_x(tA * mu * sin_beta); 77 const double siBt = tB * sas_sinx_x(tB * mu * cos_beta); 76 78 77 // inner integral (with gauss points), integration limits = 0, 1 78 double inner_total = 0.0;79 double inner_total_crim = 0.0;80 for(int j=0; j<76; j++) {81 const double uu = 0.5 * ( Gauss76Z[j] + 1.0);82 double sin_uu, cos_uu; 83 SINCOS(M_PI_2*uu, sin_uu, cos_uu);84 const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled);85 const double si2 = sas_sinx_x(mu_proj * cos_uu );86 const double si3 = sas_sinx_x(mu_proj * sin_uu * ta);87 const double si4 = sas_sinx_x(mu_proj * cos_uu * tb); 79 #if OVERLAPPING 80 const double f = dr0*siA*siB*siC 81 + drA*(siAt-siA)*siB*siC 82 + drB*siAt*(siBt-siB)*siC 83 + drC*siAt*siBt*(siCt-siC); 84 #else 85 const double f = dr0*siA*siB*siC 86 + drA*(siAt-siA)*siB*siC 87 + drB*siA*(siBt-siB)*siC 88 + drC*siA*siB*(siCt-siC); 89 #endif 88 90 89 // Expression in libCylinder.c (neither drC nor Vot are used) 90 const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; 91 const double form_crim = scale11*si1*si2; 92 93 // correct FF : sum of square of phase factors 94 inner_total += Gauss76Wt[j] * form * form; 95 inner_total_crim += Gauss76Wt[j] * form_crim * form_crim; 91 inner_sum += GAUSS_W[j] * f * f; 96 92 } 97 inner_total *= 0.5; 98 inner_total_crim *= 0.5; 93 inner_sum *= 0.5; 99 94 // now sum up the outer integral 100 const double si = sas_sinx_x(mu * c_scaled * sigma); 101 const double si_crim = sas_sinx_x(mu * tc * sigma); 102 outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim); 95 outer_sum += GAUSS_W[i] * inner_sum; 103 96 } 104 outer_ total*= 0.5;97 outer_sum *= 0.5; 105 98 106 99 //convert from [1e-12 A-1] to [cm-1] 107 return 1.0e-4 * outer_ total;100 return 1.0e-4 * outer_sum; 108 101 } 109 102 110 103 static double 111 Iq xy(double qa, double qb, double qc,104 Iqabc(double qa, double qb, double qc, 112 105 double core_sld, 113 106 double arim_sld, … … 128 121 const double drC = crim_sld-solvent_sld; 129 122 130 double Vin = length_a * length_b * length_c;131 double V1 = 2.0 * thick_rim_a * length_b * length_c; // incorrect V1(aa*bb*cc+2*ta*bb*cc)132 double V2 = 2.0 * length_a * thick_rim_b * length_c; // incorrect V2(aa*bb*cc+2*aa*tb*cc)133 double V3 = 2.0 * length_a * length_b * thick_rim_c;134 // As for the 1D case, Vot is not used135 //double Vot = (length_a * length_b * length_c +136 // 2.0 * thick_rim_a * length_b * length_c +137 // 2.0 * length_a * thick_rim_b * length_c +138 // 2.0 * length_a * length_b * thick_rim_c);139 140 123 // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 141 124 // the scaling by B. 142 double ta = length_a + 2.0*thick_rim_a; 143 double tb = length_b + 2.0*thick_rim_b; 144 double tc = length_c + 2.0*thick_rim_c; 145 //handle arg=0 separately, as sin(t)/t -> 1 as t->0 146 double siA = sas_sinx_x(0.5*length_a*qa); 147 double siB = sas_sinx_x(0.5*length_b*qb); 148 double siC = sas_sinx_x(0.5*length_c*qc); 149 double siAt = sas_sinx_x(0.5*ta*qa); 150 double siBt = sas_sinx_x(0.5*tb*qb); 151 double siCt = sas_sinx_x(0.5*tc*qc); 125 const double tA = length_a + 2.0*thick_rim_a; 126 const double tB = length_b + 2.0*thick_rim_b; 127 const double tC = length_c + 2.0*thick_rim_c; 128 const double siA = length_a*sas_sinx_x(0.5*length_a*qa); 129 const double siB = length_b*sas_sinx_x(0.5*length_b*qb); 130 const double siC = length_c*sas_sinx_x(0.5*length_c*qc); 131 const double siAt = tA*sas_sinx_x(0.5*tA*qa); 132 const double siBt = tB*sas_sinx_x(0.5*tB*qb); 133 const double siCt = tC*sas_sinx_x(0.5*tC*qc); 152 134 153 154 // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 155 // in the 1D code, but should be checked! 156 double f = ( dr0*siA*siB*siC*Vin 157 + drA*(siAt-siA)*siB*siC*V1 158 + drB*siA*(siBt-siB)*siC*V2 159 + drC*siA*siB*(siCt-siC)*V3); 135 #if OVERLAPPING 136 const double f = dr0*siA*siB*siC 137 + drA*(siAt-siA)*siB*siC 138 + drB*siAt*(siBt-siB)*siC 139 + drC*siAt*siBt*(siCt-siC); 140 #else 141 const double f = dr0*siA*siB*siC 142 + drA*(siAt-siA)*siB*siC 143 + drB*siA*(siBt-siB)*siC 144 + drC*siA*siB*(siCt-siC); 145 #endif 160 146 161 147 return 1.0e-4 * f * f; -
sasmodels/models/core_shell_parallelepiped.py
r2d81cfe r10ee838 5 5 Calculates the form factor for a rectangular solid with a core-shell structure. 6 6 The thickness and the scattering length density of the shell or 7 "rim" can be different on each (pair) of faces. However at this time the 1D 8 calculation does **NOT** actually calculate a c face rim despite the presence 9 of the parameter. Some other aspects of the 1D calculation may be wrong. 10 11 .. note:: 12 This model was originally ported from NIST IGOR macros. However, it is not 13 yet fully understood by the SasView developers and is currently under review. 7 "rim" can be different on each (pair) of faces. 14 8 15 9 The form factor is normalized by the particle volume $V$ such that … … 21 15 where $\langle \ldots \rangle$ is an average over all possible orientations 22 16 of the rectangular solid. 23 24 17 25 18 The function calculated is the form factor of the rectangular solid below. … … 41 34 V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 42 35 43 **meaning that there are "gaps" at the corners of the solid.** Again note that 44 $t_C = 0$ currently. 36 **meaning that there are "gaps" at the corners of the solid.** 45 37 46 38 The intensity calculated follows the :ref:`parallelepiped` model, with the 47 39 core-shell intensity being calculated as the square of the sum of the 48 amplitudes of the core and shell, in the same manner as a core-shell model. 49 50 .. math:: 51 52 F_{a}(Q,\alpha,\beta)= 53 \left[\frac{\sin(\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha \sin\beta) 54 }{\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha\sin\beta} 55 - \frac{\sin(\tfrac{1}{2}QL_A\sin\alpha \sin\beta) 56 }{\tfrac{1}{2}QL_A\sin\alpha \sin\beta} \right] 57 \left[\frac{\sin(\tfrac{1}{2}QL_B\sin\alpha \sin\beta) 58 }{\tfrac{1}{2}QL_B\sin\alpha \sin\beta} \right] 59 \left[\frac{\sin(\tfrac{1}{2}QL_C\sin\alpha \sin\beta) 60 }{\tfrac{1}{2}QL_C\sin\alpha \sin\beta} \right] 61 62 .. note:: 63 64 Why does t_B not appear in the above equation? 65 For the calculation of the form factor to be valid, the sides of the solid 66 MUST (perhaps not any more?) be chosen such that** $A < B < C$. 67 If this inequality is not satisfied, the model will not report an error, 68 but the calculation will not be correct and thus the result wrong. 40 amplitudes of the core and the slabs on the edges. 41 42 the scattering amplitude is computed for a particular orientation of the 43 core-shell parallelepiped with respect to the scattering vector and then 44 averaged over all possible orientations, where $\alpha$ is the angle between 45 the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 46 the angle between projection of the particle in the $xy$ detector plane 47 and the $y$ axis. 48 49 .. math:: 50 51 F(Q) 52 &= (\rho_\text{core}-\rho_\text{solvent}) 53 S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 54 &+ (\rho_\text{A}-\rho_\text{solvent}) 55 \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\ 56 &+ (\rho_\text{B}-\rho_\text{solvent}) 57 S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ 58 &+ (\rho_\text{C}-\rho_\text{solvent}) 59 S(Q_A, A) S(Q_B, B) \left[S(Q_C, C+2t_C) - S(Q_C, C)\right] 60 61 with 62 63 .. math:: 64 65 S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 66 67 and 68 69 .. math:: 70 71 Q_A &= \sin\alpha \sin\beta \\ 72 Q_B &= \sin\alpha \cos\beta \\ 73 Q_C &= \cos\alpha 74 75 76 where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 77 are the scattering length of the parallelepiped core, and the rectangular 78 slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 79 is the scattering length of the solvent. 69 80 70 81 FITTING NOTES 82 ~~~~~~~~~~~~~ 83 71 84 If the scale is set equal to the particle volume fraction, $\phi$, the returned 72 value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 73 However, **no interparticle interference effects are included in this 74 calculation.** 85 value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However, 86 **no interparticle interference effects are included in this calculation.** 75 87 76 88 There are many parameters in this model. Hold as many fixed as possible with 77 89 known values, or you will certainly end up at a solution that is unphysical. 78 90 79 Constraints must be applied during fitting to ensure that the inequality80 $A < B < C$ is not violated. The calculation will not report an error,81 but the results will not be correct.82 83 91 The returned value is in units of |cm^-1|, on absolute scale. 84 92 85 93 NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 86 94 based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 87 and length $(C+2t_C)$ values, after appropriately 88 sorting the three dimensions to give an oblate or prolate particle, to give an 89 effective radius,for $S(Q)$ when $P(Q) * S(Q)$ is applied.95 and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 96 to give an oblate or prolate particle, to give an effective radius, 97 for $S(Q)$ when $P(Q) * S(Q)$ is applied. 90 98 91 99 For 2d data the orientation of the particle is required, described using 92 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details93 of the calculation and angular dispersions see :ref:`orientation`.100 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 101 details of the calculation and angular dispersions see :ref:`orientation`. 94 102 The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 95 103 $\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 104 105 For 2d, constraints must be applied during fitting to ensure that the 106 inequality $A < B < C$ is not violated, and hence the correct definition 107 of angles is preserved. The calculation will not report an error, 108 but the results may be not correct. 96 109 97 110 .. figure:: img/parallelepiped_angle_definition.png … … 114 127 Equations (1), (13-14). (in German) 115 128 .. [#] D Singh (2009). *Small angle scattering studies of self assembly in 116 lipid mixtures*, John 's Hopkins University Thesis (2009) 223-225. `Available129 lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available 117 130 from Proquest <http://search.proquest.com/docview/304915826?accountid 118 131 =26379>`_ … … 175 188 Return equivalent radius (ER) 176 189 """ 177 178 # surface average radius (rough approximation) 179 surf_rad = sqrt((length_a + 2.0*thick_rim_a) * (length_b + 2.0*thick_rim_b) / pi) 180 181 height = length_c + 2.0*thick_rim_c 182 183 ddd = 0.75 * surf_rad * (2 * surf_rad * height + (height + surf_rad) * (height + pi * surf_rad)) 184 return 0.5 * (ddd) ** (1. / 3.) 190 from .parallelepiped import ER as ER_p 191 192 a = length_a + 2*thick_rim_a 193 b = length_b + 2*thick_rim_b 194 c = length_c + 2*thick_rim_c 195 return ER_p(a, b, c) 185 196 186 197 # VR defaults to 1.0 … … 216 227 psi_pd=10, psi_pd_n=1) 217 228 218 # rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models 229 # rkh 7/4/17 add random unit test for 2d, note make all params different, 230 # 2d values not tested against other codes or models 219 231 if 0: # pak: model rewrite; need to update tests 220 232 qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) -
sasmodels/models/cylinder.c
rbecded3 r108e70e 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] … … 45 45 46 46 static double 47 Iq xy(double qab, double qc,47 Iqac(double qab, double qc, 48 48 double sld, 49 49 double solvent_sld, -
sasmodels/models/ellipsoid.c
rbecded3 r108e70e 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] … … 39 39 40 40 static double 41 Iq xy(double qab, double qc,41 Iqac(double qab, double qc, 42 42 double sld, 43 43 double sld_solvent, -
sasmodels/models/elliptical_cylinder.c
r82592da r108e70e 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); … … 55 54 56 55 static double 57 Iq xy(double qa, double qb, double qc,56 Iqabc(double qa, double qb, double qc, 58 57 double radius_minor, double r_ratio, double length, 59 58 double sld, double solvent_sld) -
sasmodels/models/elliptical_cylinder.py
r2d81cfe r2d81cfe 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 r108e70e 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; … … 80 80 81 81 82 static double Iq xy(double qa, double qb, double qc,82 static double Iqabc(double qa, double qb, double qc, 83 83 double dnn, double d_factor, double radius, 84 84 double sld, double solvent_sld) -
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 r108e70e 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 … … 52 52 53 53 static double 54 Iq xy(double qab, double qc,54 Iqac(double qab, double qc, 55 55 double radius, double thickness, double length, 56 56 double sld, double solvent_sld) -
sasmodels/models/hollow_rectangular_prism.c
r1e7b0db0 r108e70e 1 1 double form_volume(double length_a, double b2a_ratio, double c2a_ratio, double thickness); 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, double thickness); 4 4 … … 37 37 const double v2a = 0.0; 38 38 const double v2b = M_PI_2; //phi integration limits 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); … … 84 84 return 1.0e-4 * delrho * delrho * form; 85 85 } 86 87 double Iqabc(double qa, double qb, double qc, 88 double sld, 89 double solvent_sld, 90 double length_a, 91 double b2a_ratio, 92 double c2a_ratio, 93 double thickness) 94 { 95 const double length_b = length_a * b2a_ratio; 96 const double length_c = length_a * c2a_ratio; 97 const double a_half = 0.5 * length_a; 98 const double b_half = 0.5 * length_b; 99 const double c_half = 0.5 * length_c; 100 const double vol_total = length_a * length_b * length_c; 101 const double vol_core = 8.0 * (a_half-thickness) * (b_half-thickness) * (c_half-thickness); 102 103 // Amplitude AP from eqn. (13) 104 105 const double termA1 = sas_sinx_x(qa * a_half); 106 const double termA2 = sas_sinx_x(qa * (a_half-thickness)); 107 108 const double termB1 = sas_sinx_x(qb * b_half); 109 const double termB2 = sas_sinx_x(qb * (b_half-thickness)); 110 111 const double termC1 = sas_sinx_x(qc * c_half); 112 const double termC2 = sas_sinx_x(qc * (c_half-thickness)); 113 114 const double AP1 = vol_total * termA1 * termB1 * termC1; 115 const double AP2 = vol_core * termA2 * termB2 * termC2; 116 117 // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 118 const double delrho = sld - solvent_sld; 119 120 // Convert from [1e-12 A-1] to [cm-1] 121 return 1.0e-4 * square(delrho * (AP1-AP2)); 122 } -
sasmodels/models/hollow_rectangular_prism.py
r2d81cfe r2d81cfe 5 5 This model provides the form factor, $P(q)$, for a hollow rectangular 6 6 parallelepiped with a wall of thickness $\Delta$. 7 It computes only the 1D scattering, not the 2D. 7 8 8 9 9 Definition … … 66 66 (which is unitless). 67 67 68 **The 2D scattering intensity is not computed by this model.** 68 For 2d data the orientation of the particle is required, described using 69 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 70 of the calculation and angular dispersions see :ref:`orientation` . 71 The angle $\Psi$ is the rotational angle around the long *C* axis. For example, 72 $\Psi = 0$ when the *B* axis is parallel to the *x*-axis of the detector. 73 74 For 2d, constraints must be applied during fitting to ensure that the inequality 75 $A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error, 76 but the results may be not correct. 77 78 .. figure:: img/parallelepiped_angle_definition.png 79 80 Definition of the angles for oriented hollow rectangular prism. 81 Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then 82 rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the prism. 83 The neutron or X-ray beam is along the $z$ axis. 84 85 .. figure:: img/parallelepiped_angle_projection.png 86 87 Examples of the angles for oriented hollow rectangular prisms against the 88 detector plane. 69 89 70 90 … … 113 133 ["thickness", "Ang", 1, [0, inf], "volume", 114 134 "Thickness of parallelepiped"], 135 ["theta", "degrees", 0, [-360, 360], "orientation", 136 "c axis to beam angle"], 137 ["phi", "degrees", 0, [-360, 360], "orientation", 138 "rotation about beam"], 139 ["psi", "degrees", 0, [-360, 360], "orientation", 140 "rotation about c axis"], 115 141 ] 116 142 -
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/line.py
r2d81cfe r108e70e 57 57 Iq.vectorized = True # Iq accepts an array of q values 58 58 59 59 60 def Iqxy(qx, qy, *args): 60 61 """ … … 69 70 70 71 Iqxy.vectorized = True # Iqxy accepts an array of qx qy values 72 73 # uncomment the following to test Iqxy in C models 74 #del Iq, Iqxy 75 #c_code = """ 76 #static double Iq(double q, double b, double m) { return m*q+b; } 77 #static double Iqxy(double qx, double qy, double b, double m) 78 #{ return (m*qx+b)*(m*qy+b); } 79 #""" 71 80 72 81 def random(): -
sasmodels/models/parallelepiped.c
r9b7b23f r108e70e 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; … … 53 53 54 54 static double 55 Iq xy(double qa, double qb, double qc,55 Iqabc(double qa, double qb, double qc, 56 56 double sld, 57 57 double solvent_sld, -
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
r1e7b0db0 r108e70e 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 … … 26 26 const double v2a = 0.0; 27 27 const double v2b = M_PI_2; //phi integration limits 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 53 53 double answer = 0.5*(v1b-v1a)*outer_sum; 54 54 55 // Normalize by Pi (Eqn. 16). 56 // The term (ABC)^2 does not appear because it was introduced before on 55 // Normalize by Pi (Eqn. 16). 56 // The term (ABC)^2 does not appear because it was introduced before on 57 57 // the definitions of termA, termB, termC. 58 // The factor 2 appears because the theta integral has been defined between 58 // The factor 2 appears because the theta integral has been defined between 59 59 // 0 and pi/2, instead of 0 to pi. 60 60 answer /= M_PI_2; //Form factor P(q) … … 64 64 answer *= square((sld-solvent_sld)*volume); 65 65 66 // Convert from [1e-12 A-1] to [cm-1] 66 // Convert from [1e-12 A-1] to [cm-1] 67 67 answer *= 1.0e-4; 68 68 69 69 return answer; 70 70 } 71 72 73 double Iqabc(double qa, double qb, double qc, 74 double sld, 75 double solvent_sld, 76 double length_a, 77 double b2a_ratio, 78 double c2a_ratio) 79 { 80 const double length_b = length_a * b2a_ratio; 81 const double length_c = length_a * c2a_ratio; 82 const double a_half = 0.5 * length_a; 83 const double b_half = 0.5 * length_b; 84 const double c_half = 0.5 * length_c; 85 const double volume = length_a * length_b * length_c; 86 87 // Amplitude AP from eqn. (13) 88 89 const double termA = sas_sinx_x(qa * a_half); 90 const double termB = sas_sinx_x(qb * b_half); 91 const double termC = sas_sinx_x(qc * c_half); 92 93 const double AP = termA * termB * termC; 94 95 // Multiply by contrast^2. Factor corresponding to volume^2 cancels with previous normalization. 96 const double delrho = sld - solvent_sld; 97 98 // Convert from [1e-12 A-1] to [cm-1] 99 return 1.0e-4 * square(volume * delrho * AP); 100 } -
sasmodels/models/rectangular_prism.py
r2d81cfe r2d81cfe 12 12 the prism (e.g. setting $b/a = 1$ and $c/a = 1$ and applying polydispersity 13 13 to *a* will generate a distribution of cubes of different sizes). 14 Note also that, contrary to :ref:`parallelepiped`, it does not compute15 the 2D scattering.16 14 17 15 … … 26 24 that reference), with $\theta$ corresponding to $\alpha$ in that paper, 27 25 and not to the usual convention used for example in the 28 :ref:`parallelepiped` model. As the present model does not compute 29 the 2D scattering, this has no further consequences. 26 :ref:`parallelepiped` model. 30 27 31 28 In this model the scattering from a massive parallelepiped with an … … 65 62 units) *scale* represents the volume fraction (which is unitless). 66 63 67 **The 2D scattering intensity is not computed by this model.** 64 For 2d data the orientation of the particle is required, described using 65 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 66 of the calculation and angular dispersions see :ref:`orientation` . 67 The angle $\Psi$ is the rotational angle around the long *C* axis. For example, 68 $\Psi = 0$ when the *B* axis is parallel to the *x*-axis of the detector. 69 70 For 2d, constraints must be applied during fitting to ensure that the inequality 71 $A < B < C$ is not violated, and hence the correct definition of angles is preserved. The calculation will not report an error, 72 but the results may be not correct. 73 74 .. figure:: img/parallelepiped_angle_definition.png 75 76 Definition of the angles for oriented core-shell parallelepipeds. 77 Note that rotation $\theta$, initially in the $xz$ plane, is carried out first, then 78 rotation $\phi$ about the $z$ axis, finally rotation $\Psi$ is now around the axis of the cylinder. 79 The neutron or X-ray beam is along the $z$ axis. 80 81 .. figure:: img/parallelepiped_angle_projection.png 82 83 Examples of the angles for oriented rectangular prisms against the 84 detector plane. 85 68 86 69 87 … … 108 126 ["c2a_ratio", "", 1, [0, inf], "volume", 109 127 "Ratio sides c/a"], 128 ["theta", "degrees", 0, [-360, 360], "orientation", 129 "c axis to beam angle"], 130 ["phi", "degrees", 0, [-360, 360], "orientation", 131 "rotation about beam"], 132 ["psi", "degrees", 0, [-360, 360], "orientation", 133 "rotation about c axis"], 110 134 ] 111 135 -
sasmodels/models/sc_paracrystal.c
rf728001 r108e70e 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; … … 82 82 83 83 static double 84 Iq xy(double qa, double qb, double qc,84 Iqabc(double qa, double qb, double qc, 85 85 double dnn, double d_factor, double radius, 86 86 double sld, double solvent_sld) -
sasmodels/models/stacked_disks.c
rbecded3 r108e70e 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 … … 142 142 143 143 static double 144 Iq xy(double qab, double qc,144 Iqac(double qab, double qc, 145 145 double thick_core, 146 146 double thick_layer, -
sasmodels/models/triaxial_ellipsoid.c
rbecded3 r108e70e 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 … … 46 46 47 47 static double 48 Iq xy(double qa, double qb, double qc,48 Iqabc(double qa, double qb, double qc, 49 49 double sld, 50 50 double sld_solvent, -
sasmodels/py2c.py
rddfdb16 r7b1dcf9 81 81 BOOLOP_SYMBOLS = {} 82 82 BOOLOP_SYMBOLS[ast.And] = '&&' 83 BOOLOP_SYMBOLS[ast.Or] 83 BOOLOP_SYMBOLS[ast.Or] = '||' 84 84 85 85 CMPOP_SYMBOLS = {} 86 CMPOP_SYMBOLS[ast.Eq] 86 CMPOP_SYMBOLS[ast.Eq] = '==' 87 87 CMPOP_SYMBOLS[ast.NotEq] = '!=' 88 88 CMPOP_SYMBOLS[ast.Lt] = '<' … … 145 145 self.new_lines = 0 146 146 self.c_proc = [] 147 # for C147 # for C 148 148 self.signature_line = 0 149 149 self.arguments = [] … … 181 181 def add_c_line(self, x): 182 182 string = '' 183 for iin range(self.indentation):183 for _ in range(self.indentation): 184 184 string += (" ") 185 185 string += str(x) … … 188 188 189 189 def add_current_line(self): 190 if (len(self.current_statement) > 0):190 if self.current_statement: 191 191 self.add_c_line(self.current_statement) 192 192 self.current_statement = '' 193 193 194 194 def AddUniqueVar(self, new_var): 195 if ((new_var not in self.C_Vars)):195 if new_var not in self.C_Vars: 196 196 self.C_Vars.append(str(new_var)) 197 197 … … 210 210 self.write_c('# line: %s' % node.lineno) 211 211 self.new_lines = 1 212 if (len(self.current_statement)):212 if self.current_statement: 213 213 self.Statements.append(self.current_statement) 214 214 self.current_statement = '' 215 215 216 216 def body(self, statements): 217 if (len(self.current_statement)):217 if self.current_statement: 218 218 self.add_current_line() 219 219 self.new_line = True 220 220 self.indentation += 1 221 221 for stmt in statements: 222 target_name = '' 223 if(hasattr(stmt, 'targets')): 224 if(hasattr(stmt.targets[0], 'id')): 225 target_name = stmt.targets[0].id # target name needed for debug only 222 #if hasattr(stmt, 'targets') and hasattr(stmt.targets[0], 'id'): 223 # target_name = stmt.targets[0].id # target name needed for debug only 226 224 self.visit(stmt) 227 225 self.add_current_line() # just for breaking point. to be deleted. … … 242 240 else: 243 241 want_comma.append(True) 244 # for C 242 243 # for C 245 244 for arg in node.args: 246 self.arguments.append(arg.arg) 245 # CRUFT: 2.7 uses arg.id, 3.x uses arg.arg 246 try: 247 arg_name = arg.arg 248 except AttributeError: 249 arg_name = arg.id 250 self.arguments.append(arg_name) 247 251 248 252 padding = [None] *(len(node.args) - len(node.defaults)) 249 253 for arg, default in zip(node.args, padding + node.defaults): 250 254 if default is not None: 251 self.warnings.append("Default Parameter unknown to C") 252 w_str = "Default Parameters are unknown to C: '" + arg.arg + \ 253 " = " + str(default.n) + "'" 255 # CRUFT: 2.7 uses arg.id, 3.x uses arg.arg 256 try: 257 arg_name = arg.arg 258 except AttributeError: 259 arg_name = arg.id 260 w_str = ("Default Parameters are unknown to C: '%s = %s" 261 % arg_name, str(default.n)) 254 262 self.warnings.append(w_str) 255 263 … … 271 279 272 280 def define_C_Vars(self, target): 273 if (hasattr(target, 'id')):274 # a variable is considered an array if it apears in the agrument list275 # and being assigned to. For example, the variable p in the following276 # sniplet is a pointer, while q is not277 # def somefunc(p, q):278 # p = q + 1279 # return280 #281 if (target.id not in self.C_Vars):282 if (target.id in self.arguments):281 if hasattr(target, 'id'): 282 # a variable is considered an array if it apears in the agrument list 283 # and being assigned to. For example, the variable p in the following 284 # sniplet is a pointer, while q is not 285 # def somefunc(p, q): 286 # p = q + 1 287 # return 288 # 289 if target.id not in self.C_Vars: 290 if target.id in self.arguments: 283 291 idx = self.arguments.index(target.id) 284 292 new_target = self.arguments[idx] + "[0]" 285 if (new_target not in self.C_Pointers):293 if new_target not in self.C_Pointers: 286 294 target.id = new_target 287 295 self.C_Pointers.append(self.arguments[idx]) … … 291 299 def add_semi_colon(self): 292 300 semi_pos = self.current_statement.find(';') 293 if (semi_pos > 0.0):294 self.current_statement = self.current_statement.replace(';', '')301 if semi_pos > 0.0: 302 self.current_statement = self.current_statement.replace(';', '') 295 303 self.write_c(';') 296 304 … … 302 310 self.define_C_Vars(target) 303 311 self.visit(target) 304 if (len(self.Tuples) > 0):312 if self.Tuples: 305 313 tplTargets = list(self.Tuples) 306 self.Tuples.clear()314 del self.Tuples[:] 307 315 self.write_c(' = ') 308 316 self.is_sequence = False … … 317 325 self.add_semi_colon() 318 326 self.add_current_line() 319 if ((self.is_sequence) and (not self.visited_args)):327 if self.is_sequence and not self.visited_args: 320 328 for target in node.targets: 321 if (hasattr(target, 'id')):322 if ((target.id in self.C_Vars) and(target.id not in self.C_DclPointers)):323 if (target.id not in self.C_DclPointers):329 if hasattr(target, 'id'): 330 if target.id in self.C_Vars and target.id not in self.C_DclPointers: 331 if target.id not in self.C_DclPointers: 324 332 self.C_DclPointers.append(target.id) 325 if (target.id in self.C_Vars):333 if target.id in self.C_Vars: 326 334 self.C_Vars.remove(target.id) 327 335 self.current_statement = '' 328 336 329 337 def visit_AugAssign(self, node): 330 if (node.target.id not in self.C_Vars):331 if (node.target.id not in self.arguments):338 if node.target.id not in self.C_Vars: 339 if node.target.id not in self.arguments: 332 340 self.C_Vars.append(node.target.id) 333 341 self.visit(node.target) … … 355 363 self.generic_visit(node) 356 364 357 def listToDeclare(self, Vars): 358 s = '' 359 if(len(Vars) > 0): 360 s = ",".join(Vars) 361 return(s) 365 def listToDeclare(self, vars): 366 return ", ".join(vars) 362 367 363 368 def write_C_Pointers(self, start_var): 364 if (len(self.C_DclPointers) > 0):365 var s = ""369 if self.C_DclPointers: 370 var_list = [] 366 371 for c_ptr in self.C_DclPointers: 367 372 if(len(vars) > 0): 368 373 vars += ", " 369 if(c_ptr not in self.arguments): 370 vars += "*" + c_ptr 371 if(c_ptr in self.C_Vars): 372 if(c_ptr in self.C_Vars): 373 self.C_Vars.remove(c_ptr) 374 if(len(vars) > 0): 375 c_dcl = " double " + vars + ";" 376 self.c_proc.insert(start_var, c_dcl + "\n") 374 if c_ptr not in self.arguments: 375 var_list.append("*" + c_ptr) 376 if c_ptr in self.C_Vars: 377 self.C_Vars.remove(c_ptr) 378 if var_list: 379 c_dcl = " double " + ", ".join(var_list) + ";\n" 380 self.c_proc.insert(start_var, c_dcl) 377 381 start_var += 1 378 382 return start_var … … 381 385 fLine = False 382 386 start_var = self.write_C_Pointers(start_var) 383 if (len(self.C_IntVars) > 0):387 if self.C_IntVars: 384 388 for var in self.C_IntVars: 385 if (var in self.C_Vars):389 if var in self.C_Vars: 386 390 self.C_Vars.remove(var) 387 391 s = self.listToDeclare(self.C_IntVars) … … 390 394 start_var += 1 391 395 392 if (len(self.C_Vars) > 0):396 if self.C_Vars: 393 397 s = self.listToDeclare(self.C_Vars) 394 398 self.c_proc.insert(start_var, " double " + s + ";\n") 395 399 fLine = True 396 400 start_var += 1 397 if(len(self.C_Vectors) > 0): 401 402 if self.C_Vectors: 398 403 s = self.listToDeclare(self.C_Vectors) 399 404 for n in range(len(self.C_Vectors)): … … 402 407 self.c_proc.insert(start_var, c_dcl + "\n") 403 408 start_var += 1 404 self.C_Vars.clear() 405 self.C_IntVars.clear() 406 self.C_Vectors.clear() 407 self.C_Pointers.clear() 409 410 del self.C_Vars[:] 411 del self.C_IntVars[:] 412 del self.C_Vectors[:] 413 del self.C_Pointers[:] 408 414 self.C_DclPointers 409 if (fLine == True):415 if fLine: 410 416 self.c_proc.insert(start_var, "\n") 411 return412 s = ''413 for n in range(len(self.C_Vars)):414 s += str(self.C_Vars[n])415 if n < len(self.C_Vars) - 1:416 s += ", "417 if(len(s) > 0):418 self.c_proc.insert(start_var, " double " + s + ";\n")419 self.c_proc.insert(start_var + 1, "\n")420 421 def writeInclude(self):422 if(self.MathIncludeed == False):423 self.add_c_line("#include <math.h>\n")424 self.add_c_line("static double pi = 3.14159265359;\n")425 self.MathIncludeed = True426 427 def ListToString(self, strings):428 s = ''429 for n in range(len(strings)):430 s += strings[n]431 if(n < (len(strings) - 1)):432 s += ", "433 return(s)434 435 def getMethodSignature(self):436 args_str = ''437 for n in range(len(self.arguments)):438 args_str += "double " + self.arguments[n]439 if(n < (len(self.arguments) - 1)):440 args_str += ", "441 return(args_str)442 417 443 418 def InsertSignature(self): 444 arg s_str = ''445 for n in range(len(self.arguments)):446 args_str += "double " + self.arguments[n]447 if (self.arguments[n] in self.C_Pointers):448 args_str+= "[]"449 if(n < (len(self.arguments) - 1)):450 args_str += ", "419 arg_decls = [] 420 for arg in self.arguments: 421 decl = "double " + arg 422 if arg in self.C_Pointers: 423 decl += "[]" 424 arg_decls.append(decl) 425 args_str = ", ".join(arg_decls) 451 426 self.strMethodSignature = 'double ' + self.name + '(' + args_str + ")" 452 if (self.signature_line >= 0):427 if self.signature_line >= 0: 453 428 self.c_proc.insert(self.signature_line, self.strMethodSignature) 454 429 … … 459 434 self.arguments = [] 460 435 self.name = node.name 461 print("Parsing '" + self.name + "'") 462 args_str = "" 436 #if self.name not in self.required_functions[0]: 437 # return 438 #print("Parsing '" + self.name + "'") 463 439 464 440 self.visit(node.args) 465 self.getMethodSignature()441 # for C 466 442 self.signature_line = len(self.c_proc) 443 #self.add_c_line(self.strMethodSignature) 467 444 self.add_c_line("\n{") 468 445 start_vars = len(self.c_proc) + 1 … … 489 466 paren_or_comma() 490 467 self.visit(base) 491 # XXX: the if here is used to keep this module compatible 492 # with python 2.6. 468 # CRUFT: python 2.6 does not have "keywords" attribute 493 469 if hasattr(node, 'keywords'): 494 470 for keyword in node.keywords: … … 517 493 if len(else_) == 0: 518 494 break 519 #elif hasattr(else_, 'orelse'):495 #elif hasattr(else_, 'orelse'): 520 496 elif len(else_) == 1 and isinstance(else_[0], ast.If): 521 497 node = else_[0] 522 #self.newline()498 #self.newline() 523 499 self.write_c('else if ') 524 500 self.visit(node.test) … … 527 503 self.add_current_line() 528 504 self.add_c_line('}') 529 #break505 #break 530 506 else: 531 507 self.newline() … … 537 513 def getNodeLineNo(self, node): 538 514 line_number = -1 539 if (hasattr(node,'value')):515 if hasattr(node, 'value'): 540 516 line_number = node.value.lineno 541 517 elif hasattr(node, 'iter'): 542 518 if hasattr(node.iter, 'lineno'): 543 519 line_number = node.iter.lineno 544 return (line_number)520 return line_number 545 521 546 522 def GetNodeAsString(self, node): 547 523 res = '' 548 if (hasattr(node, 'n')):524 if hasattr(node, 'n'): 549 525 res = str(node.n) 550 elif (hasattr(node, 'id')):526 elif hasattr(node, 'id'): 551 527 res = node.id 552 return (res)528 return res 553 529 554 530 def GetForRange(self, node): … … 564 540 self.current_statement = '' 565 541 self.current_statement = temp_statement 566 if (len(for_args) == 1):542 if len(for_args) == 1: 567 543 stop = for_args[0] 568 elif (len(for_args) == 2):544 elif len(for_args) == 2: 569 545 start = for_args[0] 570 546 stop = for_args[1] 571 elif (len(for_args) == 3):547 elif len(for_args) == 3: 572 548 start = for_args[0] 573 549 stop = for_args[1] … … 575 551 else: 576 552 raise("Ilegal for loop parameters") 577 return (start, stop, step)553 return start, stop, step 578 554 579 555 def visit_For(self, node): 580 # node: for iterator is stored in node.target.581 # Iterator name is in node.target.id.556 # node: for iterator is stored in node.target. 557 # Iterator name is in node.target.id. 582 558 self.add_current_line() 583 559 fForDone = False 584 560 self.current_statement = '' 585 if (hasattr(node.iter, 'func')):586 if (hasattr(node.iter.func, 'id')):587 if (node.iter.func.id == 'range'):561 if hasattr(node.iter, 'func'): 562 if hasattr(node.iter.func, 'id'): 563 if node.iter.func.id == 'range': 588 564 self.visit(node.target) 589 565 iterator = self.current_statement 590 566 self.current_statement = '' 591 if (iterator not in self.C_IntVars):567 if iterator not in self.C_IntVars: 592 568 self.C_IntVars.append(iterator) 593 569 start, stop, step = self.GetForRange(node) 594 self.write_c("for(" + iterator + "=" + str(start) + \595 " ; " + iterator + " < " + str(stop) + \596 570 self.write_c("for(" + iterator + "=" + str(start) + 571 " ; " + iterator + " < " + str(stop) + 572 " ; " + iterator + " += " + str(step) + ") {") 597 573 self.body_or_else(node) 598 574 self.write_c("}") 599 575 fForDone = True 600 if (fForDone == False):576 if not fForDone: 601 577 line_number = self.getNodeLineNo(node) 602 578 self.current_statement = '' … … 632 608 633 609 def visit_Print(self, node): 634 # XXX: python 2.6 only610 # CRUFT: python 2.6 only 635 611 self.newline(node) 636 612 self.write_c('print ') … … 700 676 701 677 def visit_Raise(self, node): 702 # XXX: Python 2.6 / 3.0 compatibility678 # CRUFT: Python 2.6 / 3.0 compatibility 703 679 self.newline(node) 704 680 self.write_python('raise') … … 734 710 else: 735 711 want_comma.append(True) 736 if (hasattr(node.func, 'id')):737 if (node.func.id not in self.C_Functions):712 if hasattr(node.func, 'id'): 713 if node.func.id not in self.C_Functions: 738 714 self.C_Functions.append(node.func.id) 739 if (node.func.id == 'abs'):715 if node.func.id == 'abs': 740 716 self.write_c("fabs ") 741 elif (node.func.id == 'int'):717 elif node.func.id == 'int': 742 718 self.write_c('(int) ') 743 elif (node.func.id == "SINCOS"):719 elif node.func.id == "SINCOS": 744 720 self.WriteSincos(node) 745 721 return … … 772 748 def visit_Name(self, node): 773 749 self.write_c(node.id) 774 if ((node.id in self.C_Pointers) and(not self.SubRef)):750 if node.id in self.C_Pointers and not self.SubRef: 775 751 self.write_c("[0]") 776 752 name = "" 777 753 sub = node.id.find("[") 778 if (sub > 0):754 if sub > 0: 779 755 name = node.id[0:sub].strip() 780 756 else: 781 757 name = node.id 782 #add variable to C_Vars if it ins't there yet, not an argument and not a number783 if ( (name not in self.C_Functions) and (name not in self.C_Vars) and \784 (name not in self.C_IntVars) and (name not in self.arguments) and \785 (name not in self.C_Constants) and (name.isnumeric() == False)):786 if (self.InSubscript):758 # add variable to C_Vars if it ins't there yet, not an argument and not a number 759 if (name not in self.C_Functions and name not in self.C_Vars and 760 name not in self.C_IntVars and name not in self.arguments and 761 name not in self.C_Constants and not name.isdigit()): 762 if self.InSubscript: 787 763 self.C_IntVars.append(node.id) 788 764 else: … … 810 786 s = "" 811 787 for idx, item in enumerate(node.elts): 812 if ((idx > 0) and(len(s) > 0)):788 if idx > 0 and s: 813 789 s += ', ' 814 if (hasattr(item, 'id')):790 if hasattr(item, 'id'): 815 791 s += item.id 816 elif (hasattr(item, 'n')):792 elif hasattr(item, 'n'): 817 793 s += str(item.n) 818 if (len(s) > 0):794 if s: 819 795 self.C_Vectors.append(s) 820 796 vec_name = "vec" + str(len(self.C_Vectors)) 821 797 self.write_c(vec_name) 822 vec_name += "#"823 798 return visit 824 799 … … 840 815 function_name = '' 841 816 is_negative_exp = False 842 if (isevaluable(str(self.current_statement))):817 if isevaluable(str(self.current_statement)): 843 818 exponent = eval(string) 844 819 is_negative_exp = exponent < 0 845 820 abs_exponent = abs(exponent) 846 if (abs_exponent == 2):821 if abs_exponent == 2: 847 822 function_name = "square" 848 elif (abs_exponent == 3):823 elif abs_exponent == 3: 849 824 function_name = "cube" 850 elif (abs_exponent == 0.5):825 elif abs_exponent == 0.5: 851 826 function_name = "sqrt" 852 elif (abs_exponent == 1.0/3.0):827 elif abs_exponent == 1.0/3.0: 853 828 function_name = "cbrt" 854 if (function_name == ''):829 if function_name == '': 855 830 function_name = "pow" 856 831 return function_name, is_negative_exp 857 832 858 833 def translate_power(self, node): 859 # get exponent by visiting the right hand argument.834 # get exponent by visiting the right hand argument. 860 835 function_name = "pow" 861 836 temp_statement = self.current_statement 862 # 'visit' functions write the results to the 'current_statement' class memnber863 # Here, a temporary variable, 'temp_statement', is used, that enables the864 # use of the 'visit' function837 # 'visit' functions write the results to the 'current_statement' class memnber 838 # Here, a temporary variable, 'temp_statement', is used, that enables the 839 # use of the 'visit' function 865 840 self.current_statement = '' 866 841 self.visit(node.right) … … 868 843 function_name, is_negative_exp = self.get_special_power(self.current_statement) 869 844 self.current_statement = temp_statement 870 if (is_negative_exp):845 if is_negative_exp: 871 846 self.write_c("1.0 /(") 872 847 self.write_c(function_name + "(") 873 848 self.visit(node.left) 874 if (function_name == "pow"):849 if function_name == "pow": 875 850 self.write_c(", ") 876 851 self.visit(node.right) 877 852 self.write_c(")") 878 if (is_negative_exp):853 if is_negative_exp: 879 854 self.write_c(")") 880 855 self.write_c(" ") … … 889 864 def visit_BinOp(self, node): 890 865 self.write_c("(") 891 if ('%s' % BINOP_SYMBOLS[type(node.op)] == BINOP_SYMBOLS[ast.Pow]):866 if '%s' % BINOP_SYMBOLS[type(node.op)] == BINOP_SYMBOLS[ast.Pow]: 892 867 self.translate_power(node) 893 elif ('%s' % BINOP_SYMBOLS[type(node.op)] == BINOP_SYMBOLS[ast.FloorDiv]):868 elif '%s' % BINOP_SYMBOLS[type(node.op)] == BINOP_SYMBOLS[ast.FloorDiv]: 894 869 self.translate_integer_divide(node) 895 870 else: … … 899 874 self.write_c(")") 900 875 901 #for C876 # for C 902 877 def visit_BoolOp(self, node): 903 878 self.write_c('(') … … 926 901 927 902 def visit_Subscript(self, node): 928 if (node.value.id not in self.C_Constants):929 if (node.value.id not in self.C_Pointers):903 if node.value.id not in self.C_Constants: 904 if node.value.id not in self.C_Pointers: 930 905 self.C_Pointers.append(node.value.id) 931 906 self.SubRef = True … … 976 951 self.visit(comprehension) 977 952 self.write_c(right) 978 #self.write_python(right)953 #self.write_python(right) 979 954 return visit 980 955 … … 1005 980 1006 981 def visit_Repr(self, node): 1007 # XXX: python 2.6 only982 # CRUFT: python 2.6 only 1008 983 self.write_c('`') 1009 984 self.visit(node.value) … … 1021 996 self.visit(node.target) 1022 997 self.write_C(' in ') 1023 #self.write_python(' in ')998 #self.write_python(' in ') 1024 999 self.visit(node.iter) 1025 1000 if node.ifs: … … 1043 1018 print(tree_source) 1044 1019 1045 def add_constants(sniplets, c_constants):1046 sniplets.append("#include <math.h>")1047 sniplets.append("")1048 vars = c_constants.keys()1049 for c_var in vars:1050 c_values = c_constants[c_var]1051 if isinstance(c_values, (int, float)):1052 parts = ["double ", c_var, " = ", "%.15g"%c_values, ";"]1053 else:1054 elements = ["%.15g"%v for v in c_values]1055 parts = ["double ", c_var, "[]", " = ", "{\n ", ", ".join(elements), "\n};"]1056 sniplets.append("".join(parts))1057 1058 1020 def translate(functions, constants=0): 1059 snip lets = []1060 add_constants (sniplets, constants)1061 for source,fname,line_no in functions:1062 line_directive = '#line %d "%s"' %(line_no,fname)1063 line_directive = line_directive.replace('\\','\\\\')1064 # sniplets.append(line_directive)1021 snippets = [] 1022 #snippets.append("#include <math.h>") 1023 #snippets.append("") 1024 for source, fname, line_no in functions: 1025 line_directive = '#line %d "%s"'%(line_no, fname.replace('\\', '\\\\')) 1026 snippets.append(line_directive) 1065 1027 tree = ast.parse(source) 1066 sniplet = to_source(tree, functions, constants) # in the future add filename, offset, constants 1067 sniplets.append(sniplet) 1068 c_code = "\n".join(sniplets) 1069 f_out = open ("xlate.c", "w+") 1070 f_out.write (c_code) 1071 f_out.close() 1072 return("\n".join(sniplets)) 1073 1074 def get_file_names(): 1075 fname_in = "" 1076 fname_out = "" 1077 if(len(sys.argv) > 1): 1078 fname_in = sys.argv[1] 1079 fname_base = os.path.splitext(fname_in) 1080 if(len(sys.argv) == 2): 1081 fname_out = str(fname_base[0]) + '.c' 1082 else: 1083 fname_out = sys.argv[2] 1084 if(len(fname_in) > 0): 1085 python_file = open(sys.argv[1], "r") 1086 if(len(fname_out) > 0): 1087 file_out = open(fname_out, "w+") 1088 return len(sys.argv), fname_in, fname_out 1089 1090 if __name__ == "__main__": 1028 # in the future add filename, offset, constants 1029 c_code = to_source(tree, functions, constants) 1030 snippets.append(c_code) 1031 return snippets 1032 1033 def main(): 1091 1034 import os 1092 1035 print("Parsing...using Python" + sys.version) 1093 try: 1094 fname_in = "" 1095 fname_out = "" 1096 if(len(sys.argv) == 1): 1097 print("Usage:\npython parse01.py <infile> [<outfile>](if omitted, output file is '<infile>.c'") 1098 else: 1099 fname_in = sys.argv[1] 1100 fname_base = os.path.splitext(fname_in) 1101 if(len(sys.argv) == 2): 1102 fname_out = str(fname_base[0]) + '.c' 1103 else: 1104 fname_out = sys.argv[2] 1105 if(len(fname_in) > 0): 1106 python_file = open(sys.argv[1], "r") 1107 if(len(fname_out) > 0): 1108 file_out = open(fname_out, "w+") 1109 functions = ["MultAsgn", "Iq41", "Iq2"] 1110 tpls = [functions, fname_in, 0] 1111 c_txt = translate(tpls) 1112 file_out.write(c_txt) 1113 file_out.close() 1114 except Exception as excp: 1115 print("Error:\n" + str(excp.args)) 1036 if len(sys.argv) == 1: 1037 print("""\ 1038 Usage: python py2c.py <infile> [<outfile>] 1039 1040 if outfile is omitted, output file is '<infile>.c' 1041 """) 1042 return 1043 1044 fname_in = sys.argv[1] 1045 if len(sys.argv) == 2: 1046 fname_base = os.path.splitext(fname_in)[0] 1047 fname_out = str(fname_base) + '.c' 1048 else: 1049 fname_out = sys.argv[2] 1050 1051 with open(fname_in, "r") as python_file: 1052 code = python_file.read() 1053 1054 translation = translate([code, fname_in, 1])[0] 1055 1056 with open(fname_out, "w") as file_out: 1057 file_out.write(translation) 1116 1058 print("...Done") 1059 1060 if __name__ == "__main__": 1061 main()
Note: See TracChangeset
for help on using the changeset viewer.