Changeset 1941ec6 in sasmodels


Ignore:
Timestamp:
Dec 18, 2017 3:52:27 PM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Children:
67cc0ff
Parents:
8224d24 (diff), 93fbc34 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
git-author:
Paul Kienzle <pkienzle@…> (12/18/17 15:51:43)
git-committer:
Paul Kienzle <pkienzle@…> (12/18/17 15:52:27)
Message:

Merge remote-tracking branch 'upstream/ticket-1043'

Files:
10 added
44 edited

Legend:

Unmodified
Added
Removed

  • .gitignore

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

  • doc/guide/plugin.rst

    r0a9fcab rc654160  
    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 Iqxy() and form_volume(). And** *scale* **and** *background* **parameters are 
    295 implicit to all models, so they do not need to be included in the parameter table.** 
     294Iqac(), Iqabc() and form_volume(). And** *scale* **and** *background* 
     295**parameters are implicit to all models, so they do not need to be included 
     296in the parameter table.** 
    296297 
    297298- **"name"** is the name of the parameter shown on the FitPage. 
    … …  
    362363    scattered intensity. 
    363364 
    364   - "volume" parameters are passed to Iq(), Iqxy(), and form_volume(), and 
    365     have polydispersity loops generated automatically. 
    366  
    367   - "orientation" parameters are only passed to Iqxy(), and have angular 
    368     dispersion. 
     365  - "volume" parameters are passed to Iq(), Iqac(), Iqabc() and form_volume(), 
     366    and have polydispersity loops generated automatically. 
     367 
     368  - "orientation" parameters are not passed, but instead are combined with 
     369    orientation dispersity to translate *qx* and *qy* to *qa*, *qb* and *qc*. 
     370    These parameters should appear at the end of the table with the specific 
     371    names *theta*, *phi* and for asymmetric shapes *psi*, in that order. 
    369372 
    370373Some models will have integer parameters, such as number of pearls in the 
    … …  
    419422That is, the individual models do not need to include polydispersity 
    420423calculations, but instead rely on numerical integration to compute the 
    421 appropriately smeared pattern.   Angular dispersion values over polar angle 
    422 $\theta$ requires an additional $\cos \theta$ weighting due to decreased 
    423 arc length for the equatorial angle $\phi$ with increasing latitude. 
     424appropriately smeared pattern. 
    424425 
    425426Python Models 
    … …  
    468469barbell).  If I(q; pars) is NaN for any $q$, then those parameters will be 
    469470ignored, and not included in the calculation of the weighted polydispersity. 
    470  
    471 Similar to *Iq*, you can define *Iqxy(qx, qy, par1, par2, ...)* where the 
    472 parameter list includes any orientation parameters.  If *Iqxy* is not defined, 
    473 then it will default to *Iqxy = Iq(sqrt(qx**2+qy**2), par1, par2, ...)*. 
    474471 
    475472Models should define *form_volume(par1, par2, ...)* where the parameter 
    … …  
    497494    } 
    498495 
    499 *Iqxy* is similar to *Iq*, except it uses parameters *qx, qy* instead of *q*, 
    500 and it includes orientation parameters. 
    501  
    502496*form_volume* defines the volume of the shape. As in python models, it 
    503497includes only the volume parameters. 
    504498 
    505 *Iqxy* will default to *Iq(sqrt(qx**2 + qy**2), par1, ...)* and 
    506 *form_volume* will default to 1.0. 
    507  
    508499**source=['fn.c', ...]** includes the listed C source files in the 
    509 program before *Iq* and *Iqxy* are defined. This allows you to extend the 
    510 library of C functions available to your model. 
     500program before *Iq* and *form_volume* are defined. This allows you to 
     501extend the library of C functions available to your model. 
     502 
     503*c_code* includes arbitrary C code into your kernel, which can be 
     504handy for defining helper functions for *Iq* and *form_volume*. Note that 
     505you 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 
     507external C file and add that file to the list of sources. 
    511508 
    512509Models are defined using double precision declarations for the 
    … …  
    532529 
    533530    #define INVALID(v) (v.bell_radius < v.radius) 
     531 
     532The INVALID define can go into *Iq*, or *c_code*, or an external C file 
     533listed in *source*. 
     534 
     535Oriented Shapes 
     536............... 
     537 
     538If the scattering is dependent on the orientation of the shape, then you 
     539will need to include *orientation* parameters *theta*, *phi* and *psi* 
     540at the end of the parameter table.  Shape orientation uses *a*, *b* and *c* 
     541axes, corresponding to the *x*, *y* and *z* axes in the laboratory coordinate 
     542system, with *z* along the beam and *x*-*y* in the detector plane, with *x* 
     543horizontal and *y* vertical.  The *psi* parameter rotates the shape 
     544about its *c* axis, the *theta* parameter then rotates the *c* axis toward 
     545the *x* axis of the detector, then *phi* rotates the shape in the detector 
     546plane.  (Prior to these rotations, orientation dispersity will be applied 
     547as roll-pitch-yaw, rotating *c*, then *b* then *a* in the shape coordinate 
     548system.)  A particular *qx*, *qy* point on the detector, then corresponds 
     549to *qa*, *qb*, *qc* with respect to the shape. 
     550 
     551The 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 
     553symmetric about *c* then *psi* is not needed, and the model is called 
     554as *Iqac(qab, qc, par1, par2, ...)*.  In either case, the orientation 
     555parameters are not included in the function call. 
     556 
     557For 1D oriented shapes, an integral over all angles is usually needed for 
     558the *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 
     570for 
     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 
     578Using the $z, w$ values for Gauss-Legendre integration in "lib/gauss76.c", the 
     579numerical 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 
     600The *z* values for the Gauss-Legendre integration extends from -1 to 1, so 
     601the double sum of *w[i]w[j]* explains the factor of 4.  Correcting for the 
     602average *dz[i]dz[j]* gives $(1-0) \cdot (\pi/2-0) = \pi/2$.  The $8/(4 \pi)$ 
     603factor comes from the integral over the quadrant.  With less symmetry (eg., 
     604in the bcc and fcc paracrystal models), then an integral over the entire 
     605sphere may be necessary. 
     606 
     607For simpler models which are rotationally symmetric a single integral 
     608suffices: 
     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 
     616for 
     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 
     624with 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 
     637Magnetism 
     638......... 
     639 
     640Magnetism is supported automatically for all shapes by modifying the 
     641effective SLD of particle according to the Halpern-Johnson vector 
     642describing the interaction between neutron spin and magnetic field.  All 
     643parameters marked as type *sld* in the parameter table are treated as 
     644possibly magnetic particles with magnitude *M0* and direction 
     645*mtheta* and *mphi*.  Polarization parameters are also provided 
     646automatically for magnetic models to set the spin state of the measurement. 
     647 
     648For more complicated systems where magnetism is not uniform throughout 
     649the individual particles, you will need to write your own models. 
     650You should not mark the nuclear sld as type *sld*, but instead leave 
     651them unmarked and provide your own magnetism and polarization parameters. 
     652For 2D measurements you will need $(q_x, q_y)$ values for the measurement 
     653to compute the proper magnetism and orientation, which you can implement 
     654using *Iqxy(qx, qy, par1, par2, ...)*. 
    534655 
    535656Special Functions 
    … …  
    796917show a 50x improvement or more over the equivalent pure python model. 
    797918 
    798 External C Models 
    799 ................. 
    800  
    801 External C models are very much like embedded C models, except that 
    802 *Iq*, *Iqxy* and *form_volume* are defined in an external source file 
    803 loaded using the *source=[...]* statement. You need to supply the function 
    804 declarations for each of these that you need instead of building them 
    805 automatically from the parameter table. 
    806  
    807919 
    808920.. _Form_Factors: 
    … …  
    10061118variable name *Rg* for example because $R_g$ is the right name for the model 
    10071119parameter then ignore the lint errors.  Also, ignore *missing-docstring* 
    1008 for standard model functions *Iq*, *Iqxy*, etc. 
     1120for standard model functions *Iq*, *Iqac*, etc. 
    10091121 
    10101122We will have delinting sessions at the SasView Code Camps, where we can 

  • explore/asymint.py

    r1820208 ra1c32c2  
    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(0.5*a*qa/2) 
    89         siB = env.sas_sinx_x(0.5*b*qb/2) 
    90         siC = env.sas_sinx_x(0.5*c*qc/2) 
     88        siA = env.sas_sinx_x(a*qa/2) 
     89        siB = env.sas_sinx_x(b*qb/2) 
     90        siC = env.sas_sinx_x(c*qc/2) 
    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 
     97def 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) 
    95128    return norm, Fq 
    96129 
    … …  
    184217    NORM, KERNEL = make_parallelepiped(A, B, C) 
    185218    NORM_MP, KERNEL_MP = make_parallelepiped(A, B, C, env=MPenv) 
     219elif 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) 
    186237elif shape == 'paracrystal': 
    187238    LATTICE = 'bcc' 
    … …  
    342393    print("gauss-150", *gauss_quad_2d(Q, n=150)) 
    343394    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)) 
    344397    #gridded_2d(Q, n=2**8+1) 
    345398    gridded_2d(Q, n=2**10+1) 
    346     #gridded_2d(Q, n=2**13+1) 
     399    #gridded_2d(Q, n=2**12+1) 
    347400    #gridded_2d(Q, n=2**15+1) 
    348     if shape != 'paracrystal':  # adaptive forms are too slow! 
     401    if shape not in ('paracrystal', 'core_shell_parallelepiped'): 
     402        # adaptive forms on models for which the calculations are fast enough 
    349403        print("dblquad", *scipy_dblquad_2d(Q)) 
    350404        print("semi-romberg-100", *semi_romberg_2d(Q, n=100)) 

  • sasmodels/compare.py

    r2d81cfe ra261a83  
    4242from .data import plot_theory, empty_data1D, empty_data2D, load_data 
    4343from .direct_model import DirectModel, get_mesh 
    44 from .generate import FLOAT_RE 
     44from .generate import FLOAT_RE, set_integration_size 
    4545from .weights import plot_weights 
    4646 
    … …  
    693693        data = empty_data2D(q, resolution=res) 
    694694        data.accuracy = opts['accuracy'] 
    695         set_beam_stop(data, 0.0004) 
     695        set_beam_stop(data, qmin) 
    696696        index = ~data.mask 
    697697    else: 
    … …  
    706706    return data, index 
    707707 
    708 def make_engine(model_info, data, dtype, cutoff): 
     708def make_engine(model_info, data, dtype, cutoff, ngauss=0): 
    709709    # type: (ModelInfo, Data, str, float) -> Calculator 
    710710    """ 
    … …  
    714714    than OpenCL. 
    715715    """ 
     716    if ngauss: 
     717        set_integration_size(model_info, ngauss) 
     718 
    716719    if dtype is None or not dtype.endswith('!'): 
    717720        return eval_opencl(model_info, data, dtype=dtype, cutoff=cutoff) 
    … …  
    954957    'poly', 'mono', 'cutoff=', 
    955958    'magnetic', 'nonmagnetic', 
    956     'accuracy=', 
     959    'accuracy=', 'ngauss=', 
    957960    'neval=',  # for timing... 
    958961 
    … …  
    10891092        'show_weights' : False, 
    10901093        'sphere'    : 0, 
     1094        'ngauss'    : '0', 
    10911095    } 
    10921096    for arg in flags: 
    … …  
    11151119        elif arg.startswith('-engine='):   opts['engine'] = arg[8:] 
    11161120        elif arg.startswith('-neval='):    opts['count'] = arg[7:] 
     1121        elif arg.startswith('-ngauss='):   opts['ngauss'] = arg[8:] 
    11171122        elif arg.startswith('-random='): 
    11181123            opts['seed'] = int(arg[8:]) 
    … …  
    11691174 
    11701175    comparison = any(PAR_SPLIT in v for v in values) 
     1176 
    11711177    if PAR_SPLIT in name: 
    11721178        names = name.split(PAR_SPLIT, 2) 
    … …  
    11811187        return None 
    11821188 
     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 
    11831195    if PAR_SPLIT in opts['engine']: 
    11841196        opts['engine'] = opts['engine'].split(PAR_SPLIT, 2) 
    … …  
    11991211        opts['cutoff'] = [float(opts['cutoff'])]*2 
    12001212 
    1201     base = make_engine(model_info[0], data, opts['engine'][0], opts['cutoff'][0]) 
     1213    base = make_engine(model_info[0], data, opts['engine'][0], 
     1214                       opts['cutoff'][0], opts['ngauss'][0]) 
    12021215    if comparison: 
    1203         comp = make_engine(model_info[1], data, opts['engine'][1], opts['cutoff'][1]) 
     1216        comp = make_engine(model_info[1], data, opts['engine'][1], 
     1217                           opts['cutoff'][1], opts['ngauss'][1]) 
    12041218    else: 
    12051219        comp = None 
    … …  
    12741288        if model_info != model_info2: 
    12751289            pars2 = randomize_pars(model_info2, pars2) 
    1276             limit_dimensions(model_info, pars2, maxdim) 
     1290            limit_dimensions(model_info2, pars2, maxdim) 
    12771291            # Share values for parameters with the same name 
    12781292            for k, v in pars.items(): 

  • sasmodels/details.py

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

  • sasmodels/generate.py

    rdb03406 r108e70e  
    77    particular dimensions averaged over all orientations. 
    88 
    9     *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy for a form 
    10     with particular dimensions for a single orientation. 
    11  
    12     *Imagnetic(qx, qy, result[], p1, p2, ...)* returns the scattering for the 
    13     polarized neutron spin states (up-up, up-down, down-up, down-down) for 
    14     a form with particular dimensions for a single orientation. 
     9    *Iqac(qab, qc, p1, p2, ...)* returns the scattering at qab, qc 
     10    for a rotationally symmetric form with particular dimensions. 
     11    qab, qc are determined from shape orientation and scattering angles. 
     12    This call is used if the shape has orientation parameters theta and phi. 
     13 
     14    *Iqabc(qa, qb, qc, p1, p2, ...)* returns the scattering at qa, qb, qc 
     15    for a form with particular dimensions.  qa, qb, qc are determined from 
     16    shape orientation and scattering angles. This call is used if the shape 
     17    has orientation parameters theta, phi and psi. 
     18 
     19    *Iqxy(qx, qy, p1, p2, ...)* returns the scattering at qx, qy.  Use this 
     20    to create an arbitrary 2D theory function, needed for q-dependent 
     21    background functions and for models with non-uniform magnetism. 
    1522 
    1623    *form_volume(p1, p2, ...)* returns the volume of the form with particular 
    … …  
    3138scale and background parameters for each model. 
    3239 
    33 *Iq*, *Iqxy*, *Imagnetic* and *form_volume* should be stylized C-99 
    34 functions written for OpenCL.  All functions need prototype declarations 
    35 even if the are defined before they are used.  OpenCL does not support 
    36 *#include* preprocessor directives, so instead the list of includes needs 
    37 to be given as part of the metadata in the kernel module definition. 
    38 The included files should be listed using a path relative to the kernel 
    39 module, or if using "lib/file.c" if it is one of the standard includes 
    40 provided with the sasmodels source.  The includes need to be listed in 
    41 order so that functions are defined before they are used. 
     40C code should be stylized C-99 functions written for OpenCL. All functions 
     41need prototype declarations even if the are defined before they are used. 
     42Although OpenCL supports *#include* preprocessor directives, the list of 
     43includes should be given as part of the metadata in the kernel module 
     44definition. The included files should be listed using a path relative to the 
     45kernel module, or if using "lib/file.c" if it is one of the standard includes 
     46provided with the sasmodels source. The includes need to be listed in order 
     47so that functions are defined before they are used. 
    4248 
    4349Floating point values should be declared as *double*.  For single precision 
    … …  
    107113    present, the volume ratio is 1. 
    108114 
    109     *form_volume*, *Iq*, *Iqxy*, *Imagnetic* are strings containing the 
    110     C source code for the body of the volume, Iq, and Iqxy functions 
     115    *form_volume*, *Iq*, *Iqac*, *Iqabc* are strings containing 
     116    the C source code for the body of the volume, Iq, and Iqac functions 
    111117    respectively.  These can also be defined in the last source file. 
    112118 
    113     *Iq* and *Iqxy* also be instead be python functions defining the 
     119    *Iq*, *Iqac*, *Iqabc* also be instead be python functions defining the 
    114120    kernel.  If they are marked as *Iq.vectorized = True* then the 
    115121    kernel is passed the entire *q* vector at once, otherwise it is 
    … …  
    168174from zlib import crc32 
    169175from inspect import currentframe, getframeinfo 
     176import logging 
    170177 
    171178import numpy as np  # type: ignore 
    … …  
    181188    pass 
    182189# pylint: enable=unused-import 
     190 
     191logger = logging.getLogger(__name__) 
    183192 
    184193# jitter projection to use in the kernel code.  See explore/jitter.py 
    … …  
    270279""" 
    271280 
     281 
     282def 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] 
    272299 
    273300def format_units(units): 
    … …  
    608635 
    609636""" 
    610 def _gen_fn(name, pars, body, filename, line): 
    611     # type: (str, List[Parameter], str, str, int) -> str 
     637def _gen_fn(model_info, name, pars): 
     638    # type: (ModelInfo, str, List[Parameter]) -> str 
    612639    """ 
    613640    Generate a function given pars and body. 
    … …  
    621648    """ 
    622649    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) 
    623654    return _FN_TEMPLATE % { 
    624655        'name': name, 'pars': par_decl, 'body': body, 
    625         'filename': filename.replace('\\', '\\\\'), 'line': line, 
     656        'filename': filename.replace('\\', '\\\\'), 'line': lineno, 
    626657    } 
    627658 
    … …  
    638669 
    639670# type in IQXY pattern could be single, float, double, long double, ... 
    640 _IQXY_PATTERN = re.compile("^((inline|static) )? *([a-z ]+ )? *Iqxy *([(]|$)", 
     671_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]", 
    641672                           flags=re.MULTILINE) 
    642 def _have_Iqxy(sources): 
     673def find_xy_mode(source): 
    643674    # type: (List[str]) -> bool 
    644675    """ 
    645     Return true if any file defines Iqxy. 
     676    Return the xy mode as qa, qac, qabc or qxy. 
    646677 
    647678    Note this is not a C parser, and so can be easily confused by 
    648679    non-standard syntax.  Also, it will incorrectly identify the following 
    649     as having Iqxy:: 
     680    as having 2D models:: 
    650681 
    651682        /* 
    652         double Iqxy(qx, qy, ...) { ... fill this in later ... } 
     683        double Iqac(qab, qc, ...) { ... fill this in later ... } 
    653684        */ 
    654685 
    655     If you want to comment out an Iqxy function, use // on the front of the 
    656     line instead. 
    657     """ 
    658     for _path, code in sources: 
    659         if _IQXY_PATTERN.search(code): 
    660             return True 
    661     return False 
    662  
    663  
    664 def _add_source(source, code, path): 
     686    If you want to comment out the function, use // on the front of the 
     687    line:: 
     688 
     689        /* 
     690        // double Iqac(qab, qc, ...) { ... fill this in later ... } 
     691        */ 
     692 
     693    """ 
     694    for code in source: 
     695        m = _IQXY_PATTERN.search(code) 
     696        if m is not None: 
     697            return m.group('mode') 
     698    return 'qa' 
     699 
     700 
     701def _add_source(source, code, path, lineno=1): 
    665702    """ 
    666703    Add a file to the list of source code chunks, tagged with path and line. 
    667704    """ 
    668705    path = path.replace('\\', '\\\\') 
    669     source.append('#line 1 "%s"' % path) 
     706    source.append('#line %d "%s"' % (lineno, path)) 
    670707    source.append(code) 
    671708 
    … …  
    698735    user_code = [(f, open(f).read()) for f in model_sources(model_info)] 
    699736 
    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  
    705737    # Build initial sources 
    706738    source = [] 
    … …  
    710742 
    711743    if model_info.c_code: 
    712         source.append(model_info.c_code) 
     744        _add_source(source, model_info.c_code, model_info.filename, 
     745                    lineno=model_info.lineno.get('c_code', 1)) 
    713746 
    714747    # Make parameters for q, qx, qy so that we can use them in declarations 
    715     q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 
     748    q, qx, qy, qab, qa, qb, qc \ 
     749        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()] 
    716750    # Generate form_volume function, etc. from body only 
    717751    if isinstance(model_info.form_volume, str): 
    718752        pars = partable.form_volume_parameters 
    719         source.append(_gen_fn('form_volume', pars, model_info.form_volume, 
    720                               model_info.filename, model_info._form_volume_line)) 
     753        source.append(_gen_fn(model_info, 'form_volume', pars)) 
    721754    if isinstance(model_info.Iq, str): 
    722755        pars = [q] + partable.iq_parameters 
    723         source.append(_gen_fn('Iq', pars, model_info.Iq, 
    724                               model_info.filename, model_info._Iq_line)) 
     756        source.append(_gen_fn(model_info, 'Iq', pars)) 
    725757    if isinstance(model_info.Iqxy, str): 
    726         pars = [qx, qy] + partable.iqxy_parameters 
    727         source.append(_gen_fn('Iqxy', pars, model_info.Iqxy, 
    728                               model_info.filename, model_info._Iqxy_line)) 
     758        pars = [qx, qy] + partable.iq_parameters + partable.orientation_parameters 
     759        source.append(_gen_fn(model_info, 'Iqxy', pars)) 
     760    if isinstance(model_info.Iqac, str): 
     761        pars = [qab, qc] + partable.iq_parameters 
     762        source.append(_gen_fn(model_info, 'Iqac', pars)) 
     763    if isinstance(model_info.Iqabc, str): 
     764        pars = [qa, qb, qc] + partable.iq_parameters 
     765        source.append(_gen_fn(model_info, 'Iqabc', pars)) 
     766 
     767    # What kind of 2D model do we need?  Is it consistent with the parameters? 
     768    xy_mode = find_xy_mode(source) 
     769    if xy_mode == 'qabc' and not partable.is_asymmetric: 
     770        raise ValueError("asymmetric oriented models need to define Iqabc") 
     771    elif xy_mode == 'qac' and partable.is_asymmetric: 
     772        raise ValueError("symmetric oriented models need to define Iqac") 
     773    elif not partable.orientation_parameters and xy_mode in ('qac', 'qabc'): 
     774        raise ValueError("Unexpected function I%s for unoriented shape"%xy_mode) 
     775    elif partable.orientation_parameters and xy_mode not in ('qac', 'qabc'): 
     776        if xy_mode == 'qxy': 
     777            logger.warn("oriented shapes should define Iqac or Iqabc") 
     778        else: 
     779            raise ValueError("Expected function Iqac or Iqabc for oriented shape") 
    729780 
    730781    # Define the parameter table 
    … …  
    752803    if xy_mode == 'qabc': 
    753804        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
    754         call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqxy(%s)" % pars 
     805        call_iqxy = "#define CALL_IQ_ABC(_qa,_qb,_qc,_v) Iqabc(%s)" % pars 
    755806        clear_iqxy = "#undef CALL_IQ_ABC" 
    756807    elif xy_mode == 'qac': 
    757808        pars = ",".join(["_qa", "_qc"] + model_refs) 
    758         call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqxy(%s)" % pars 
     809        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars 
    759810        clear_iqxy = "#undef CALL_IQ_AC" 
    760     else:  # xy_mode == 'qa' 
     811    elif xy_mode == 'qa': 
    761812        pars = ",".join(["_qa"] + model_refs) 
    762813        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
    763814        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 
    764827 
    765828    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 

  • sasmodels/kernel_header.c

    r8698a0d r108e70e  
    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

    r6aee3ab r108e70e  
    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 
    3334//  INVALID(table) : test if the current point is feesible to calculate.  This 
    3435//      will be defined in the kernel definition file. 
    … …  
    469470  #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc) 
    470471  #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) 
    471479#endif 
    472480 
    … …  
    477485#if defined(CALL_IQ) || defined(CALL_IQ_A) 
    478486  #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 
    479517#else // !spherosymmetric projection 
    480518  // Grab the "view" angles (theta, phi, psi) from the initial parameter table. 

  • sasmodels/kernelpy.py

    r2d81cfe r108e70e  
    2626# pylint: enable=unused-import 
    2727 
     28logger = logging.getLogger(__name__) 
     29 
    2830class PyModel(KernelModel): 
    2931    """ 
    … …  
    3133    """ 
    3234    def __init__(self, model_info): 
    33         # Make sure Iq and Iqxy are available and vectorized 
     35        # Make sure Iq is available and vectorized 
    3436        _create_default_functions(model_info) 
    3537        self.info = model_info 
    3638        self.dtype = np.dtype('d') 
    37         logging.info("load python model " + self.info.name) 
     39        logger.info("load python model " + self.info.name) 
    3840 
    3941    def make_kernel(self, q_vectors): 
    4042        q_input = PyInput(q_vectors, dtype=F64) 
    41         kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq 
    42         return PyKernel(kernel, self.info, q_input) 
     43        return PyKernel(self.info, q_input) 
    4344 
    4445    def release(self): 
    … …  
    8990    Callable SAS kernel. 
    9091 
    91     *kernel* is the DllKernel object to call. 
     92    *kernel* is the kernel object to call. 
    9293 
    9394    *model_info* is the module information 
    … …  
    104105    Call :meth:`release` when done with the kernel instance. 
    105106    """ 
    106     def __init__(self, kernel, model_info, q_input): 
     107    def __init__(self, model_info, q_input): 
    107108        # type: (callable, ModelInfo, List[np.ndarray]) -> None 
    108109        self.dtype = np.dtype('d') 
    … …  
    110111        self.q_input = q_input 
    111112        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 
    162                         else (lambda: 1.0)) 
     161        self._volume = ((lambda: form_volume(*volume_args)) if form_volume else 
     162                        (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 
    263264    _create_vector_Iq(model_info) 
    264     _create_vector_Iqxy(model_info)  # call create_vector_Iq() first 
     265    _create_vector_Iqxy(model_info) 
    265266 
    266267 
    … …  
    280281        model_info.Iq = vector_Iq 
    281282 
     283 
    282284def _create_vector_Iqxy(model_info): 
    283285    """ 
    284286    Define Iqxy as a vector function if it exists, or default it from Iq(). 
    285287    """ 
    286     Iq, Iqxy = model_info.Iq, model_info.Iqxy 
     288    Iqxy = getattr(model_info, 'Iqxy', None) 
    287289    if callable(Iqxy): 
    288290        if not getattr(Iqxy, 'vectorized', False): 
    … …  
    295297            vector_Iqxy.vectorized = True 
    296298            model_info.Iqxy = vector_Iqxy 
    297     elif callable(Iq): 
     299    else: 
    298300        #print("defaulting Iqxy") 
    299301        # Iq is vectorized because create_vector_Iq was already called. 
     302        Iq = model_info.Iq 
    300303        def default_Iqxy(qx, qy, *args): 
    301304            """ 

  • sasmodels/modelinfo.py

    rdb03406 r1941ec6  
    4242 
    4343# assumptions about common parameters exist throughout the code, such as: 
    44 # (1) kernel functions Iq, Iqxy, form_volume, ... don't see them 
     44# (1) kernel functions Iq, Iqxy, Iqac, Iqabc, 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 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 
     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 
    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  
    396393    * *form_volume_parameters* is the list of parameters to the form_volume(...) 
    397394      function, with vector parameter p sent as p[]. 
    … …  
    448445        self.iq_parameters = [p for p in self.kernel_parameters 
    449446                              if p.type not in ('orientation', 'magnetic')] 
    450         # note: orientation no longer sent to Iqxy, so its the same as 
    451         #self.iqxy_parameters = [p for p in self.kernel_parameters 
    452         #                        if p.type != 'magnetic'] 
     447        self.orientation_parameters = [p for p in self.kernel_parameters 
     448                                       if p.type == 'orientation'] 
    453449        self.form_volume_parameters = [p for p in self.kernel_parameters 
    454450                                       if p.type == 'volume'] 
    … …  
    495491                if p.type != 'orientation': 
    496492                    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") 
    497495        if theta >= 0 and phi >= 0: 
     496            last_par = len(self.kernel_parameters) - 1 
    498497            if phi != theta+1: 
    499498                raise TypeError("phi must follow theta") 
    500499            if psi >= 0 and psi != phi+1: 
    501500                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") 
    502504        elif theta >= 0 or phi >= 0 or psi >= 0: 
    503505            raise TypeError("oriented shapes must have both theta and phi and maybe psi") 
    … …  
    720722 
    721723 
     724#: Set of variables defined in the model that might contain C code 
     725C_SYMBOLS = ['Imagnetic', 'Iq', 'Iqxy', 'Iqac', 'Iqabc', 'form_volume', 'c_code'] 
     726 
    722727def _find_source_lines(model_info, kernel_module): 
    723728    # type: (ModelInfo, ModuleType) -> None 
    … …  
    725730    Identify the location of the C source inside the model definition file. 
    726731 
    727     This code runs through the source of the kernel module looking for 
    728     lines that start with 'Iq', 'Iqxy' or 'form_volume'.  Clearly there are 
    729     all sorts of reasons why this might not work (e.g., code commented out 
    730     in a triple-quoted line block, code built using string concatenation, 
    731     or code defined in the branch of an 'if' block), but it should work 
    732     properly in the 95% case, and getting the incorrect line number will 
    733     be harmless. 
    734     """ 
    735     # Check if we need line numbers at all 
    736     if callable(model_info.Iq): 
    737         return None 
    738  
    739     if (model_info.Iq is None 
    740             and model_info.Iqxy is None 
    741             and model_info.Imagnetic is None 
    742             and model_info.form_volume is None): 
     732    This code runs through the source of the kernel module looking for lines 
     733    that contain C code (because it is a c function definition).  Clearly 
     734    there are all sorts of reasons why this might not work (e.g., code 
     735    commented out in a triple-quoted line block, code built using string 
     736    concatenation, code defined in the branch of an 'if' block, code imported 
     737    from another file), but it should work properly in the 95% case, and for 
     738    the remainder, getting the incorrect line number will merely be 
     739    inconvenient. 
     740    """ 
     741    # Only need line numbers if we are creating a C module and the C symbols 
     742    # are defined. 
     743    if (callable(model_info.Iq) 
     744            or not any(hasattr(model_info, s) for s in C_SYMBOLS)): 
    743745        return 
    744746 
    745     # find the defintion lines for the different code blocks 
     747    # load the module source if we can 
    746748    try: 
    747749        source = inspect.getsource(kernel_module) 
    748750    except IOError: 
    749751        return 
    750     for k, v in enumerate(source.split('\n')): 
    751         if v.startswith('Imagnetic'): 
    752             model_info._Imagnetic_line = k+1 
    753         elif v.startswith('Iqxy'): 
    754             model_info._Iqxy_line = k+1 
    755         elif v.startswith('Iq'): 
    756             model_info._Iq_line = k+1 
    757         elif v.startswith('form_volume'): 
    758             model_info._form_volume_line = k+1 
    759  
     752 
     753    # look for symbol at the start of the line 
     754    for lineno, line in enumerate(source.split('\n')): 
     755        for name in C_SYMBOLS: 
     756            if line.startswith(name): 
     757                # Add 1 since some compilers complain about "#line 0" 
     758                model_info.lineno[name] = lineno + 1 
     759                break 
    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 Iqxy functions will be created for python 
     768    Note: vectorized Iq and Iqac/Iqabc 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) 
    798799    # TODO: check the structure of the tests 
    799800    info.tests = getattr(kernel_module, 'tests', []) 
    … …  
    803804    info.Iq = getattr(kernel_module, 'Iq', None) # type: ignore 
    804805    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 
    805808    info.Imagnetic = getattr(kernel_module, 'Imagnetic', None) # type: ignore 
    806809    info.profile = getattr(kernel_module, 'profile', None) # type: ignore 
    … …  
    817820    info.hidden = getattr(kernel_module, 'hidden', None) # type: ignore 
    818821 
     822    if callable(info.Iq) and parameters.has_2d: 
     823        raise ValueError("oriented python models not supported") 
     824 
     825    info.lineno = {} 
    819826    _find_source_lines(info, kernel_module) 
    820827    try: 
    … …  
    835842 
    836843    The structure should be mostly static, other than the delayed definition 
    837     of *Iq* and *Iqxy* if they need to be defined. 
     844    of *Iq*, *Iqac* and *Iqabc* if they need to be defined. 
    838845    """ 
    839846    #: Full path to the file defining the kernel, if any. 
    … …  
    917924    structure_factor = None # type: bool 
    918925    #: List of C source files used to define the model.  The source files 
    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. 
     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. 
    924931    source = None           # type: List[str] 
    925932    #: The set of tests that must pass.  The format of the tests is described 
    … …  
    970977    #: include the decimal point. See :mod:`generate` for more details. 
    971978    Iq = None               # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
    972     #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`. 
    973     Iqxy = None             # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
     979    #: Returns *I(qab, qc, a, b, ...)*.  The interface follows :attr:`Iq`. 
     980    Iqac = None             # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
     981    #: Returns *I(qa, qb, qc, a, b, ...)*.  The interface follows :attr:`Iq`. 
     982    Iqabc = None            # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
    974983    #: Returns *I(qx, qy, a, b, ...)*.  The interface follows :attr:`Iq`. 
    975984    Imagnetic = None        # type: Union[None, str, Callable[[np.ndarray], np.ndarray]] 
    … …  
    987996    #: Returns a random parameter set for the model 
    988997    random = None           # type: Optional[Callable[[], Dict[str, float]]] 
    989  
    990     # line numbers within the python file for bits of C source, if defined 
    991     # NB: some compilers fail with a "#line 0" directive, so default to 1. 
    992     _Imagnetic_line = 1 
    993     _Iqxy_line = 1 
    994     _Iq_line = 1 
    995     _form_volume_line = 1 
    996  
     998    #: Line numbers for symbols defining C code 
     999    lineno = None           # type: Dict[str, int] 
    9971000 
    9981001    def __init__(self): 

  • sasmodels/models/_spherepy.py

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

  • sasmodels/models/barbell.c

    rbecded3 r108e70e  
    2323    const double qab_r = radius_bell*qab; // Q*R*sin(theta) 
    2424    double total = 0.0; 
    25     for (int i = 0; i < 76; i++){ 
    26         const double t = Gauss76Z[i]*zm + zb; 
     25    for (int i = 0; i < GAUSS_N; i++){ 
     26        const double t = GAUSS_Z[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 += Gauss76Wt[i] * Fq; 
     30        total += GAUSS_W[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 < 76; i++){ 
    76         const double alpha= Gauss76Z[i]*zm + zb; 
     75    for (int i = 0; i < GAUSS_N; i++){ 
     76        const double alpha= GAUSS_Z[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 += Gauss76Wt[i] * Aq * Aq * sin_alpha; 
     80        total += GAUSS_W[i] * Aq * Aq * sin_alpha; 
    8181    } 
    8282    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    9090 
    9191static double 
    92 Iqxy(double qab, double qc, 
     92Iqac(double qab, double qc, 
    9393    double sld, double solvent_sld, 
    9494    double radius_bell, double radius, double length) 

  • sasmodels/models/bcc_paracrystal.c

    rea60e08 r108e70e  
    8181 
    8282    double outer_sum = 0.0; 
    83     for(int i=0; i<150; i++) { 
     83    for(int i=0; i<GAUSS_N; i++) { 
    8484        double inner_sum = 0.0; 
    85         const double theta = Gauss150Z[i]*theta_m + theta_b; 
     85        const double theta = GAUSS_Z[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<150;j++) { 
    91             const double phi = Gauss150Z[j]*phi_m + phi_b; 
     90        for(int j=0;j<GAUSS_N;j++) { 
     91            const double phi = GAUSS_Z[j]*phi_m + phi_b; 
    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 += Gauss150Wt[j] * form; 
     97            inner_sum += GAUSS_W[j] * form; 
    9898        } 
    9999        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    100         outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     100        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    101101    } 
    102102    outer_sum *= theta_m; 
    … …  
    107107 
    108108 
    109 static double Iqxy(double qa, double qb, double qc, 
     109static double Iqabc(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

    rbecded3 r108e70e  
    3030    const double qab_r = radius_cap*qab; // Q*R*sin(theta) 
    3131    double total = 0.0; 
    32     for (int i=0; i<76 ;i++) { 
    33         const double t = Gauss76Z[i]*zm + zb; 
     32    for (int i=0; i<GAUSS_N; i++) { 
     33        const double t = GAUSS_Z[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 += Gauss76Wt[i] * Fq; 
     37        total += GAUSS_W[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<76 ;i++) { 
    98         const double theta = Gauss76Z[i]*zm + zb; 
     97    for (int i=0; i<GAUSS_N ;i++) { 
     98        const double theta = GAUSS_Z[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 += Gauss76Wt[i] * Aq * Aq * sin_theta; 
     105        total += GAUSS_W[i] * Aq * Aq * sin_theta; 
    106106    } 
    107107    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    115115 
    116116static double 
    117 Iqxy(double qab, double qc, 
     117Iqac(double qab, double qc, 
    118118    double sld, double solvent_sld, double radius, 
    119119    double radius_cap, double length) 

  • sasmodels/models/core_shell_bicelle.c

    rbecded3 r108e70e  
    5252 
    5353    double total = 0.0; 
    54     for(int i=0;i<N_POINTS_76;i++) { 
    55         double theta = (Gauss76Z[i] + 1.0)*uplim; 
     54    for(int i=0;i<GAUSS_N;i++) { 
     55        double theta = (GAUSS_Z[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 += Gauss76Wt[i]*fq*fq*sin_theta; 
     60        total += GAUSS_W[i]*fq*fq*sin_theta; 
    6161    } 
    6262 
    … …  
    6767 
    6868static double 
    69 Iqxy(double qab, double qc, 
     69Iqac(double qab, double qc, 
    7070    double radius, 
    7171    double thick_rim, 

  • sasmodels/models/core_shell_bicelle_elliptical.c

    r82592da r108e70e  
    3737    //initialize integral 
    3838    double outer_sum = 0.0; 
    39     for(int i=0;i<76;i++) { 
     39    for(int i=0;i<GAUSS_N;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 = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
    44         const double cos_theta = ( Gauss76Z[i] + 1.0 )/2.0; 
     43        //const double cos_theta = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
     44        const double cos_theta = ( GAUSS_Z[i] + 1.0 )/2.0; 
    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<76;j++) { 
     51        for(int j=0;j<GAUSS_N;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 = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    57             const double phi = ( Gauss76Z[j] +1.0)*M_PI_2; 
     56            //const double phi = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     57            const double phi = ( GAUSS_Z[j] +1.0)*M_PI_2; 
    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 += Gauss76Wt[j] * fq * fq; 
     63            inner_sum += GAUSS_W[j] * fq * fq; 
    6464        } 
    6565        //now calculate outer integral 
    66         outer_sum += Gauss76Wt[i] * inner_sum; 
     66        outer_sum += GAUSS_W[i] * inner_sum; 
    6767    } 
    6868 
    … …  
    7171 
    7272static double 
    73 Iqxy(double qa, double qb, double qc, 
     73Iqabc(double qa, double qb, double qc, 
    7474    double r_minor, 
    7575    double x_core, 

  • sasmodels/models/core_shell_bicelle_elliptical_belt_rough.c

    r82592da r108e70e  
    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<76;i++) { 
     49    for(int i=0;i<GAUSS_N;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 = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
    53         const double cos_alpha = ( Gauss76Z[i] + 1.0 )/2.0; 
     52        //const double cos_alpha = ( GAUSS_Z[i]*(vb-va) + va + vb )/2.0; 
     53        const double cos_alpha = ( GAUSS_Z[i] + 1.0 )/2.0; 
    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<76;j++) { 
     60        for(int j=0;j<GAUSS_N;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 = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    63             const double beta = ( Gauss76Z[j] +1.0)*M_PI_2; 
     62            //const double beta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     63            const double beta = ( GAUSS_Z[j] +1.0)*M_PI_2; 
    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 += Gauss76Wt[j] *square(dr1*si1*be1 + 
     69            inner_sum += GAUSS_W[j] *square(dr1*si1*be1 + 
    7070                                              dr2*si1*be2 + 
    7171                                              dr3*si2*be1); 
    7272        } 
    7373        //now calculate outer integral 
    74         outer_sum += Gauss76Wt[i] * inner_sum; 
     74        outer_sum += GAUSS_W[i] * inner_sum; 
    7575    } 
    7676 
    … …  
    7979 
    8080static double 
    81 Iqxy(double qa, double qb, double qc, 
     81Iqabc(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

    r110f69c r108e70e  
    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"], 
    151152    ["theta",       "degrees",    90.0, [-360, 360], "orientation", "cylinder axis to beam angle"], 
    152153    ["phi",         "degrees",    0,    [-360, 360], "orientation", "rotation about beam"], 
    153154    ["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

    rbecded3 r108e70e  
    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<76 ;i++) { 
     32    for (int i=0; i<GAUSS_N ;i++) { 
    3333        // translate a point in [-1,1] to a point in [0, pi/2] 
    34         //const double theta = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
     34        //const double theta = ( GAUSS_Z[i]*(upper-lower) + upper + lower )/2.0; 
    3535        double sin_theta, cos_theta; 
    36         const double theta = Gauss76Z[i]*M_PI_4 + M_PI_4; 
     36        const double theta = GAUSS_Z[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 += Gauss76Wt[i] * fq * fq * sin_theta; 
     42        total += GAUSS_W[i] * fq * fq * sin_theta; 
    4343    } 
    4444    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    4848 
    4949 
    50 double Iqxy(double qab, double qc, 
     50double Iqac(double qab, double qc, 
    5151    double core_sld, 
    5252    double shell_sld, 

  • sasmodels/models/core_shell_ellipsoid.c

    rbecded3 r108e70e  
    5959    const double b = 0.5; 
    6060    double total = 0.0;     //initialize intergral 
    61     for(int i=0;i<76;i++) { 
    62         const double cos_theta = Gauss76Z[i]*m + b; 
     61    for(int i=0;i<GAUSS_N;i++) { 
     62        const double cos_theta = GAUSS_Z[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 += Gauss76Wt[i] * fq * fq; 
     68        total += GAUSS_W[i] * fq * fq; 
    6969    } 
    7070    total *= m; 
    … …  
    7575 
    7676static double 
    77 Iqxy(double qab, double qc, 
     77Iqac(double qab, double qc, 
    7878    double radius_equat_core, 
    7979    double x_core, 

  • sasmodels/models/core_shell_parallelepiped.c

    r904cd9c r108e70e  
     1// Set OVERLAPPING to 1 in order to fill in the edges of the box, with 
     2// c endcaps and b overlapping a.  With the proper choice of parameters, 
     3// (setting rim slds to sld, core sld to solvent, rim thickness to thickness 
     4// and subtracting 2*thickness from length, this should match the hollow 
     5// rectangular prism.)  Set it to 0 for the documented behaviour. 
     6#define OVERLAPPING 0 
    17static double 
    28form_volume(double length_a, double length_b, double length_c, 
    39    double thick_rim_a, double thick_rim_b, double thick_rim_c) 
    410{ 
    5     //return length_a * length_b * length_c; 
    6     return length_a * length_b * length_c + 
    7            2.0 * thick_rim_a * length_b * length_c + 
    8            2.0 * thick_rim_b * length_a * length_c + 
    9            2.0 * thick_rim_c * length_a * length_b; 
     11    return 
     12#if OVERLAPPING 
     13        // Hollow rectangular prism only includes the volume of the shell 
     14        // so uncomment the next line when comparing.  Solid rectangular 
     15        // prism, or parallelepiped want filled cores, so comment when 
     16        // comparing. 
     17        //-length_a * length_b * length_c + 
     18        (length_a + 2.0*thick_rim_a) * 
     19        (length_b + 2.0*thick_rim_b) * 
     20        (length_c + 2.0*thick_rim_c); 
     21#else 
     22        length_a * length_b * length_c + 
     23        2.0 * thick_rim_a * length_b * length_c + 
     24        2.0 * length_a * thick_rim_b * length_c + 
     25        2.0 * length_a * length_b * thick_rim_c; 
     26#endif 
    1027} 
    1128 
    … …  
    2441    double thick_rim_c) 
    2542{ 
    26     // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c_scaled 
     43    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c 
    2744    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    28     //Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
     45    // Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
     46    // Code rewritten (PAK) 
    2947 
    30     const double mu = 0.5 * q * length_b; 
     48    const double half_q = 0.5*q; 
    3149 
    32     //calculate volume before rescaling (in original code, but not used) 
    33     //double vol = form_volume(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); 
    34     //double vol = length_a * length_b * length_c + 
    35     //       2.0 * thick_rim_a * length_b * length_c + 
    36     //       2.0 * thick_rim_b * length_a * length_c + 
    37     //       2.0 * thick_rim_c * length_a * length_b; 
     50    const double tA = length_a + 2.0*thick_rim_a; 
     51    const double tB = length_b + 2.0*thick_rim_b; 
     52    const double tC = length_c + 2.0*thick_rim_c; 
    3853 
    39     // Scale sides by B 
    40     const double a_scaled = length_a / length_b; 
    41     const double c_scaled = length_c / length_b; 
    42  
    43     double ta = a_scaled + 2.0*thick_rim_a/length_b; // incorrect ta = (a_scaled + 2.0*thick_rim_a)/length_b; 
    44     double tb = 1+ 2.0*thick_rim_b/length_b; // incorrect tb = (a_scaled + 2.0*thick_rim_b)/length_b; 
    45     double tc = c_scaled + 2.0*thick_rim_c/length_b; //not present 
    46  
    47     double Vin = length_a * length_b * length_c; 
    48     //double Vot = (length_a * length_b * length_c + 
    49     //            2.0 * thick_rim_a * length_b * length_c + 
    50     //            2.0 * length_a * thick_rim_b * length_c + 
    51     //            2.0 * length_a * length_b * thick_rim_c); 
    52     double V1 = (2.0 * thick_rim_a * length_b * length_c);    // incorrect V1 (aa*bb*cc+2*ta*bb*cc) 
    53     double V2 = (2.0 * length_a * thick_rim_b * length_c);    // incorrect V2(aa*bb*cc+2*aa*tb*cc) 
    54     double V3 = (2.0 * length_a * length_b * thick_rim_c);    //not present 
    55  
    56     // Scale factors (note that drC is not used later) 
    57     const double drho0 = (core_sld-solvent_sld); 
    58     const double drhoA = (arim_sld-solvent_sld); 
    59     const double drhoB = (brim_sld-solvent_sld); 
    60     const double drhoC = (crim_sld-solvent_sld);  // incorrect const double drC_Vot = (crim_sld-solvent_sld)*Vot; 
    61  
    62  
    63     // Precompute scale factors for combining cross terms from the shape 
    64     const double scale23 = drhoA*V1; 
    65     const double scale14 = drhoB*V2; 
    66     const double scale24 = drhoC*V3; 
    67     const double scale11 = drho0*Vin; 
    68     const double scale12 = drho0*Vin - scale23 - scale14 - scale24; 
     54    // Scale factors 
     55    const double dr0 = (core_sld-solvent_sld); 
     56    const double drA = (arim_sld-solvent_sld); 
     57    const double drB = (brim_sld-solvent_sld); 
     58    const double drC = (crim_sld-solvent_sld); 
    6959 
    7060    // outer integral (with gauss points), integration limits = 0, 1 
    71     double outer_total = 0; //initialize integral 
     61    double outer_sum = 0; //initialize integral 
     62    for( int i=0; i<GAUSS_N; i++) { 
     63        const double cos_alpha = 0.5 * ( GAUSS_Z[i] + 1.0 ); 
     64        const double mu = half_q * sqrt(1.0-cos_alpha*cos_alpha); 
    7265 
    73     for( int i=0; i<76; i++) { 
    74         double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
    75         double mu_proj = mu * sqrt(1.0-sigma*sigma); 
     66        // inner integral (with gauss points), integration limits = 0, pi/2 
     67        const double siC = length_c * sas_sinx_x(length_c * cos_alpha * half_q); 
     68        const double siCt = tC * sas_sinx_x(tC * cos_alpha * half_q); 
     69        double inner_sum = 0.0; 
     70        for(int j=0; j<GAUSS_N; j++) { 
     71            const double beta = 0.5 * ( GAUSS_Z[j] + 1.0 ); 
     72            double sin_beta, cos_beta; 
     73            SINCOS(M_PI_2*beta, sin_beta, cos_beta); 
     74            const double siA = length_a * sas_sinx_x(length_a * mu * sin_beta); 
     75            const double siB = length_b * sas_sinx_x(length_b * mu * cos_beta); 
     76            const double siAt = tA * sas_sinx_x(tA * mu * sin_beta); 
     77            const double siBt = tB * sas_sinx_x(tB * mu * cos_beta); 
    7678 
    77         // inner integral (with gauss points), integration limits = 0, 1 
    78         double inner_total = 0.0; 
    79         double inner_total_crim = 0.0; 
    80         for(int j=0; j<76; j++) { 
    81             const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
    82             double sin_uu, cos_uu; 
    83             SINCOS(M_PI_2*uu, sin_uu, cos_uu); 
    84             const double si1 = sas_sinx_x(mu_proj * sin_uu * a_scaled); 
    85             const double si2 = sas_sinx_x(mu_proj * cos_uu ); 
    86             const double si3 = sas_sinx_x(mu_proj * sin_uu * ta); 
    87             const double si4 = sas_sinx_x(mu_proj * cos_uu * tb); 
     79#if OVERLAPPING 
     80            const double f = dr0*siA*siB*siC 
     81                + drA*(siAt-siA)*siB*siC 
     82                + drB*siAt*(siBt-siB)*siC 
     83                + drC*siAt*siBt*(siCt-siC); 
     84#else 
     85            const double f = dr0*siA*siB*siC 
     86                + drA*(siAt-siA)*siB*siC 
     87                + drB*siA*(siBt-siB)*siC 
     88                + drC*siA*siB*(siCt-siC); 
     89#endif 
    8890 
    89             // Expression in libCylinder.c (neither drC nor Vot are used) 
    90             const double form = scale12*si1*si2 + scale23*si2*si3 + scale14*si1*si4; 
    91             const double form_crim = scale11*si1*si2; 
    92  
    93             //  correct FF : sum of square of phase factors 
    94             inner_total += Gauss76Wt[j] * form * form; 
    95             inner_total_crim += Gauss76Wt[j] * form_crim * form_crim; 
     91            inner_sum += GAUSS_W[j] * f * f; 
    9692        } 
    97         inner_total *= 0.5; 
    98         inner_total_crim *= 0.5; 
     93        inner_sum *= 0.5; 
    9994        // now sum up the outer integral 
    100         const double si = sas_sinx_x(mu * c_scaled * sigma); 
    101         const double si_crim = sas_sinx_x(mu * tc * sigma); 
    102         outer_total += Gauss76Wt[i] * (inner_total * si * si + inner_total_crim * si_crim * si_crim); 
     95        outer_sum += GAUSS_W[i] * inner_sum; 
    10396    } 
    104     outer_total *= 0.5; 
     97    outer_sum *= 0.5; 
    10598 
    10699    //convert from [1e-12 A-1] to [cm-1] 
    107     return 1.0e-4 * outer_total; 
     100    return 1.0e-4 * outer_sum; 
    108101} 
    109102 
    110103static double 
    111 Iqxy(double qa, double qb, double qc, 
     104Iqabc(double qa, double qb, double qc, 
    112105    double core_sld, 
    113106    double arim_sld, 
    … …  
    128121    const double drC = crim_sld-solvent_sld; 
    129122 
    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  
    140123    // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    141124    // the scaling by B. 
    142     double ta = length_a + 2.0*thick_rim_a; 
    143     double tb = length_b + 2.0*thick_rim_b; 
    144     double tc = length_c + 2.0*thick_rim_c; 
    145     //handle arg=0 separately, as sin(t)/t -> 1 as t->0 
    146     double siA = sas_sinx_x(0.5*length_a*qa); 
    147     double siB = sas_sinx_x(0.5*length_b*qb); 
    148     double siC = sas_sinx_x(0.5*length_c*qc); 
    149     double siAt = sas_sinx_x(0.5*ta*qa); 
    150     double siBt = sas_sinx_x(0.5*tb*qb); 
    151     double siCt = sas_sinx_x(0.5*tc*qc); 
     125    const double tA = length_a + 2.0*thick_rim_a; 
     126    const double tB = length_b + 2.0*thick_rim_b; 
     127    const double tC = length_c + 2.0*thick_rim_c; 
     128    const double siA = length_a*sas_sinx_x(0.5*length_a*qa); 
     129    const double siB = length_b*sas_sinx_x(0.5*length_b*qb); 
     130    const double siC = length_c*sas_sinx_x(0.5*length_c*qc); 
     131    const double siAt = tA*sas_sinx_x(0.5*tA*qa); 
     132    const double siBt = tB*sas_sinx_x(0.5*tB*qb); 
     133    const double siCt = tC*sas_sinx_x(0.5*tC*qc); 
    152134 
    153  
    154     // f uses Vin, V1, V2, and V3 and it seems to have more sense than the value computed 
    155     // in the 1D code, but should be checked! 
    156     double f = ( dr0*siA*siB*siC*Vin 
    157                + drA*(siAt-siA)*siB*siC*V1 
    158                + drB*siA*(siBt-siB)*siC*V2 
    159                + drC*siA*siB*(siCt-siC)*V3); 
     135#if OVERLAPPING 
     136    const double f = dr0*siA*siB*siC 
     137        + drA*(siAt-siA)*siB*siC 
     138        + drB*siAt*(siBt-siB)*siC 
     139        + drC*siAt*siBt*(siCt-siC); 
     140#else 
     141    const double f = dr0*siA*siB*siC 
     142        + drA*(siAt-siA)*siB*siC 
     143        + drB*siA*(siBt-siB)*siC 
     144        + drC*siA*siB*(siCt-siC); 
     145#endif 
    160146 
    161147    return 1.0e-4 * f * f; 

  • sasmodels/models/core_shell_parallelepiped.py

    r2d81cfe r10ee838  
    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. However at this time the 1D 
    8 calculation does **NOT** actually calculate a c face rim despite the presence 
    9 of the parameter. Some other aspects of the 1D calculation may be wrong. 
    10  
    11 .. note:: 
    12    This model was originally ported from NIST IGOR macros. However, it is not 
    13    yet fully understood by the SasView developers and is currently under review. 
     7"rim" can be different on each (pair) of faces. 
    148 
    159The form factor is normalized by the particle volume $V$ such that 
    … …  
    2115where $\langle \ldots \rangle$ is an average over all possible orientations 
    2216of the rectangular solid. 
    23  
    2417 
    2518The function calculated is the form factor of the rectangular solid below. 
    … …  
    4134    V = ABC + 2t_ABC + 2t_BAC + 2t_CAB 
    4235 
    43 **meaning that there are "gaps" at the corners of the solid.**  Again note that 
    44 $t_C = 0$ currently. 
     36**meaning that there are "gaps" at the corners of the solid.** 
    4537 
    4638The intensity calculated follows the :ref:`parallelepiped` model, with the 
    4739core-shell intensity being calculated as the square of the sum of the 
    48 amplitudes of the core and shell, in the same manner as a core-shell model. 
    49  
    50 .. math:: 
    51  
    52     F_{a}(Q,\alpha,\beta)= 
    53     \left[\frac{\sin(\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha \sin\beta) 
    54                }{\tfrac{1}{2}Q(L_A+2t_A)\sin\alpha\sin\beta} 
    55     - \frac{\sin(\tfrac{1}{2}QL_A\sin\alpha \sin\beta) 
    56            }{\tfrac{1}{2}QL_A\sin\alpha \sin\beta} \right] 
    57     \left[\frac{\sin(\tfrac{1}{2}QL_B\sin\alpha \sin\beta) 
    58                }{\tfrac{1}{2}QL_B\sin\alpha \sin\beta} \right] 
    59     \left[\frac{\sin(\tfrac{1}{2}QL_C\sin\alpha \sin\beta) 
    60                }{\tfrac{1}{2}QL_C\sin\alpha \sin\beta} \right] 
    61  
    62 .. note:: 
    63  
    64     Why does t_B not appear in the above equation? 
    65     For the calculation of the form factor to be valid, the sides of the solid 
    66     MUST (perhaps not any more?) be chosen such that** $A < B < C$. 
    67     If this inequality is not satisfied, the model will not report an error, 
    68     but the calculation will not be correct and thus the result wrong. 
     40amplitudes of the core and the slabs on the edges. 
     41 
     42the scattering amplitude is computed for a particular orientation of the 
     43core-shell parallelepiped with respect to the scattering vector and then 
     44averaged over all possible orientations, where $\alpha$ is the angle between 
     45the $z$ axis and the $C$ axis of the parallelepiped, $\beta$ is 
     46the angle between projection of the particle in the $xy$ detector plane 
     47and the $y$ axis. 
     48 
     49.. math:: 
     50 
     51    F(Q) 
     52    &= (\rho_\text{core}-\rho_\text{solvent}) 
     53       S(Q_A, A) S(Q_B, B) S(Q_C, C) \\ 
     54    &+ (\rho_\text{A}-\rho_\text{solvent}) 
     55        \left[S(Q_A, A+2t_A) - S(Q_A, Q)\right] S(Q_B, B) S(Q_C, C) \\ 
     56    &+ (\rho_\text{B}-\rho_\text{solvent}) 
     57        S(Q_A, A) \left[S(Q_B, B+2t_B) - S(Q_B, B)\right] S(Q_C, C) \\ 
     58    &+ (\rho_\text{C}-\rho_\text{solvent}) 
     59        S(Q_A, A) S(Q_B, B) \left[S(Q_C, C+2t_C) - S(Q_C, C)\right] 
     60 
     61with 
     62 
     63.. math:: 
     64 
     65    S(Q, L) = L \frac{\sin \tfrac{1}{2} Q L}{\tfrac{1}{2} Q L} 
     66 
     67and 
     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 
     76where $\rho_\text{core}$, $\rho_\text{A}$, $\rho_\text{B}$ and $\rho_\text{C}$ 
     77are the scattering length of the parallelepiped core, and the rectangular 
     78slabs of thickness $t_A$, $t_B$ and $t_C$, respectively. $\rho_\text{solvent}$ 
     79is the scattering length of the solvent. 
    6980 
    7081FITTING NOTES 
     82~~~~~~~~~~~~~ 
     83 
    7184If the scale is set equal to the particle volume fraction, $\phi$, the returned 
    72 value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. 
    73 However, **no interparticle interference effects are included in this 
    74 calculation.** 
     85value is the scattered intensity per unit volume, $I(q) = \phi P(q)$. However, 
     86**no interparticle interference effects are included in this calculation.** 
    7587 
    7688There are many parameters in this model. Hold as many fixed as possible with 
    7789known values, or you will certainly end up at a solution that is unphysical. 
    7890 
    79 Constraints must be applied during fitting to ensure that the inequality 
    80 $A < B < C$ is not violated. The calculation will not report an error, 
    81 but the results will not be correct. 
    82  
    8391The returned value is in units of |cm^-1|, on absolute scale. 
    8492 
    8593NB: The 2nd virial coefficient of the core_shell_parallelepiped is calculated 
    8694based on the the averaged effective radius $(=\sqrt{(A+2t_A)(B+2t_B)/\pi})$ 
    87 and length $(C+2t_C)$ values, after appropriately 
    88 sorting the three dimensions to give an oblate or prolate particle, to give an 
    89 effective radius, for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
     95and length $(C+2t_C)$ values, after appropriately sorting the three dimensions 
     96to give an oblate or prolate particle, to give an effective radius, 
     97for $S(Q)$ when $P(Q) * S(Q)$ is applied. 
    9098 
    9199For 2d data the orientation of the particle is required, described using 
    92 angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
    93 of the calculation and angular dispersions see :ref:`orientation` . 
     100angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further 
     101details of the calculation and angular dispersions see :ref:`orientation`. 
    94102The angle $\Psi$ is the rotational angle around the *long_c* axis. For example, 
    95103$\Psi = 0$ when the *short_b* axis is parallel to the *x*-axis of the detector. 
     104 
     105For 2d, constraints must be applied during fitting to ensure that the 
     106inequality $A < B < C$ is not violated, and hence the correct definition 
     107of angles is preserved. The calculation will not report an error, 
     108but the results may be not correct. 
    96109 
    97110.. figure:: img/parallelepiped_angle_definition.png 
    … …  
    114127    Equations (1), (13-14). (in German) 
    115128.. [#] D Singh (2009). *Small angle scattering studies of self assembly in 
    116    lipid mixtures*, John's Hopkins University Thesis (2009) 223-225. `Available 
     129   lipid mixtures*, Johns Hopkins University Thesis (2009) 223-225. `Available 
    117130   from Proquest <http://search.proquest.com/docview/304915826?accountid 
    118131   =26379>`_ 
    … …  
    175188        Return equivalent radius (ER) 
    176189    """ 
    177  
    178     # surface average radius (rough approximation) 
    179     surf_rad = sqrt((length_a + 2.0*thick_rim_a) * (length_b + 2.0*thick_rim_b) / pi) 
    180  
    181     height = length_c + 2.0*thick_rim_c 
    182  
    183     ddd = 0.75 * surf_rad * (2 * surf_rad * height + (height + surf_rad) * (height + pi * surf_rad)) 
    184     return 0.5 * (ddd) ** (1. / 3.) 
     190    from .parallelepiped import ER as ER_p 
     191 
     192    a = length_a + 2*thick_rim_a 
     193    b = length_b + 2*thick_rim_b 
     194    c = length_c + 2*thick_rim_c 
     195    return ER_p(a, b, c) 
    185196 
    186197# VR defaults to 1.0 
    … …  
    216227            psi_pd=10, psi_pd_n=1) 
    217228 
    218 # rkh 7/4/17 add random unit test for 2d, note make all params different, 2d values not tested against other codes or models 
     229# rkh 7/4/17 add random unit test for 2d, note make all params different, 
     230# 2d values not tested against other codes or models 
    219231if 0:  # pak: model rewrite; need to update tests 
    220232    qx, qy = 0.2 * cos(pi/6.), 0.2 * sin(pi/6.) 

  • sasmodels/models/cylinder.c

    rbecded3 r108e70e  
    2121 
    2222    double total = 0.0; 
    23     for (int i=0; i<76 ;i++) { 
    24         const double theta = Gauss76Z[i]*zm + zb; 
     23    for (int i=0; i<GAUSS_N ;i++) { 
     24        const double theta = GAUSS_Z[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 += Gauss76Wt[i] * form * form * sin_theta; 
     29        total += GAUSS_W[i] * form * form * sin_theta; 
    3030    } 
    3131    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    4545 
    4646static double 
    47 Iqxy(double qab, double qc, 
     47Iqac(double qab, double qc, 
    4848    double sld, 
    4949    double solvent_sld, 

  • sasmodels/models/ellipsoid.c

    rbecded3 r108e70e  
    2222 
    2323    // translate a point in [-1,1] to a point in [0, 1] 
    24     // const double u = Gauss76Z[i]*(upper-lower)/2 + (upper+lower)/2; 
     24    // const double u = GAUSS_Z[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<76;i++) { 
    29         const double u = Gauss76Z[i]*zm + zb; 
     28    for (int i=0;i<GAUSS_N;i++) { 
     29        const double u = GAUSS_Z[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 += Gauss76Wt[i] * f * f; 
     32        total += GAUSS_W[i] * f * f; 
    3333    } 
    3434    // translate dx in [-1,1] to dx in [lower,upper] 
    … …  
    3939 
    4040static double 
    41 Iqxy(double qab, double qc, 
     41Iqac(double qab, double qc, 
    4242    double sld, 
    4343    double sld_solvent, 

  • sasmodels/models/elliptical_cylinder.c

    r82592da r108e70e  
    2222    //initialize integral 
    2323    double outer_sum = 0.0; 
    24     for(int i=0;i<76;i++) { 
     24    for(int i=0;i<GAUSS_N;i++) { 
    2525        //setup inner integral over the ellipsoidal cross-section 
    26         const double cos_val = ( Gauss76Z[i]*(vb-va) + va + vb )/2.0; 
     26        const double cos_val = ( GAUSS_Z[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<76;j++) { 
    31             //20 gauss points for the inner integral, increase to 76, RKH 6Nov2017 
    32             const double theta = ( Gauss76Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
     30        for(int j=0;j<GAUSS_N;j++) { 
     31            const double theta = ( GAUSS_Z[j]*(vbj-vaj) + vaj + vbj )/2.0; 
    3332            const double r = sin_val*sqrt(rA - rB*cos(theta)); 
    3433            const double be = sas_2J1x_x(q*r); 
    35             inner_sum += Gauss76Wt[j] * be * be; 
     34            inner_sum += GAUSS_W[j] * be * be; 
    3635        } 
    3736        //now calculate the value of the inner integral 
    … …  
    4039        //now calculate outer integral 
    4140        const double si = sas_sinx_x(q*0.5*length*cos_val); 
    42         outer_sum += Gauss76Wt[i] * inner_sum * si * si; 
     41        outer_sum += GAUSS_W[i] * inner_sum * si * si; 
    4342    } 
    4443    outer_sum *= 0.5*(vb-va); 
    … …  
    5554 
    5655static double 
    57 Iqxy(double qa, double qb, double qc, 
     56Iqabc(double qa, double qb, double qc, 
    5857     double radius_minor, double r_ratio, double length, 
    5958     double sld, double solvent_sld) 

  • sasmodels/models/elliptical_cylinder.py

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

  • sasmodels/models/fcc_paracrystal.c

    rf728001 r108e70e  
    5353 
    5454    double outer_sum = 0.0; 
    55     for(int i=0; i<150; i++) { 
     55    for(int i=0; i<GAUSS_N; i++) { 
    5656        double inner_sum = 0.0; 
    57         const double theta = Gauss150Z[i]*theta_m + theta_b; 
     57        const double theta = GAUSS_Z[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<150;j++) { 
    63             const double phi = Gauss150Z[j]*phi_m + phi_b; 
     62        for(int j=0;j<GAUSS_N;j++) { 
     63            const double phi = GAUSS_Z[j]*phi_m + phi_b; 
    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 += Gauss150Wt[j] * form; 
     69            inner_sum += GAUSS_W[j] * form; 
    7070        } 
    7171        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    72         outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     72        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    7373    } 
    7474    outer_sum *= theta_m; 
    … …  
    8080 
    8181 
    82 static double Iqxy(double qa, double qb, double qc, 
     82static double Iqabc(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

    r592343f r74768cb  
    1717    double sum=0.0; 
    1818 
    19     for(int i=0;i<N_POINTS_76;i++) { 
    20         const double zi = ( Gauss76Z[i] + 1.0 )*M_PI_4; 
     19    for(int i=0;i<GAUSS_N;i++) { 
     20        const double zi = ( GAUSS_Z[i] + 1.0 )*M_PI_4; 
    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 += Gauss76Wt[i] * yyy * yyy; 
     25        sum += GAUSS_W[i] * yyy * yyy; 
    2626    } 
    2727    sum *= 0.5; 

  • sasmodels/models/hollow_cylinder.c

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

  • sasmodels/models/hollow_rectangular_prism.c

    r1e7b0db0 r108e70e  
    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      
     39 
    4040    double outer_sum = 0.0; 
    41     for(int i=0; i<76; i++) { 
     41    for(int i=0; i<GAUSS_N; i++) { 
    4242 
    43         const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
     43        const double theta = 0.5 * ( GAUSS_Z[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<76; j++) { 
     51        for(int j=0; j<GAUSS_N; j++) { 
    5252 
    53             const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
     53            const double phi = 0.5 * ( GAUSS_Z[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 += Gauss76Wt[j] * square(AP1-AP2); 
     68            inner_sum += GAUSS_W[j] * square(AP1-AP2); 
    6969        } 
    7070        inner_sum *= 0.5 * (v2b-v2a); 
    7171 
    72         outer_sum += Gauss76Wt[i] * inner_sum * sin(theta); 
     72        outer_sum += GAUSS_W[i] * inner_sum * sin(theta); 
    7373    } 
    7474    outer_sum *= 0.5*(v1b-v1a); 
    … …  
    8484    return 1.0e-4 * delrho * delrho * form; 
    8585} 
     86 
     87double 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 r0e55afe  
    55This model provides the form factor, $P(q)$, for a hollow rectangular 
    66parallelepiped with a wall of thickness $\Delta$. 
    7 It computes only the 1D scattering, not the 2D. 
     7 
    88 
    99Definition 
    … …  
    6666(which is unitless). 
    6767 
    68 **The 2D scattering intensity is not computed by this model.** 
     68For 2d data the orientation of the particle is required, described using 
     69angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
     70of the calculation and angular dispersions see :ref:`orientation` . 
     71The 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 
     74For 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, 
     76but 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. 
    6989 
    7090 
    … …  
    113133              ["thickness", "Ang", 1, [0, inf], "volume", 
    114134               "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"], 
    115141             ] 
    116142 

  • sasmodels/models/hollow_rectangular_prism_thin_walls.c

    rab2aea8 r74768cb  
    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<76; i++) { 
    34         const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
     33    for(int i=0; i<GAUSS_N; i++) { 
     34        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
    3535 
    3636        double sin_theta, cos_theta; 
    … …  
    4444 
    4545        double inner_sum = 0.0; 
    46         for(int j=0; j<76; j++) { 
    47             const double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
     46        for(int j=0; j<GAUSS_N; j++) { 
     47            const double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
    4848 
    4949            double sin_phi, cos_phi; 
    … …  
    6262                * ( cos_a*sin_b/cos_phi + cos_b*sin_a/sin_phi ); 
    6363 
    64             inner_sum += Gauss76Wt[j] * square(AL+AT); 
     64            inner_sum += GAUSS_W[j] * square(AL+AT); 
    6565        } 
    6666 
    6767        inner_sum *= 0.5 * (v2b-v2a); 
    68         outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 
     68        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    6969    } 
    7070 

  • sasmodels/models/lib/gauss150.c

    r994d77f r74768cb  
    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 
    919 
    1020// Gaussians 
    11 constant double Gauss150Z[150]={ 
     21constant double Gauss150Z[152]={ 
    1222        -0.9998723404457334, 
    1323        -0.9993274305065947, 
    … …  
    159169        0.9983473449340834, 
    160170        0.9993274305065947, 
    161         0.9998723404457334 
     171        0.9998723404457334, 
     172        0., 
     173        0. 
    162174}; 
    163175 
    164 constant double Gauss150Wt[150]={ 
     176constant double Gauss150Wt[152]={ 
    165177        0.0003276086705538, 
    166178        0.0007624720924706, 
    … …  
    312324        0.0011976474864367, 
    313325        0.0007624720924706, 
    314         0.0003276086705538 
     326        0.0003276086705538, 
     327        0., 
     328        0. 
    315329}; 

  • sasmodels/models/lib/gauss20.c

    r994d77f r74768cb  
    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 
    917 
    1018// Gaussians 

  • sasmodels/models/lib/gauss76.c

    r66d119f r74768cb  
    77 * 
    88 */ 
    9 #define N_POINTS_76 76 
     9 #ifdef GAUSS_N 
     10 # undef GAUSS_N 
     11 # undef GAUSS_Z 
     12 # undef GAUSS_W 
     13 #endif 
     14 #define GAUSS_N 76 
     15 #define GAUSS_Z Gauss76Z 
     16 #define GAUSS_W Gauss76Wt 
    1017 
    1118// Gaussians 
    12 constant double Gauss76Wt[N_POINTS_76]={ 
     19constant double Gauss76Wt[76]={ 
    1320        .00126779163408536,             //0 
    1421        .00294910295364247, 
    … …  
    8996}; 
    9097 
    91 constant double Gauss76Z[N_POINTS_76]={ 
     98constant double Gauss76Z[76]={ 
    9299        -.999505948362153,              //0 
    93100        -.997397786355355, 

  • sasmodels/models/line.py

    r2d81cfe r108e70e  
    5757Iq.vectorized = True # Iq accepts an array of q values 
    5858 
     59 
    5960def Iqxy(qx, qy, *args): 
    6061    """ 
    … …  
    6970 
    7071Iqxy.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#""" 
    7180 
    7281def random(): 

  • sasmodels/models/parallelepiped.c

    r9b7b23f r108e70e  
    2323    double outer_total = 0; //initialize integral 
    2424 
    25     for( int i=0; i<76; i++) { 
    26         const double sigma = 0.5 * ( Gauss76Z[i] + 1.0 ); 
     25    for( int i=0; i<GAUSS_N; i++) { 
     26        const double sigma = 0.5 * ( GAUSS_Z[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<76; j++) { 
    33             const double uu = 0.5 * ( Gauss76Z[j] + 1.0 ); 
     32        for(int j=0; j<GAUSS_N; j++) { 
     33            const double uu = 0.5 * ( GAUSS_Z[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 += Gauss76Wt[j] * square(si1 * si2); 
     38            inner_total += GAUSS_W[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 += Gauss76Wt[i] * inner_total * si * si; 
     43        outer_total += GAUSS_W[i] * inner_total * si * si; 
    4444    } 
    4545    outer_total *= 0.5; 
    … …  
    5353 
    5454static double 
    55 Iqxy(double qa, double qb, double qc, 
     55Iqabc(double qa, double qb, double qc, 
    5656    double sld, 
    5757    double solvent_sld, 

  • sasmodels/models/pringle.c

    r1e7b0db0 r74768cb  
    2929    double sumC = 0.0;          // initialize integral 
    3030    double r; 
    31     for (int i=0; i < 76; i++) { 
    32         r = Gauss76Z[i]*zm + zb; 
     31    for (int i=0; i < GAUSS_N; i++) { 
     32        r = GAUSS_Z[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 = Gauss76Wt[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs); 
     37        double y = GAUSS_W[i] * r * sas_JN(n, beta*qrrc) * sas_JN(2*n, qrs); 
    3838        double S, C; 
    3939        SINCOS(alpha*qrrc, S, C); 
    … …  
    8686 
    8787    double sum = 0.0; 
    88     for (int i = 0; i < 76; i++) { 
    89         double psi = Gauss76Z[i]*zm + zb; 
     88    for (int i = 0; i < GAUSS_N; i++) { 
     89        double psi = GAUSS_Z[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 += Gauss76Wt[i] * pringle_kernel; 
     95        sum += GAUSS_W[i] * pringle_kernel; 
    9696    } 
    9797 

  • sasmodels/models/rectangular_prism.c

    r1e7b0db0 r108e70e  
    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<76; i++) { 
    31         const double theta = 0.5 * ( Gauss76Z[i]*(v1b-v1a) + v1a + v1b ); 
     30    for(int i=0; i<GAUSS_N; i++) { 
     31        const double theta = 0.5 * ( GAUSS_Z[i]*(v1b-v1a) + v1a + v1b ); 
    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<76; j++) { 
    39             double phi = 0.5 * ( Gauss76Z[j]*(v2b-v2a) + v2a + v2b ); 
     38        for(int j=0; j<GAUSS_N; j++) { 
     39            double phi = 0.5 * ( GAUSS_Z[j]*(v2b-v2a) + v2a + v2b ); 
    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 += Gauss76Wt[j] * AP * AP; 
     47            inner_sum += GAUSS_W[j] * AP * AP; 
    4848        } 
    4949        inner_sum = 0.5 * (v2b-v2a) * inner_sum; 
    50         outer_sum += Gauss76Wt[i] * inner_sum * sin_theta; 
     50        outer_sum += GAUSS_W[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 
     73double 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 r0e55afe  
    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). 
    14 Note also that, contrary to :ref:`parallelepiped`, it does not compute 
    15 the 2D scattering. 
    1614 
    1715 
    … …  
    2624that reference), with $\theta$ corresponding to $\alpha$ in that paper, 
    2725and not to the usual convention used for example in the 
    28 :ref:`parallelepiped` model. As the present model does not compute 
    29 the 2D scattering, this has no further consequences. 
     26:ref:`parallelepiped` model. 
    3027 
    3128In this model the scattering from a massive parallelepiped with an 
    … …  
    6562units) *scale* represents the volume fraction (which is unitless). 
    6663 
    67 **The 2D scattering intensity is not computed by this model.** 
     64For 2d data the orientation of the particle is required, described using 
     65angles $\theta$, $\phi$ and $\Psi$ as in the diagrams below, for further details 
     66of the calculation and angular dispersions see :ref:`orientation` . 
     67The 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 
     70For 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, 
     72but 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 
    6886 
    6987 
    … …  
    108126              ["c2a_ratio", "", 1, [0, inf], "volume", 
    109127               "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"], 
    110134             ] 
    111135 

  • sasmodels/models/sc_paracrystal.c

    rf728001 r108e70e  
    5454 
    5555    double outer_sum = 0.0; 
    56     for(int i=0; i<150; i++) { 
     56    for(int i=0; i<GAUSS_N; i++) { 
    5757        double inner_sum = 0.0; 
    58         const double theta = Gauss150Z[i]*theta_m + theta_b; 
     58        const double theta = GAUSS_Z[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<150;j++) { 
    64             const double phi = Gauss150Z[j]*phi_m + phi_b; 
     63        for(int j=0;j<GAUSS_N;j++) { 
     64            const double phi = GAUSS_Z[j]*phi_m + phi_b; 
    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 += Gauss150Wt[j] * form; 
     70            inner_sum += GAUSS_W[j] * form; 
    7171        } 
    7272        inner_sum *= phi_m;  // sum(f(x)dx) = sum(f(x)) dx 
    73         outer_sum += Gauss150Wt[i] * inner_sum * sin_theta; 
     73        outer_sum += GAUSS_W[i] * inner_sum * sin_theta; 
    7474    } 
    7575    outer_sum *= theta_m; 
    … …  
    8282 
    8383static double 
    84 Iqxy(double qa, double qb, double qc, 
     84Iqabc(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

    rbecded3 r108e70e  
    8181    double halfheight = 0.5*thick_core; 
    8282 
    83     for(int i=0; i<N_POINTS_76; i++) { 
    84         double zi = (Gauss76Z[i] + 1.0)*M_PI_4; 
     83    for(int i=0; i<GAUSS_N; i++) { 
     84        double zi = (GAUSS_Z[i] + 1.0)*M_PI_4; 
    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 += Gauss76Wt[i] * yyy * sin_alpha; 
     97        summ += GAUSS_W[i] * yyy * sin_alpha; 
    9898    } 
    9999 
    … …  
    142142 
    143143static double 
    144 Iqxy(double qab, double qc, 
     144Iqac(double qab, double qc, 
    145145    double thick_core, 
    146146    double thick_layer, 

  • sasmodels/models/triaxial_ellipsoid.c

    rbecded3 r108e70e  
    2121    const double zb = M_PI_4; 
    2222    double outer = 0.0; 
    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; 
     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; 
    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<76;j++) { 
     31        for (int j=0;j<GAUSS_N;j++) { 
    3232            // translate a point in [-1,1] to a point in [0, 1] 
    33             const double usq = square(Gauss76Z[j]*um + ub); 
     33            const double usq = square(GAUSS_Z[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 += Gauss76Wt[j] * fq * fq; 
     36            inner += GAUSS_W[j] * fq * fq; 
    3737        } 
    38         outer += Gauss76Wt[i] * inner;  // correcting for dx later 
     38        outer += GAUSS_W[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 Iqxy(double qa, double qb, double qc, 
     48Iqabc(double qa, double qb, double qc, 
    4949    double sld, 
    5050    double sld_solvent, 

  • sasmodels/kerneldll.py

    r2d81cfe r1ddb794  
    185185        subprocess.check_output(command, shell=shell, stderr=subprocess.STDOUT) 
    186186    except subprocess.CalledProcessError as exc: 
    187         raise RuntimeError("compile failed.\n%s\n%s"%(command_str, exc.output)) 
     187        raise RuntimeError("compile failed.\n%s\n%s" 
     188                           % (command_str, exc.output.decode())) 
    188189    if not os.path.exists(output): 
    189190        raise RuntimeError("compile failed.  File is in %r"%source) 
Note: See TracChangeset for help on using the changeset viewer.