Changes in / [1941ec6:8224d24] in sasmodels


Ignore:
Files:
2 deleted
43 edited

Legend:

Unmodified
Added
Removed

  • .gitignore

    r2badeca r2badeca  
    2525/example/Fit_*/ 
    2626/example/batch_fit.csv 
    27 /sasmodels/models/lib/gauss*.c 

  • doc/guide/plugin.rst

    rc654160 r0a9fcab  
    292292**Note: The order of the parameters in the definition will be the order of the 
    293293parameters in the user interface and the order of the parameters in Iq(), 
    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.** 
     294Iqxy() and form_volume(). And** *scale* **and** *background* **parameters are 
     295implicit to all models, so they do not need to be included in the parameter table.** 
    297296 
    298297- **"name"** is the name of the parameter shown on the FitPage. 
    … …  
    363362    scattered intensity. 
    364363 
    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. 
     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. 
    372369 
    373370Some models will have integer parameters, such as number of pearls in the 
    … …  
    422419That is, the individual models do not need to include polydispersity 
    423420calculations, but instead rely on numerical integration to compute the 
    424 appropriately smeared pattern. 
     421appropriately smeared pattern.   Angular dispersion values over polar angle 
     422$\theta$ requires an additional $\cos \theta$ weighting due to decreased 
     423arc length for the equatorial angle $\phi$ with increasing latitude. 
    425424 
    426425Python Models 
    … …  
    469468barbell).  If I(q; pars) is NaN for any $q$, then those parameters will be 
    470469ignored, and not included in the calculation of the weighted polydispersity. 
     470 
     471Similar to *Iq*, you can define *Iqxy(qx, qy, par1, par2, ...)* where the 
     472parameter list includes any orientation parameters.  If *Iqxy* is not defined, 
     473then it will default to *Iqxy = Iq(sqrt(qx**2+qy**2), par1, par2, ...)*. 
    471474 
    472475Models should define *form_volume(par1, par2, ...)* where the parameter 
    … …  
    494497    } 
    495498 
     499*Iqxy* is similar to *Iq*, except it uses parameters *qx, qy* instead of *q*, 
     500and it includes orientation parameters. 
     501 
    496502*form_volume* defines the volume of the shape. As in python models, it 
    497503includes only the volume parameters. 
    498504 
     505*Iqxy* will default to *Iq(sqrt(qx**2 + qy**2), par1, ...)* and 
     506*form_volume* will default to 1.0. 
     507 
    499508**source=['fn.c', ...]** includes the listed C source files in the 
    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. 
     509program before *Iq* and *Iqxy* are defined. This allows you to extend the 
     510library of C functions available to your model. 
    508511 
    509512Models are defined using double precision declarations for the 
    … …  
    529532 
    530533    #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, ...)*. 
    655534 
    656535Special Functions 
    … …  
    917796show a 50x improvement or more over the equivalent pure python model. 
    918797 
     798External C Models 
     799................. 
     800 
     801External C models are very much like embedded C models, except that 
     802*Iq*, *Iqxy* and *form_volume* are defined in an external source file 
     803loaded using the *source=[...]* statement. You need to supply the function 
     804declarations for each of these that you need instead of building them 
     805automatically from the parameter table. 
     806 
    919807 
    920808.. _Form_Factors: 
    … …  
    11181006variable name *Rg* for example because $R_g$ is the right name for the model 
    11191007parameter then ignore the lint errors.  Also, ignore *missing-docstring* 
    1120 for standard model functions *Iq*, *Iqac*, etc. 
     1008for standard model functions *Iq*, *Iqxy*, etc. 
    11211009 
    11221010We will have delinting sessions at the SasView Code Camps, where we can 

  • explore/asymint.py

    ra1c32c2 r1820208  
    8686    a, b, c = env.mpf(a), env.mpf(b), env.mpf(c) 
    8787    def Fq(qa, qb, qc): 
    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) 
     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) 
    9191        return siA * siB * siC 
    9292    Fq.__doc__ = "parallelepiped a=%g, b=%g c=%g"%(a, b, c) 
    9393    volume = a*b*c 
    9494    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) 
    12895    return norm, Fq 
    12996 
    … …  
    217184    NORM, KERNEL = make_parallelepiped(A, B, C) 
    218185    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) 
    237186elif shape == 'paracrystal': 
    238187    LATTICE = 'bcc' 
    … …  
    393342    print("gauss-150", *gauss_quad_2d(Q, n=150)) 
    394343    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)) 
    397344    #gridded_2d(Q, n=2**8+1) 
    398345    gridded_2d(Q, n=2**10+1) 
    399     #gridded_2d(Q, n=2**12+1) 
     346    #gridded_2d(Q, n=2**13+1) 
    400347    #gridded_2d(Q, n=2**15+1) 
    401     if shape not in ('paracrystal', 'core_shell_parallelepiped'): 
    402         # adaptive forms on models for which the calculations are fast enough 
     348    if shape != 'paracrystal':  # adaptive forms are too slow! 
    403349        print("dblquad", *scipy_dblquad_2d(Q)) 
    404350        print("semi-romberg-100", *semi_romberg_2d(Q, n=100)) 

  • sasmodels/compare.py

    r2d81cfe r2d81cfe  
    4242from .data import plot_theory, empty_data1D, empty_data2D, load_data 
    4343from .direct_model import DirectModel, get_mesh 
    44 from .generate import FLOAT_RE, set_integration_size 
     44from .generate import FLOAT_RE 
    4545from .weights import plot_weights 
    4646 
    … …  
    693693        data = empty_data2D(q, resolution=res) 
    694694        data.accuracy = opts['accuracy'] 
    695         set_beam_stop(data, qmin) 
     695        set_beam_stop(data, 0.0004) 
    696696        index = ~data.mask 
    697697    else: 
    … …  
    706706    return data, index 
    707707 
    708 def make_engine(model_info, data, dtype, cutoff, ngauss=0): 
     708def make_engine(model_info, data, dtype, cutoff): 
    709709    # type: (ModelInfo, Data, str, float) -> Calculator 
    710710    """ 
    … …  
    714714    than OpenCL. 
    715715    """ 
    716     if ngauss: 
    717         set_integration_size(model_info, ngauss) 
    718  
    719716    if dtype is None or not dtype.endswith('!'): 
    720717        return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff) 
    … …  
    957954    'poly', 'mono', 'cutoff=', 
    958955    'magnetic', 'nonmagnetic', 
    959     'accuracy=', 'ngauss=', 
     956    'accuracy=', 
    960957    'neval=',  # for timing... 
    961958 
    … …  
    10921089        'show_weights' : False, 
    10931090        'sphere'    : 0, 
    1094         'ngauss'    : '0', 
    10951091    } 
    10961092    for arg in flags: 
    … …  
    11191115        elif arg.startswith('-engine='):   opts['engine'] = arg[8:] 
    11201116        elif arg.startswith('-neval='):    opts['count'] = arg[7:] 
    1121         elif arg.startswith('-ngauss='):   opts['ngauss'] = arg[8:] 
    11221117        elif arg.startswith('-random='): 
    11231118            opts['seed'] = int(arg[8:]) 
    … …  
    11741169 
    11751170    comparison = any(PAR_SPLIT in v for v in values) 
    1176  
    11771171    if PAR_SPLIT in name: 
    11781172        names = name.split(PAR_SPLIT, 2) 
    … …  
    11871181        return None 
    11881182 
    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  
    11951183    if PAR_SPLIT in opts['engine']: 
    11961184        opts['engine'] = opts['engine'].split(PAR_SPLIT, 2) 
    … …  
    12111199        opts['cutoff'] = [float(opts['cutoff'])]*2 
    12121200 
    1213     base = make_engine(model_info[0], data, opts['engine'][0], 
    1214                        opts['cutoff'][0], opts['ngauss'][0]) 
     1201    base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0]) 
    12151202    if comparison: 
    1216         comp = make_engine(model_info[1], data, opts['engine'][1], 
    1217                            opts['cutoff'][1], opts['ngauss'][1]) 
     1203        comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1]) 
    12181204    else: 
    12191205        comp = None 
    … …  
    12881274        if model_info != model_info2: 
    12891275            pars2 = randomize_pars(model_info2, pars2) 
    1290             limit_dimensions(model_info2, pars2, maxdim) 
     1276            limit_dimensions(model_info, pars2, maxdim) 
    12911277            # Share values for parameters with the same name 
    12921278            for k, v in pars.items(): 

  • sasmodels/details.py

    r108e70e r2d81cfe  
    258258    # type: (...) -> Sequence[np.ndarray] 
    259259    """ 
    260     **Deprecated** Theta weights will be computed in the kernel wrapper if 
    261     they are needed. 
    262  
    263260    If there is a theta parameter, update the weights of that parameter so that 
    264261    the cosine weighting required for polar integration is preserved. 
    … …  
    275272    Returns updated weights vectors 
    276273    """ 
     274    # TODO: explain in a comment why scale and background are missing 
    277275    # Apparently the parameters.theta_offset similarly skips scale and 
    278276    # and background, so the indexing works out, but they are still shipped 
    … …  
    281279        index = parameters.theta_offset 
    282280        theta = dispersity[index] 
     281        # TODO: modify the dispersity vector to avoid the theta=-90,90,270,... 
    283282        theta_weight = abs(cos(radians(theta))) 
    284283        weights = tuple(theta_weight*w if k == index else w 

  • sasmodels/generate.py

    r108e70e rdb03406  
    77    particular dimensions averaged over all orientations. 
    88 
    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. 
     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. 
    2215 
    2316    *form_volume(p1, p2, ...)* returns the volume of the form with particular 
    … …  
    3831scale and background parameters for each model. 
    3932 
    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. 
     33*Iq*, *Iqxy*, *Imagnetic* and *form_volume* should be stylized C-99 
     34functions written for OpenCL.  All functions need prototype declarations 
     35even if the are defined before they are used.  OpenCL does not support 
     36*#include* preprocessor directives, so instead the list of includes needs 
     37to be given as part of the metadata in the kernel module definition. 
     38The included files should be listed using a path relative to the kernel 
     39module, or if using "lib/file.c" if it is one of the standard includes 
     40provided with the sasmodels source.  The includes need to be listed in 
     41order so that functions are defined before they are used. 
    4842 
    4943Floating point values should be declared as *double*.  For single precision 
    … …  
    113107    present, the volume ratio is 1. 
    114108 
    115     *form_volume*, *Iq*, *Iqac*, *Iqabc* are strings containing 
    116     the C source code for the body of the volume, Iq, and Iqac functions 
     109    *form_volume*, *Iq*, *Iqxy*, *Imagnetic* are strings containing the 
     110    C source code for the body of the volume, Iq, and Iqxy functions 
    117111    respectively.  These can also be defined in the last source file. 
    118112 
    119     *Iq*, *Iqac*, *Iqabc* also be instead be python functions defining the 
     113    *Iq* and *Iqxy* also be instead be python functions defining the 
    120114    kernel.  If they are marked as *Iq.vectorized = True* then the 
    121115    kernel is passed the entire *q* vector at once, otherwise it is 
    … …  
    174168from zlib import crc32 
    175169from inspect import currentframe, getframeinfo 
    176 import logging 
    177170 
    178171import numpy as np  # type: ignore 
    … …  
    188181    pass 
    189182# pylint: enable=unused-import 
    190  
    191 logger = logging.getLogger(__name__) 
    192183 
    193184# jitter projection to use in the kernel code.  See explore/jitter.py 
    … …  
    279270""" 
    280271 
    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] 
    299272 
    300273def format_units(units): 
    … …  
    635608 
    636609""" 
    637 def _gen_fn(model_info, name, pars): 
    638     # type: (ModelInfo, str, List[Parameter]) -> str 
     610def _gen_fn(name, pars, body, filename, line): 
     611    # type: (str, List[Parameter], str, str, int) -> str 
    639612    """ 
    640613    Generate a function given pars and body. 
    … …  
    648621    """ 
    649622    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) 
    654623    return _FN_TEMPLATE % { 
    655624        'name': name, 'pars': par_decl, 'body': body, 
    656         'filename': filename.replace('\\', '\\\\'), 'line': lineno, 
     625        'filename': filename.replace('\\', '\\\\'), 'line': line, 
    657626    } 
    658627 
    … …  
    669638 
    670639# type in IQXY pattern could be single, float, double, long double, ... 
    671 _IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]", 
     640_IQXY_PATTERN = re.compile("^((inline|static) )? *([a-z ]+ )? *Iqxy *([(]|$)", 
    672641                           flags=re.MULTILINE) 
    673 def find_xy_mode(source): 
     642def _have_Iqxy(sources): 
    674643    # type: (List[str]) -> bool 
    675644    """ 
    676     Return the xy mode as qa, qac, qabc or qxy. 
     645    Return true if any file defines Iqxy. 
    677646 
    678647    Note this is not a C parser, and so can be easily confused by 
    679648    non-standard syntax.  Also, it will incorrectly identify the following 
    680     as having 2D models:: 
     649    as having Iqxy:: 
    681650 
    682651        /* 
    683         double Iqac(qab, qc, ...) { ... fill this in later ... } 
     652        double Iqxy(qx, qy, ...) { ... fill this in later ... } 
    684653        */ 
    685654 
    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): 
     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 
     664def _add_source(source, code, path): 
    702665    """ 
    703666    Add a file to the list of source code chunks, tagged with path and line. 
    704667    """ 
    705668    path = path.replace('\\', '\\\\') 
    706     source.append('#line %d "%s"' % (lineno, path)) 
     669    source.append('#line 1 "%s"' % path) 
    707670    source.append(code) 
    708671 
    … …  
    735698    user_code = [(f, open(f).read()) for f in model_sources(model_info)] 
    736699 
     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_asymmetric 
     703               else 'qabc') 
     704 
    737705    # Build initial sources 
    738706    source = [] 
    … …  
    742710 
    743711    if model_info.c_code: 
    744         _add_source(source, model_info.c_code, model_info.filename, 
    745                     lineno=model_info.lineno.get('c_code', 1)) 
     712        source.append(model_info.c_code) 
    746713 
    747714    # Make parameters for q, qx, qy so that we can use them in declarations 
    748     q, qx, qy, qab, qa, qb, qc \ 
    749         = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()] 
     715    q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 
    750716    # Generate form_volume function, etc. from body only 
    751717    if isinstance(model_info.form_volume, str): 
    752718        pars = partable.form_volume_parameters 
    753         source.append(_gen_fn(model_info, 'form_volume', pars)) 
     719        source.append(_gen_fn('form_volume', pars, model_info.form_volume, 
     720                              model_info.filename, model_info._form_volume_line)) 
    754721    if isinstance(model_info.Iq, str): 
    755722        pars = [q] + partable.iq_parameters 
    756         source.append(_gen_fn(model_info, 'Iq', pars)) 
     723        source.append(_gen_fn('Iq', pars, model_info.Iq, 
     724                              model_info.filename, model_info._Iq_line)) 
    757725    if isinstance(model_info.Iqxy, str): 
    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") 
     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)) 
    780729 
    781730    # Define the parameter table 
    … …  
    803752    if xy_mode == 'qabc': 
    804753        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
    805         call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqabc(%s)" % pars 
     754        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars 
    806755        clear_iqxy = "#undef CALL_IQ_ABC" 
    807756    elif xy_mode == 'qac': 
    808757        pars = ",".join(["_qa", "_qc"] + model_refs) 
    809         call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars 
     758        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars 
    810759        clear_iqxy = "#undef CALL_IQ_AC" 
    811     elif xy_mode == 'qa': 
     760    else:  # xy_mode == 'qa' 
    812761        pars = ",".join(["_qa"] + model_refs) 
    813762        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
    814763        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  
    827764 
    828765    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 

  • sasmodels/kernel_header.c

    r108e70e r8698a0d  
    150150inline double cube(double x) { return x*x*x; } 
    151151inline 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

    r108e70e r6aee3ab  
    3131//  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes 
    3232//  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 
    3433//  INVALID(table) : test if the current point is feesible to calculate.  This 
    3534//      will be defined in the kernel definition file. 
    … …  
    470469  #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc) 
    471470  #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) 
    479471#endif 
    480472 
    … …  
    485477#if defined(CALL_IQ) || defined(CALL_IQ_A) 
    486478  #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 
    517479#else // !spherosymmetric projection 
    518480  // Grab the "view" angles (theta, phi, psi) from the initial parameter table. 

  • sasmodels/kernelpy.py

    r108e70e r2d81cfe  
    2626# pylint: enable=unused-import 
    2727 
    28 logger = logging.getLogger(__name__) 
    29  
    3028class PyModel(KernelModel): 
    3129    """ 
    … …  
    3331    """ 
    3432    def __init__(self, model_info): 
    35         # Make sure Iq is available and vectorized 
     33        # Make sure Iq and Iqxy are available and vectorized 
    3634        _create_default_functions(model_info) 
    3735        self.info = model_info 
    3836        self.dtype = np.dtype('d') 
    39         logger.info("load python model " + self.info.name) 
     37        logging.info("load python model " + self.info.name) 
    4038 
    4139    def make_kernel(self, q_vectors): 
    4240        q_input = PyInput(q_vectors, dtype=F64) 
    43         return PyKernel(self.info, q_input) 
     41        kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq 
     42        return PyKernel(kernel, self.info, q_input) 
    4443 
    4544    def release(self): 
    … …  
    9089    Callable SAS kernel. 
    9190 
    92     *kernel* is the kernel object to call. 
     91    *kernel* is the DllKernel object to call. 
    9392 
    9493    *model_info* is the module information 
    … …  
    105104    Call :meth:`release` when done with the kernel instance. 
    106105    """ 
    107     def __init__(self, model_info, q_input): 
     106    def __init__(self, kernel, model_info, q_input): 
    108107        # type: (callable, ModelInfo, List[np.ndarray]) -> None 
    109108        self.dtype = np.dtype('d') 
    … …  
    111110        self.q_input = q_input 
    112111        self.res = np.empty(q_input.nq, q_input.dtype) 
     112        self.kernel = kernel 
    113113        self.dim = '2d' if q_input.is_2d else '1d' 
    114114 
    … …  
    159159        # Generate a closure which calls the form_volume if it exists. 
    160160        form_volume = model_info.form_volume 
    161         self._volume = ((lambda: form_volume(*volume_args)) if form_volume else 
    162                         (lambda: 1.0)) 
     161        self._volume = ((lambda: form_volume(*volume_args)) if form_volume 
     162                        else (lambda: 1.0)) 
    163163 
    164164    def __call__(self, call_details, values, cutoff, magnetic): 
    … …  
    261261    any functions that are not already marked as vectorized. 
    262262    """ 
    263     # Note: must call create_vector_Iq before create_vector_Iqxy 
    264263    _create_vector_Iq(model_info) 
    265     _create_vector_Iqxy(model_info) 
     264    _create_vector_Iqxy(model_info)  # call create_vector_Iq() first 
    266265 
    267266 
    … …  
    281280        model_info.Iq = vector_Iq 
    282281 
    283  
    284282def _create_vector_Iqxy(model_info): 
    285283    """ 
    286284    Define Iqxy as a vector function if it exists, or default it from Iq(). 
    287285    """ 
    288     Iqxy = getattr(model_info, 'Iqxy', None) 
     286    Iq, Iqxy = model_info.Iq, model_info.Iqxy 
    289287    if callable(Iqxy): 
    290288        if not getattr(Iqxy, 'vectorized', False): 
    … …  
    297295            vector_Iqxy.vectorized = True 
    298296            model_info.Iqxy = vector_Iqxy 
    299     else: 
     297    elif callable(Iq): 
    300298        #print("defaulting Iqxy") 
    301299        # Iq is vectorized because create_vector_Iq was already called. 
    302         Iq = model_info.Iq 
    303300        def default_Iqxy(qx, qy, *args): 
    304301            """ 

  • sasmodels/modelinfo.py

    r108e70e rdb03406  
    4242 
    4343# assumptions about common parameters exist throughout the code, such as: 
    44 # (1) kernel functions Iq, Iqxy, Iqac, Iqabc, form_volume, ... don't see them 
     44# (1) kernel functions Iq, Iqxy, form_volume, ... don't see them 
    4545# (2) kernel drivers assume scale is par[0] and background is par[1] 
    4646# (3) mixture models drop the background on components and replace the scale 
    … …  
    261261 
    262262    *type* indicates how the parameter will be used.  "volume" parameters 
    263     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 
     263    will be used in all functions.  "orientation" parameters will be used 
     264    in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in 
     265    *Imagnetic* only.  If *type* is the empty string, the parameter will 
    266266    be used in all of *Iq*, *Iqxy* and *Imagnetic*.  "sld" parameters 
    267267    can automatically be promoted to magnetic parameters, each of which 
    … …  
    391391      with vector parameter p sent as p[]. 
    392392 
     393    * [removed] *iqxy_parameters* is the list of parameters to the Iqxy(qx, qy, ...) 
     394      function, with vector parameter p sent as p[]. 
     395 
    393396    * *form_volume_parameters* is the list of parameters to the form_volume(...) 
    394397      function, with vector parameter p sent as p[]. 
    … …  
    445448        self.iq_parameters = [p for p in self.kernel_parameters 
    446449                              if p.type not in ('orientation', 'magnetic')] 
    447         self.orientation_parameters = [p for p in self.kernel_parameters 
    448                                        if p.type == 'orientation'] 
     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'] 
    449453        self.form_volume_parameters = [p for p in self.kernel_parameters 
    450454                                       if p.type == 'volume'] 
    … …  
    491495                if p.type != 'orientation': 
    492496                    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") 
    495497        if theta >= 0 and phi >= 0: 
    496             last_par = len(self.kernel_parameters) - 1 
    497498            if phi != theta+1: 
    498499                raise TypeError("phi must follow theta") 
    499500            if psi >= 0 and psi != phi+1: 
    500501                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") 
    504502        elif theta >= 0 or phi >= 0 or psi >= 0: 
    505503            raise TypeError("oriented shapes must have both theta and phi and maybe psi") 
    … …  
    722720 
    723721 
    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  
    727722def _find_source_lines(model_info, kernel_module): 
    728723    # type: (ModelInfo, ModuleType) -> None 
    … …  
    730725    Identify the location of the C source inside the model definition file. 
    731726 
    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)): 
     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): 
    745743        return 
    746744 
    747     # load the module source if we can 
     745    # find the defintion lines for the different code blocks 
    748746    try: 
    749747        source = inspect.getsource(kernel_module) 
    750748    except IOError: 
    751749        return 
    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 
     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 
    760760 
    761761def make_model_info(kernel_module): 
    … …  
    766766    Fill in default values for parts of the module that are not provided. 
    767767 
    768     Note: vectorized Iq and Iqac/Iqabc functions will be created for python 
     768    Note: vectorized Iq and Iqxy functions will be created for python 
    769769    models when the model is first called, not when the model is loaded. 
    770770    """ 
    … …  
    796796    info.c_code = getattr(kernel_module, 'c_code', None) 
    797797    info.source = getattr(kernel_module, 'source', []) 
    798     info.c_code = getattr(kernel_module, 'c_code', None) 
    799798    # TODO: check the structure of the tests 
    800799    info.tests = getattr(kernel_module, 'tests', []) 
    … …  
    804803    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    805804    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 
    808805    info.Imagnetic = getattr(kernel_module, 'Imagnetic', None) # type: ignore 
    809806    info.profile = getattr(kernel_module, 'profile', None) # type: ignore 
    … …  
    820817    info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore 
    821818 
    822     if callable(info.Iq) and parameters.has_2d: 
    823         raise ValueError("oriented python models not supported") 
    824  
    825     info.lineno = {} 
    826819    _find_source_lines(info, kernel_module) 
    827820    try: 
    … …  
    842835 
    843836    The structure should be mostly static, other than the delayed definition 
    844     of *Iq*, *Iqac* and *Iqabc* if they need to be defined. 
     837    of *Iq* and *Iqxy* if they need to be defined. 
    845838    """ 
    846839    #: Full path to the file defining the kernel, if any. 
    … …  
    924917    structure_factor = None # type: bool 
    925918    #: List of C source files used to define the model.  The source files 
    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. 
     919    #: should define the *Iq* function, and possibly *Iqxy*, though a default 
     920    #: *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 are 
     923    #: indicated by providing a :attr:`ER` function. 
    931924    source = None           # type: List[str] 
    932925    #: The set of tests that must pass.  The format of the tests is described 
    … …  
    977970    #: include the decimal point. See :mod:`generate` for more details. 
    978971    Iq = 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]] 
     972    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`. 
     973    Iqxy = None             # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
    983974    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`. 
    984975    Imagnetic = None        # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
    … …  
    996987    #: Returns a random parameter set for the model 
    997988    random = None           # type: Optional[Callable[[], Dict[str, float]]] 
    998     #: Line numbers for symbols defining C code 
    999     lineno = None           # type: Dict[str, int] 
     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 
    1000997 
    1001998    def __init__(self): 

  • sasmodels/models/_spherepy.py

    r108e70e ref07e95  
    8888Iq.vectorized = True  # Iq accepts an array of q values 
    8989 
     90def Iqxy(qx, qy, sld, sld_solvent, radius): 
     91    return Iq(sqrt(qx ** 2 + qy ** 2), sld, sld_solvent, radius) 
     92Iqxy.vectorized = True  # Iqxy accepts arrays of qx, qy values 
     93 
    9094def sesans(z, sld, sld_solvent, radius): 
    9195    """ 

  • sasmodels/models/barbell.c

    r108e70e rbecded3  
    2323    const double qab_r = radius_bell*qab; // Q*R*sin(theta) 
    2424    double total = 0.0; 
    25     for (int i = 0; i < GAUSS_N; i++){ 
    26         const double t = GAUSS_Z[i]*zm + zb; 
     25    for (int i = 0; i < 76; i++){ 
     26        const double t = Gauss76Z[i]*zm + zb; 
    2727        const double radical = 1.0 - t*t; 
    2828        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    2929        const double Fq = cos(m*t + b) * radical * bj; 
    30         total += GAUSS_W[i] * Fq; 
     30        total += Gauss76Wt[i] * Fq; 
    3131    } 
    3232    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    7373    const double zb = M_PI_4; 
    7474    double total = 0.0; 
    75     for (int i = 0; i < GAUSS_N; i++){ 
    76         const double alpha= GAUSS_Z[i]*zm + zb; 
     75    for (int i = 0; i < 76; i++){ 
     76        const double alpha= Gauss76Z[i]*zm + zb; 
    7777        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    7878        SINCOS(alpha, sin_alpha, cos_alpha); 
    7979        const double Aq = _fq(q*sin_alpha, q*cos_alpha, h, radius_bell, radius, half_length); 
    80         total += GAUSS_W[i] * Aq * Aq * sin_alpha; 
     80        total += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
    8181    } 
    8282    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    9090 
    9191static double 
    92 Iqac(double qab, double qc, 
     92Iqxy(double qab, double qc, 
    9393    double sld, double solvent_sld, 
    9494    double radius_bell, double radius, double length) 

  • sasmodels/models/bcc_paracrystal.c

    r108e70e rea60e08  
    8181 
    8282    double outer_sum = 0.0; 
    83     for(int i=0; i<GAUSS_N; i++) { 
     83    for(int i=0; i<150; i++) { 
    8484        double inner_sum = 0.0; 
    85         const double theta = GAUSS_Z[i]*theta_m + theta_b; 
     85        const double theta = Gauss150Z[i]*theta_m + theta_b; 
    8686        double sin_theta, cos_theta; 
    8787        SINCOS(theta, sin_theta, cos_theta); 
    8888        const double qc = q*cos_theta; 
    8989        const double qab = q*sin_theta; 
    90         for(int j=0;j<GAUSS_N;j++) { 
    91             const double phi = GAUSS_Z[j]*phi_m + phi_b; 
     90        for(int j=0;j<150;j++) { 
     91            const double phi = Gauss150Z[j]*phi_m + phi_b; 
    9292            double sin_phi, cos_phi; 
    9393            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    9595            const double qb = qab*sin_phi; 
    9696            const double form = bcc_Zq(qa, qb, qc, dnn, d_factor); 
    97             inner_sum += GAUSS_W[j] * form; 
     97            inner_sum += Gauss150Wt[j] * form; 
    9898        } 
    9999        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    100         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     100        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
    101101    } 
    102102    outer_sum *= theta_m; 
    … …  
    107107 
    108108 
    109 static double Iqabc(double qa, double qb, double qc, 
     109static double Iqxy(double qa, double qb, double qc, 
    110110    double dnn, double d_factor, double radius, 
    111111    double sld, double solvent_sld) 

  • sasmodels/models/capped_cylinder.c

    r108e70e rbecded3  
    3030    const double qab_r = radius_cap*qab; // Q*R*sin(theta) 
    3131    double total = 0.0; 
    32     for (int i=0; i<GAUSS_N; i++) { 
    33         const double t = GAUSS_Z[i]*zm + zb; 
     32    for (int i=0; i<76 ;i++) { 
     33        const double t = Gauss76Z[i]*zm + zb; 
    3434        const double radical = 1.0 - t*t; 
    3535        const double bj = sas_2J1x_x(qab_r*sqrt(radical)); 
    3636        const double Fq = cos(m*t + b) * radical * bj; 
    37         total += GAUSS_W[i] * Fq; 
     37        total += Gauss76Wt[i] * Fq; 
    3838    } 
    3939    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    9595    const double zb = M_PI_4; 
    9696    double total = 0.0; 
    97     for (int i=0; i<GAUSS_N ;i++) { 
    98         const double theta = GAUSS_Z[i]*zm + zb; 
     97    for (int i=0; i<76 ;i++) { 
     98        const double theta = Gauss76Z[i]*zm + zb; 
    9999        double sin_theta, cos_theta; // slots to hold sincos function output 
    100100        SINCOS(theta, sin_theta, cos_theta); 
    … …  
    103103        const double Aq = _fq(qab, qc, h, radius_cap, radius, half_length); 
    104104        // scale by sin_theta for spherical coord integration 
    105         total += GAUSS_W[i] * Aq * Aq * sin_theta; 
     105        total += Gauss76Wt[i] * Aq * Aq * sin_theta; 
    106106    } 
    107107    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    115115 
    116116static double 
    117 Iqac(double qab, double qc, 
     117Iqxy(double qab, double qc, 
    118118    double sld, double solvent_sld, double radius, 
    119119    double radius_cap, double length) 

  • sasmodels/models/core_shell_bicelle.c

    r108e70e rbecded3  
    5252 
    5353    double total = 0.0; 
    54     for(int i=0;i<GAUSS_N;i++) { 
    55         double theta = (GAUSS_Z[i] + 1.0)*uplim; 
     54    for(int i=0;i<N_POINTS_76;i++) { 
     55        double theta = (Gauss76Z[i] + 1.0)*uplim; 
    5656        double sin_theta, cos_theta; // slots to hold sincos function output 
    5757        SINCOS(theta, sin_theta, cos_theta); 
    5858        double fq = bicelle_kernel(q*sin_theta, q*cos_theta, radius, thick_radius, thick_face, 
    5959                                   halflength, sld_core, sld_face, sld_rim, sld_solvent); 
    60         total += GAUSS_W[i]*fq*fq*sin_theta; 
     60        total += Gauss76Wt[i]*fq*fq*sin_theta; 
    6161    } 
    6262 
    … …  
    6767 
    6868static double 
    69 Iqac(double qab, double qc, 
     69Iqxy(double qab, double qc, 
    7070    double radius, 
    7171    double thick_rim, 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    r108e70e r82592da  
    3737    //initialize integral 
    3838    double outer_sum = 0.0; 
    39     for(int i=0;i<GAUSS_N;i++) { 
     39    for(int i=0;i<76;i++) { 
    4040        //setup inner integral over the ellipsoidal cross-section 
    4141        //const double va = 0.0; 
    4242        //const double vb = 1.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; 
     43        //const double cos_theta = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
     44        const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0; 
    4545        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
    4646        const double qab = q*sin_theta; 
    … …  
    4949        const double si2 = sas_sinx_x((halfheight+thick_face)*qc); 
    5050        double inner_sum=0.0; 
    51         for(int j=0;j<GAUSS_N;j++) { 
     51        for(int j=0;j<76;j++) { 
    5252            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    5353            // inner integral limits 
    5454            //const double vaj=0.0; 
    5555            //const double vbj=M_PI; 
    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; 
     56            //const double phi = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     57            const double phi = ( Gauss76Z[j] +1.0)*M_PI_2; 
    5858            const double rr = sqrt(r2A - r2B*cos(phi)); 
    5959            const double be1 = sas_2J1x_x(rr*qab); 
    … …  
    6161            const double fq = dr1*si1*be1 + dr2*si2*be2 + dr3*si2*be1; 
    6262 
    63             inner_sum += GAUSS_W[j] * fq * fq; 
     63            inner_sum += Gauss76Wt[j] * fq * fq; 
    6464        } 
    6565        //now calculate outer integral 
    66         outer_sum += GAUSS_W[i] * inner_sum; 
     66        outer_sum += Gauss76Wt[i] * inner_sum; 
    6767    } 
    6868 
    … …  
    7171 
    7272static double 
    73 Iqabc(double qa, double qb, double qc, 
     73Iqxy(double qa, double qb, double qc, 
    7474    double r_minor, 
    7575    double x_core, 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r108e70e r82592da  
    77        double length) 
    88{ 
    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 +  
    1010                 square(r_minor)*x_core*2.0*thick_face  ); 
    1111} 
    … …  
    4747    //initialize integral 
    4848    double outer_sum = 0.0; 
    49     for(int i=0;i<GAUSS_N;i++) { 
     49    for(int i=0;i<76;i++) { 
    5050        //setup inner integral over the ellipsoidal cross-section 
    5151        // since we generate these lots of times, why not store them somewhere? 
    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; 
     52        //const double cos_alpha = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
     53        const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0; 
    5454        const double sin_alpha = sqrt(1.0 - cos_alpha*cos_alpha); 
    5555        double inner_sum=0; 
    … …  
    5858        si1 = sas_sinx_x(sinarg1); 
    5959        si2 = sas_sinx_x(sinarg2); 
    60         for(int j=0;j<GAUSS_N;j++) { 
     60        for(int j=0;j<76;j++) { 
    6161            //76 gauss points for the inner integral (WAS 20 points,so this may make unecessarily slow, but playing safe) 
    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; 
     62            //const double beta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     63            const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 
    6464            const double rr = sqrt(r2A - r2B*cos(beta)); 
    6565            double besarg1 = q*rr*sin_alpha; 
    … …  
    6767            be1 = sas_2J1x_x(besarg1); 
    6868            be2 = sas_2J1x_x(besarg2); 
    69             inner_sum += GAUSS_W[j] *square(dr1*si1*be1 + 
     69            inner_sum += Gauss76Wt[j] *square(dr1*si1*be1 + 
    7070                                              dr2*si1*be2 + 
    7171                                              dr3*si2*be1); 
    7272        } 
    7373        //now calculate outer integral 
    74         outer_sum += GAUSS_W[i] * inner_sum; 
     74        outer_sum += Gauss76Wt[i] * inner_sum; 
    7575    } 
    7676 
    … …  
    7979 
    8080static double 
    81 Iqabc(double qa, double qb, double qc, 
     81Iqxy(double qa, double qb, double qc, 
    8282          double r_minor, 
    8383          double x_core, 
    … …  
    114114    return 1.0e-4 * Aq*exp(-0.5*(square(qa) + square(qb) + square(qc) )*square(sigma)); 
    115115} 
     116 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.py

    r108e70e r110f69c  
    149149    ["sld_rim",        "1e-6/Ang^2", 1, [-inf, inf], "sld",         "Cylinder rim scattering length density"], 
    150150    ["sld_solvent",    "1e-6/Ang^2", 6, [-inf, inf], "sld",         "Solvent scattering length density"], 
    151     ["sigma",       "Ang",        0,    [0, inf],    "",            "interfacial roughness"], 
    152151    ["theta",       "degrees",    90.0, [-360, 360], "orientation", "cylinder axis to beam angle"], 
    153152    ["phi",         "degrees",    0,    [-360, 360], "orientation", "rotation about beam"], 
    154153    ["psi",         "degrees",    0,    [-360, 360], "orientation", "rotation about cylinder axis"], 
     154    ["sigma",       "Ang",        0,    [0, inf],    "",            "interfacial roughness"] 
    155155    ] 
    156156 

  • sasmodels/models/core_shell_cylinder.c

    r108e70e rbecded3  
    3030    const double shell_vd = form_volume(radius,thickness,length) * (shell_sld-solvent_sld); 
    3131    double total = 0.0; 
    32     for (int i=0; i<GAUSS_N ;i++) { 
     32    for (int i=0; i<76 ;i++) { 
    3333        // translate a point in [-1,1] to a point in [0, pi/2] 
    34         //const double theta = ( GAUSS_Z[i]*(upper-lower) + upper + lower )/2.0; 
     34        //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    3535        double sin_theta, cos_theta; 
    36         const double theta = GAUSS_Z[i]*M_PI_4 + M_PI_4; 
     36        const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4; 
    3737        SINCOS(theta, sin_theta,  cos_theta); 
    3838        const double qab = q*sin_theta; 
    … …  
    4040        const double fq = _cyl(core_vd, core_r*qab, core_h*qc) 
    4141            + _cyl(shell_vd, shell_r*qab, shell_h*qc); 
    42         total += GAUSS_W[i] * fq * fq * sin_theta; 
     42        total += Gauss76Wt[i] * fq * fq * sin_theta; 
    4343    } 
    4444    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    4848 
    4949 
    50 double Iqac(double qab, double qc, 
     50double Iqxy(double qab, double qc, 
    5151    double core_sld, 
    5252    double shell_sld, 

  • sasmodels/models/core_shell_ellipsoid.c

    r108e70e rbecded3  
    5959    const double b = 0.5; 
    6060    double total = 0.0;     //initialize intergral 
    61     for(int i=0;i<GAUSS_N;i++) { 
    62         const double cos_theta = GAUSS_Z[i]*m + b; 
     61    for(int i=0;i<76;i++) { 
     62        const double cos_theta = Gauss76Z[i]*m + b; 
    6363        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
    6464        double fq = _cs_ellipsoid_kernel(q*sin_theta, q*cos_theta, 
    … …  
    6666            equat_shell, polar_shell, 
    6767            sld_core_shell, sld_shell_solvent); 
    68         total += GAUSS_W[i] * fq * fq; 
     68        total += Gauss76Wt[i] * fq * fq; 
    6969    } 
    7070    total *= m; 
    … …  
    7575 
    7676static double 
    77 Iqac(double qab, double qc, 
     77Iqxy(double qab, double qc, 
    7878    double radius_equat_core, 
    7979    double x_core, 

  • sasmodels/models/core_shell_parallelepiped.c

    r108e70e rc69d6d6  
    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 
    71static double 
    82form_volume(double length_a, double length_b, double length_c, 
    93    double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    104{ 
    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 
     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; 
    2710} 
    2811 
    … …  
    4124    double thick_rim_c) 
    4225{ 
    43     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c 
     26    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
    4427    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    45     // Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
    46     // Code rewritten (PAK) 
     28    //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
    4729 
    48     const double half_q = 0.5*q; 
     30    const double mu = 0.5 * q * length_b; 
    4931 
    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; 
     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; 
    5338 
    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); 
     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; 
    5969 
    6070    // outer integral (with gauss points), integration limits = 0, 1 
    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); 
     71    double outer_total = 0; //initialize integral 
    6572 
    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); 
     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); 
    7876 
    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 
     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); 
    9088 
    91             inner_sum += GAUSS_W[j] * f * f; 
     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; 
    9296        } 
    93         inner_sum *= 0.5; 
     97        inner_total *= 0.5; 
     98        inner_total_crim *= 0.5; 
    9499        // now sum up the outer integral 
    95         outer_sum += GAUSS_W[i] * inner_sum; 
     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); 
    96103    } 
    97     outer_sum *= 0.5; 
     104    outer_total *= 0.5; 
    98105 
    99106    //convert from [1e-12 A-1] to [cm-1] 
    100     return 1.0e-4 * outer_sum; 
     107    return 1.0e-4 * outer_total; 
    101108} 
    102109 
    103110static double 
    104 Iqabc(double qa, double qb, double qc, 
     111Iqxy(double qa, double qb, double qc, 
    105112    double core_sld, 
    106113    double arim_sld, 
    … …  
    121128    const double drC = crim_sld-solvent_sld; 
    122129 
     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 used 
     135    //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 
    123140    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    124141    // the scaling by B. 
    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); 
     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); 
    134152 
    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 
     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); 
    146160 
    147161    return 1.0e-4 * f * f; 

  • sasmodels/models/core_shell_parallelepiped.py

    r10ee838 r2d81cfe  
    55Calculates the form factor for a rectangular solid with a core-shell structure. 
    66The thickness and the scattering length density of the shell or 
    7 "rim" can be different on each (pair) of faces. 
     7"rim" can be different on each (pair) of faces. However at this time the 1D 
     8calculation does **NOT** actually calculate a c face rim despite the presence 
     9of 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. 
    814 
    915The form factor is normalized by the particle volume $V$ such that 
    … …  
    1521where $\langle \ldots \rangle$ is an average over all possible orientations 
    1622of the rectangular solid. 
     23 
    1724 
    1825The function calculated is the form factor of the rectangular solid below. 
    … …  
    3441    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    3542 
    36 **meaning that there are "gaps" at the corners of the solid.** 
     43**meaning that there are "gaps" at the corners of the solid.**  Again note that 
     44$t_C = 0$ currently. 
    3745 
    3846The intensity calculated follows the :ref:`parallelepiped` model, with the 
    3947core-shell intensity being calculated as the square of the sum of the 
    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. 
     48amplitudes of the core and shell, in the same manner as a core-shell model. 
    4849 
    4950.. math:: 
    5051 
    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. 
     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. 
    8069 
    8170FITTING NOTES 
    82 ~~~~~~~~~~~~~ 
    83  
    8471If the scale is set equal to the particle volume fraction, $\phi$, the returned 
    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.** 
     72value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
     73However, **no interparticle interference effects are included in this 
     74calculation.** 
    8775 
    8876There are many parameters in this model. Hold as many fixed as possible with 
    8977known values, or you will certainly end up at a solution that is unphysical. 
    9078 
     79Constraints must be applied during fitting to ensure that the inequality 
     80$A < B < C$ is not violated. The calculation will not report an error, 
     81but the results will not be correct. 
     82 
    9183The returned value is in units of |cm^-1|, on absolute scale. 
    9284 
    9385NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
    9486based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    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. 
     87and length $(C+2t_C)$ values, after appropriately 
     88sorting the three dimensions to give an oblate or prolate particle, to give an 
     89effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    9890 
    9991For 2d data the orientation of the particle is required, described using 
    100 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 
    101 details of the calculation and angular dispersions see :ref:`orientation`. 
     92angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
     93of the calculation and angular dispersions see :ref:`orientation` . 
    10294The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 
    10395$\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. 
    10996 
    11097.. figure:: img/parallelepiped_angle_definition.png 
    … …  
    127114    Equations (1), (13-14). (in German) 
    128115.. [#] D Singh (2009). *Small angle scattering studies of self assembly in 
    129    lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available 
     116   lipid mixtures*, John's Hopkins University Thesis (2009) 223-225. `Available 
    130117   from Proquest <http://search.proquest.com/docview/304915826?accountid 
    131118   =26379>`_ 
    … …  
    188175        Return equivalent radius (ER) 
    189176    """ 
    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) 
     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.) 
    196185 
    197186# VR defaults to 1.0 
    … …  
    227216            psi_pd=10, psi_pd_n=1) 
    228217 
    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 
     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 
    231219if 0:  # pak: model rewrite; need to update tests 
    232220    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 

  • sasmodels/models/cylinder.c

    r108e70e rbecded3  
    2121 
    2222    double total = 0.0; 
    23     for (int i=0; i<GAUSS_N ;i++) { 
    24         const double theta = GAUSS_Z[i]*zm + zb; 
     23    for (int i=0; i<76 ;i++) { 
     24        const double theta = Gauss76Z[i]*zm + zb; 
    2525        double sin_theta, cos_theta; // slots to hold sincos function output 
    2626        // theta (theta,phi) the projection of the cylinder on the detector plane 
    2727        SINCOS(theta , sin_theta, cos_theta); 
    2828        const double form = fq(q*sin_theta, q*cos_theta, radius, length); 
    29         total += GAUSS_W[i] * form * form * sin_theta; 
     29        total += Gauss76Wt[i] * form * form * sin_theta; 
    3030    } 
    3131    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    4545 
    4646static double 
    47 Iqac(double qab, double qc, 
     47Iqxy(double qab, double qc, 
    4848    double sld, 
    4949    double solvent_sld, 

  • sasmodels/models/ellipsoid.c

    r108e70e rbecded3  
    2222 
    2323    // translate a point in [-1,1] to a point in [0, 1] 
    24     // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 
     24    // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2; 
    2525    const double zm = 0.5; 
    2626    const double zb = 0.5; 
    2727    double total = 0.0; 
    28     for (int i=0;i<GAUSS_N;i++) { 
    29         const double u = GAUSS_Z[i]*zm + zb; 
     28    for (int i=0;i<76;i++) { 
     29        const double u = Gauss76Z[i]*zm + zb; 
    3030        const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 
    3131        const double f = sas_3j1x_x(q*r); 
    32         total += GAUSS_W[i] * f * f; 
     32        total += Gauss76Wt[i] * f * f; 
    3333    } 
    3434    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    3939 
    4040static double 
    41 Iqac(double qab, double qc, 
     41Iqxy(double qab, double qc, 
    4242    double sld, 
    4343    double sld_solvent, 

  • sasmodels/models/elliptical_cylinder.c

    r108e70e r82592da  
    2222    //initialize integral 
    2323    double outer_sum = 0.0; 
    24     for(int i=0;i<GAUSS_N;i++) { 
     24    for(int i=0;i<76;i++) { 
    2525        //setup inner integral over the ellipsoidal cross-section 
    26         const double cos_val = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
     26        const double cos_val = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
    2727        const double sin_val = sqrt(1.0 - cos_val*cos_val); 
    2828        //const double arg = radius_minor*sin_val; 
    2929        double inner_sum=0; 
    30         for(int j=0;j<GAUSS_N;j++) { 
    31             const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     30        for(int j=0;j<76;j++) { 
     31            //20 gauss points for the inner integral, increase to 76, RKH 6Nov2017 
     32            const double theta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    3233            const double r = sin_val*sqrt(rA - rB*cos(theta)); 
    3334            const double be = sas_2J1x_x(q*r); 
    34             inner_sum += GAUSS_W[j] * be * be; 
     35            inner_sum += Gauss76Wt[j] * be * be; 
    3536        } 
    3637        //now calculate the value of the inner integral 
    … …  
    3940        //now calculate outer integral 
    4041        const double si = sas_sinx_x(q*0.5*length*cos_val); 
    41         outer_sum += GAUSS_W[i] * inner_sum * si * si; 
     42        outer_sum += Gauss76Wt[i] * inner_sum * si * si; 
    4243    } 
    4344    outer_sum *= 0.5*(vb-va); 
    … …  
    5455 
    5556static double 
    56 Iqabc(double qa, double qb, double qc, 
     57Iqxy(double qa, double qb, double qc, 
    5758     double radius_minor, double r_ratio, double length, 
    5859     double sld, double solvent_sld) 

  • sasmodels/models/elliptical_cylinder.py

    r2d81cfe r2d81cfe  
    121121# pylint: enable=bad-whitespace, line-too-long 
    122122 
    123 source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "elliptical_cylinder.c"] 
     123source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "lib/gauss20.c", 
     124          "elliptical_cylinder.c"] 
    124125 
    125126demo = dict(scale=1, background=0, radius_minor=100, axis_ratio=1.5, length=400.0, 

  • sasmodels/models/fcc_paracrystal.c

    r108e70e rf728001  
    5353 
    5454    double outer_sum = 0.0; 
    55     for(int i=0; i<GAUSS_N; i++) { 
     55    for(int i=0; i<150; i++) { 
    5656        double inner_sum = 0.0; 
    57         const double theta = GAUSS_Z[i]*theta_m + theta_b; 
     57        const double theta = Gauss150Z[i]*theta_m + theta_b; 
    5858        double sin_theta, cos_theta; 
    5959        SINCOS(theta, sin_theta, cos_theta); 
    6060        const double qc = q*cos_theta; 
    6161        const double qab = q*sin_theta; 
    62         for(int j=0;j<GAUSS_N;j++) { 
    63             const double phi = GAUSS_Z[j]*phi_m + phi_b; 
     62        for(int j=0;j<150;j++) { 
     63            const double phi = Gauss150Z[j]*phi_m + phi_b; 
    6464            double sin_phi, cos_phi; 
    6565            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    6767            const double qb = qab*sin_phi; 
    6868            const double form = fcc_Zq(qa, qb, qc, dnn, d_factor); 
    69             inner_sum += GAUSS_W[j] * form; 
     69            inner_sum += Gauss150Wt[j] * form; 
    7070        } 
    7171        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    72         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     72        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
    7373    } 
    7474    outer_sum *= theta_m; 
    … …  
    8080 
    8181 
    82 static double Iqabc(double qa, double qb, double qc, 
     82static double Iqxy(double qa, double qb, double qc, 
    8383    double dnn, double d_factor, double radius, 
    8484    double sld, double solvent_sld) 

  • sasmodels/models/flexible_cylinder_elliptical.c

    r74768cb r592343f  
    1717    double sum=0.0; 
    1818 
    19     for(int i=0;i<GAUSS_N;i++) { 
    20         const double zi = ( GAUSS_Z[i] + 1.0 )*M_PI_4; 
     19    for(int i=0;i<N_POINTS_76;i++) { 
     20        const double zi = ( Gauss76Z[i] + 1.0 )*M_PI_4; 
    2121        double sn, cn; 
    2222        SINCOS(zi, sn, cn); 
    2323        const double arg = q*sqrt(a*a*sn*sn + b*b*cn*cn); 
    2424        const double yyy = sas_2J1x_x(arg); 
    25         sum += GAUSS_W[i] * yyy * yyy; 
     25        sum += Gauss76Wt[i] * yyy * yyy; 
    2626    } 
    2727    sum *= 0.5; 

  • sasmodels/models/hollow_cylinder.c

    r108e70e rbecded3  
    3838 
    3939    double summ = 0.0;            //initialize intergral 
    40     for (int i=0;i<GAUSS_N;i++) { 
    41         const double cos_theta = 0.5*( GAUSS_Z[i] * (upper-lower) + lower + upper ); 
     40    for (int i=0;i<76;i++) { 
     41        const double cos_theta = 0.5*( Gauss76Z[i] * (upper-lower) + lower + upper ); 
    4242        const double sin_theta = sqrt(1.0 - cos_theta*cos_theta); 
    4343        const double form = _fq(q*sin_theta, q*cos_theta, 
    4444                                radius, thickness, length); 
    45         summ += GAUSS_W[i] * form * form; 
     45        summ += Gauss76Wt[i] * form * form; 
    4646    } 
    4747 
    … …  
    5252 
    5353static double 
    54 Iqac(double qab, double qc, 
     54Iqxy(double qab, double qc, 
    5555    double radius, double thickness, double length, 
    5656    double sld, double solvent_sld) 

  • sasmodels/models/hollow_rectangular_prism.c

    r108e70e r1e7b0db0  
    11double 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, 
     2double Iq(double q, double sld, double solvent_sld, double length_a,  
    33          double b2a_ratio, double c2a_ratio, double thickness); 
    44 
    … …  
    3737    const double v2a = 0.0; 
    3838    const double v2b = M_PI_2;  //phi integration limits 
     39     
     40    double outer_sum = 0.0; 
     41    for(int i=0; i<76; i++) { 
    3942 
    40     double outer_sum = 0.0; 
    41     for(int i=0; i<GAUSS_N; i++) { 
    42  
    43         const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
     43        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
    4444        double sin_theta, cos_theta; 
    4545        SINCOS(theta, sin_theta, cos_theta); 
    … …  
    4949 
    5050        double inner_sum = 0.0; 
    51         for(int j=0; j<GAUSS_N; j++) { 
     51        for(int j=0; j<76; j++) { 
    5252 
    53             const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
     53            const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
    5454            double sin_phi, cos_phi; 
    5555            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    6666            const double AP2 = vol_core * termA2 * termB2 * termC2; 
    6767 
    68             inner_sum += GAUSS_W[j] * square(AP1-AP2); 
     68            inner_sum += Gauss76Wt[j] * square(AP1-AP2); 
    6969        } 
    7070        inner_sum *= 0.5 * (v2b-v2a); 
    7171 
    72         outer_sum += GAUSS_W[i] * inner_sum * sin(theta); 
     72        outer_sum += Gauss76Wt[i] * inner_sum * sin(theta); 
    7373    } 
    7474    outer_sum *= 0.5*(v1b-v1a); 
    … …  
    8484    return 1.0e-4 * delrho * delrho * form; 
    8585} 
    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  
    55This model provides the form factor, $P(q)$, for a hollow rectangular 
    66parallelepiped with a wall of thickness $\Delta$. 
    7  
     7It computes only the 1D scattering, not the 2D. 
    88 
    99Definition 
    … …  
    6666(which is unitless). 
    6767 
    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. 
     68**The 2D scattering intensity is not computed by this model.** 
    8969 
    9070 
    … …  
    133113              ["thickness", "Ang", 1, [0, inf], "volume", 
    134114               "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"], 
    141115             ] 
    142116 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    r74768cb rab2aea8  
    11double form_volume(double length_a, double b2a_ratio, double c2a_ratio); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a, 
     2double Iq(double q, double sld, double solvent_sld, double length_a,  
    33          double b2a_ratio, double c2a_ratio); 
    44 
    … …  
    2929    const double v2a = 0.0; 
    3030    const double v2b = M_PI_2;  //phi integration limits 
    31  
     31     
    3232    double outer_sum = 0.0; 
    33     for(int i=0; i<GAUSS_N; i++) { 
    34         const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
     33    for(int i=0; i<76; i++) { 
     34        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
    3535 
    3636        double sin_theta, cos_theta; 
    … …  
    4444 
    4545        double inner_sum = 0.0; 
    46         for(int j=0; j<GAUSS_N; j++) { 
    47             const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
     46        for(int j=0; j<76; j++) { 
     47            const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
    4848 
    4949            double sin_phi, cos_phi; 
    … …  
    6262                * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi ); 
    6363 
    64             inner_sum += GAUSS_W[j] * square(AL+AT); 
     64            inner_sum += Gauss76Wt[j] * square(AL+AT); 
    6565        } 
    6666 
    6767        inner_sum *= 0.5 * (v2b-v2a); 
    68         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     68        outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 
    6969    } 
    7070 

  • sasmodels/models/lib/gauss150.c

    r74768cb r994d77f  
    77 * 
    88 */ 
    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 
    199 
    2010// Gaussians 
    21 constant double Gauss150Z[152]={ 
     11constant double Gauss150Z[150]={ 
    2212        -0.9998723404457334, 
    2313        -0.9993274305065947, 
    … …  
    169159        0.9983473449340834, 
    170160        0.9993274305065947, 
    171         0.9998723404457334, 
    172         0., 
    173         0. 
     161        0.9998723404457334 
    174162}; 
    175163 
    176 constant double Gauss150Wt[152]={ 
     164constant double Gauss150Wt[150]={ 
    177165        0.0003276086705538, 
    178166        0.0007624720924706, 
    … …  
    324312        0.0011976474864367, 
    325313        0.0007624720924706, 
    326         0.0003276086705538, 
    327         0., 
    328         0. 
     314        0.0003276086705538 
    329315}; 

  • sasmodels/models/lib/gauss20.c

    r74768cb r994d77f  
    77 * 
    88 */ 
    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 
    179 
    1810// Gaussians 

  • sasmodels/models/lib/gauss76.c

    r74768cb r66d119f  
    77 * 
    88 */ 
    9  #ifdef GAUSS_N 
    10  # undef GAUSS_N 
    11  # undef GAUSS_Z 
    12  # undef GAUSS_W 
    13  #endif 
    14  #define GAUSS_N 76 
    15  #define GAUSS_Z Gauss76Z 
    16  #define GAUSS_W Gauss76Wt 
     9#define N_POINTS_76 76 
    1710 
    1811// Gaussians 
    19 constant double Gauss76Wt[76]={ 
     12constant double Gauss76Wt[N_POINTS_76]={ 
    2013        .00126779163408536,             //0 
    2114        .00294910295364247, 
    … …  
    9689}; 
    9790 
    98 constant double Gauss76Z[76]={ 
     91constant double Gauss76Z[N_POINTS_76]={ 
    9992        -.999505948362153,              //0 
    10093        -.997397786355355, 

  • sasmodels/models/line.py

    r108e70e r2d81cfe  
    5757Iq.vectorized = True # Iq accepts an array of q values 
    5858 
    59  
    6059def Iqxy(qx, qy, *args): 
    6160    """ 
    … …  
    7069 
    7170Iqxy.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 #""" 
    8071 
    8172def random(): 

  • sasmodels/models/parallelepiped.c

    r108e70e r9b7b23f  
    2323    double outer_total = 0; //initialize integral 
    2424 
    25     for( int i=0; i<GAUSS_N; i++) { 
    26         const double sigma = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
     25    for( int i=0; i<76; i++) { 
     26        const double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
    2727        const double mu_proj = mu * sqrt(1.0-sigma*sigma); 
    2828 
    … …  
    3030        // corresponding to angles from 0 to pi/2. 
    3131        double inner_total = 0.0; 
    32         for(int j=0; j<GAUSS_N; j++) { 
    33             const double uu = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     32        for(int j=0; j<76; j++) { 
     33            const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    3434            double sin_uu, cos_uu; 
    3535            SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
    3636            const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 
    3737            const double si2 = sas_sinx_x(mu_proj * cos_uu); 
    38             inner_total += GAUSS_W[j] * square(si1 * si2); 
     38            inner_total += Gauss76Wt[j] * square(si1 * si2); 
    3939        } 
    4040        inner_total *= 0.5; 
    4141 
    4242        const double si = sas_sinx_x(mu * c_scaled * sigma); 
    43         outer_total += GAUSS_W[i] * inner_total * si * si; 
     43        outer_total += Gauss76Wt[i] * inner_total * si * si; 
    4444    } 
    4545    outer_total *= 0.5; 
    … …  
    5353 
    5454static double 
    55 Iqabc(double qa, double qb, double qc, 
     55Iqxy(double qa, double qb, double qc, 
    5656    double sld, 
    5757    double solvent_sld, 

  • sasmodels/models/pringle.c

    r74768cb r1e7b0db0  
    2929    double sumC = 0.0;          // initialize integral 
    3030    double r; 
    31     for (int i=0; i < GAUSS_N; i++) { 
    32         r = GAUSS_Z[i]*zm + zb; 
     31    for (int i=0; i < 76; i++) { 
     32        r = Gauss76Z[i]*zm + zb; 
    3333 
    3434        const double qrs = r*q_sin_psi; 
    3535        const double qrrc = r*r*q_cos_psi; 
    3636 
    37         double y = GAUSS_W[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs); 
     37        double y = Gauss76Wt[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs); 
    3838        double S, C; 
    3939        SINCOS(alpha*qrrc, S, C); 
    … …  
    8686 
    8787    double sum = 0.0; 
    88     for (int i = 0; i < GAUSS_N; i++) { 
    89         double psi = GAUSS_Z[i]*zm + zb; 
     88    for (int i = 0; i < 76; i++) { 
     89        double psi = Gauss76Z[i]*zm + zb; 
    9090        double sin_psi, cos_psi; 
    9191        SINCOS(psi, sin_psi, cos_psi); 
    … …  
    9393        double sinc_term = square(sas_sinx_x(q * thickness * cos_psi / 2.0)); 
    9494        double pringle_kernel = 4.0 * sin_psi * bessel_term * sinc_term; 
    95         sum += GAUSS_W[i] * pringle_kernel; 
     95        sum += Gauss76Wt[i] * pringle_kernel; 
    9696    } 
    9797 

  • sasmodels/models/rectangular_prism.c

    r108e70e r1e7b0db0  
    11double form_volume(double length_a, double b2a_ratio, double c2a_ratio); 
    2 double Iq(double q, double sld, double solvent_sld, double length_a, 
     2double Iq(double q, double sld, double solvent_sld, double length_a,  
    33          double b2a_ratio, double c2a_ratio); 
    44 
    … …  
    2626    const double v2a = 0.0; 
    2727    const double v2b = M_PI_2;  //phi integration limits 
    28  
     28     
    2929    double outer_sum = 0.0; 
    30     for(int i=0; i<GAUSS_N; i++) { 
    31         const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
     30    for(int i=0; i<76; i++) { 
     31        const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
    3232        double sin_theta, cos_theta; 
    3333        SINCOS(theta, sin_theta, cos_theta); 
    … …  
    3636 
    3737        double inner_sum = 0.0; 
    38         for(int j=0; j<GAUSS_N; j++) { 
    39             double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
     38        for(int j=0; j<76; j++) { 
     39            double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
    4040            double sin_phi, cos_phi; 
    4141            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    4545            const double termB = sas_sinx_x(q * b_half * sin_theta * cos_phi); 
    4646            const double AP = termA * termB * termC; 
    47             inner_sum += GAUSS_W[j] * AP * AP; 
     47            inner_sum += Gauss76Wt[j] * AP * AP; 
    4848        } 
    4949        inner_sum = 0.5 * (v2b-v2a) * inner_sum; 
    50         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     50        outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 
    5151    } 
    5252 
    5353    double answer = 0.5*(v1b-v1a)*outer_sum; 
    5454 
    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  
    5757    // 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  
    5959    // 0 and pi/2, instead of 0 to pi. 
    6060    answer /= M_PI_2; //Form factor P(q) 
    … …  
    6464    answer *= square((sld-solvent_sld)*volume); 
    6565 
    66     // Convert from [1e-12 A-1] to [cm-1] 
     66    // Convert from [1e-12 A-1] to [cm-1]  
    6767    answer *= 1.0e-4; 
    6868 
    6969    return answer; 
    7070} 
    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  
    1212the prism (e.g. setting $b/a = 1$ and $c/a = 1$ and applying polydispersity 
    1313to *a* will generate a distribution of cubes of different sizes). 
     14Note also that, contrary to :ref:`parallelepiped`, it does not compute 
     15the 2D scattering. 
    1416 
    1517 
    … …  
    2426that reference), with $\theta$ corresponding to $\alpha$ in that paper, 
    2527and not to the usual convention used for example in the 
    26 :ref:`parallelepiped` model. 
     28:ref:`parallelepiped` model. As the present model does not compute 
     29the 2D scattering, this has no further consequences. 
    2730 
    2831In this model the scattering from a massive parallelepiped with an 
    … …  
    6265units) *scale* represents the volume fraction (which is unitless). 
    6366 
    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  
     67**The 2D scattering intensity is not computed by this model.** 
    8668 
    8769 
    … …  
    126108              ["c2a_ratio", "", 1, [0, inf], "volume", 
    127109               "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"], 
    134110             ] 
    135111 

  • sasmodels/models/sc_paracrystal.c

    r108e70e rf728001  
    5454 
    5555    double outer_sum = 0.0; 
    56     for(int i=0; i<GAUSS_N; i++) { 
     56    for(int i=0; i<150; i++) { 
    5757        double inner_sum = 0.0; 
    58         const double theta = GAUSS_Z[i]*theta_m + theta_b; 
     58        const double theta = Gauss150Z[i]*theta_m + theta_b; 
    5959        double sin_theta, cos_theta; 
    6060        SINCOS(theta, sin_theta, cos_theta); 
    6161        const double qc = q*cos_theta; 
    6262        const double qab = q*sin_theta; 
    63         for(int j=0;j<GAUSS_N;j++) { 
    64             const double phi = GAUSS_Z[j]*phi_m + phi_b; 
     63        for(int j=0;j<150;j++) { 
     64            const double phi = Gauss150Z[j]*phi_m + phi_b; 
    6565            double sin_phi, cos_phi; 
    6666            SINCOS(phi, sin_phi, cos_phi); 
    … …  
    6868            const double qb = qab*sin_phi; 
    6969            const double form = sc_Zq(qa, qb, qc, dnn, d_factor); 
    70             inner_sum += GAUSS_W[j] * form; 
     70            inner_sum += Gauss150Wt[j] * form; 
    7171        } 
    7272        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    73         outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
     73        outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
    7474    } 
    7575    outer_sum *= theta_m; 
    … …  
    8282 
    8383static double 
    84 Iqabc(double qa, double qb, double qc, 
     84Iqxy(double qa, double qb, double qc, 
    8585    double dnn, double d_factor, double radius, 
    8686    double sld, double solvent_sld) 

  • sasmodels/models/stacked_disks.c

    r108e70e rbecded3  
    8181    double halfheight = 0.5*thick_core; 
    8282 
    83     for(int i=0; i<GAUSS_N; i++) { 
    84         double zi = (GAUSS_Z[i] + 1.0)*M_PI_4; 
     83    for(int i=0; i<N_POINTS_76; i++) { 
     84        double zi = (Gauss76Z[i] + 1.0)*M_PI_4; 
    8585        double sin_alpha, cos_alpha; // slots to hold sincos function output 
    8686        SINCOS(zi, sin_alpha, cos_alpha); 
    … …  
    9595                           solvent_sld, 
    9696                           d); 
    97         summ += GAUSS_W[i] * yyy * sin_alpha; 
     97        summ += Gauss76Wt[i] * yyy * sin_alpha; 
    9898    } 
    9999 
    … …  
    142142 
    143143static double 
    144 Iqac(double qab, double qc, 
     144Iqxy(double qab, double qc, 
    145145    double thick_core, 
    146146    double thick_layer, 

  • sasmodels/models/triaxial_ellipsoid.c

    r108e70e rbecded3  
    2121    const double zb = M_PI_4; 
    2222    double outer = 0.0; 
    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; 
     23    for (int i=0;i<76;i++) { 
     24        //const double u = Gauss76Z[i]*(upper-lower)/2 + (upper + lower)/2; 
     25        const double phi = Gauss76Z[i]*zm + zb; 
    2626        const double pa_sinsq_phi = pa*square(sin(phi)); 
    2727 
    … …  
    2929        const double um = 0.5; 
    3030        const double ub = 0.5; 
    31         for (int j=0;j<GAUSS_N;j++) { 
     31        for (int j=0;j<76;j++) { 
    3232            // translate a point in [-1,1] to a point in [0, 1] 
    33             const double usq = square(GAUSS_Z[j]*um + ub); 
     33            const double usq = square(Gauss76Z[j]*um + ub); 
    3434            const double r = radius_equat_major*sqrt(pa_sinsq_phi*(1.0-usq) + 1.0 + pc*usq); 
    3535            const double fq = sas_3j1x_x(q*r); 
    36             inner += GAUSS_W[j] * fq * fq; 
     36            inner += Gauss76Wt[j] * fq * fq; 
    3737        } 
    38         outer += GAUSS_W[i] * inner;  // correcting for dx later 
     38        outer += Gauss76Wt[i] * inner;  // correcting for dx later 
    3939    } 
    4040    // translate integration ranges from [-1,1] to [lower,upper] and normalize by 4 pi 
    … …  
    4646 
    4747static double 
    48 Iqabc(double qa, double qb, double qc, 
     48Iqxy(double qa, double qb, double qc, 
    4949    double sld, 
    5050    double sld_solvent, 
Note: See TracChangeset for help on using the changeset viewer.