Changes in / [d5014e4:6ceca44] in sasmodels


Ignore:
Files:
13 edited

Legend:

Unmodified
Added
Removed
  • .gitignore

    r67cc0ff re9ed2de  
    2424/example/Fit_*/ 
    2525/example/batch_fit.csv 
     26# Note: changes to gauss20, gauss76 and gauss150 are still tracked since they 
     27# are part of the repo and required by some models. 
    2628/sasmodels/models/lib/gauss*.c 
  • .travis.yml

    rce8c388 r2d09df1  
    66    env: 
    77    - PY=2.7 
     8    - DEPLOY=True 
    89  - os: linux 
    910    env: 
     
    5152  on: 
    5253    branch: master 
     54    condition: $DEPLOY = True 
    5355notifications: 
    5456  slack: 
  • appveyor.yml

    rd810d96 r1258e32  
    6767    - cmd: conda install --yes --quiet obvious-ci 
    6868    - cmd: conda install --yes --quiet numpy toolchain scipy cython cffi 
    69     - cmd: conda install --yes --channel conda-forge pyopencl 
     69    #- cmd: conda install --yes --channel conda-forge pyopencl 
    7070    - cmd: pip install bumps unittest-xml-reporting tinycc 
    7171 
  • doc/developer/overview.rst

    r3d40839 r2a7e20e  
    185185jitter applied before particle orientation. 
    186186 
     187When computing the orientation dispersity integral, the weights for 
     188the individual points depends on the map projection used to translate jitter 
     189angles into latitude/longitude.  The choice of projection is set by 
     190*sasmodels.generate.PROJECTION*, with the default *PROJECTION=1* for 
     191equirectangular and *PROJECTION=2* for sinusoidal.  The more complicated 
     192Guyou and Postel projections are not implemented. See explore.jitter.draw_mesh 
     193for details. 
     194 
    187195For numerical integration within form factors etc. sasmodels is mostly using 
    188196Gaussian quadrature with 20, 76 or 150 points depending on the model. It also 
     
    199207Useful testing routines include: 
    200208 
    201 :mod:`asymint` a direct implementation of the surface integral for certain 
    202 models to get a more trusted value for the 1D integral using a 
    203 reimplementation of the 2D kernel in python and mpmath (which computes math 
    204 functions to arbitrary precision). It uses $\theta$ ranging from 0 to $\pi$ 
    205 and $\phi$ ranging from 0 to $2\pi$. It perhaps would benefit from including 
    206 the U-substitution for $\theta$. 
    207  
    208 :mod:`check1d` uses sasmodels 1D integration and compares that with a 
     209The *sascomp* utility is used to view and compare models with different 
     210parameters and calculation engines. The usual case is to simply plot a 
     211model that you are developing:: 
     212 
     213    python sascomp path/to/model.py 
     214 
     215Once the obvious problems are addressed, check the numerical precision 
     216across a variety of randomly generated inputs:: 
     217 
     218    python sascomp -engine=single,double path/to/model.py -sets=10 
     219 
     220You can compare different parameter values for the same or different models. 
     221For example when looking along the long axis of a cylinder ($\theta=0$), 
     222dispersity in $\theta$ should be equivalent to dispersity in $\phi$ 
     223when $\phi=90$:: 
     224 
     225    python sascomp -2d cylinder theta=0 phi=0,90 theta_pd_type=rectangle \\ 
     226    phi_pd_type=rectangle phi_pd=10,1 theta_pd=1,10 length=500 radius=10 
     227 
     228It turns out that they are not because the equirectangular map projection 
     229weights the points by $\cos(\theta)$ so $\Delta\theta$ is not identical 
     230to $\Delta\phi$.  Setting *PROJECTION=2* in :mod:`sasmodels.generate` helps 
     231somewhat.  Postel would help even more in this case, though leading 
     232to distortions elsewhere.  See :mod:`sasmodels.compare` for many more details. 
     233 
     234*sascomp -ngauss=n* allows you to set the number of quadrature points used 
     235for the 1D integration for any model.  For example, a carbon nanotube with 
     236length 10 $\mu$\ m and radius 1 nm is not computed correctly at high $q$:: 
     237 
     238    python sascomp cylinder length=100000 radius=10 -ngauss=76,500 -double -highq 
     239 
     240Note: ticket 702 gives some forms for long rods and thin disks that may 
     241be used in place of the integration, depending on $q$, radius and length; if 
     242the cylinder model is updated with these corrections then above call will show 
     243no difference. 
     244 
     245*explore/check1d.py* uses sasmodels 1D integration and compares that with a 
    209246rectangle distribution in $\theta$ and $\phi$, with $\theta$ limits set to 
    210247$\pm 90/\sqrt(3)$ and $\phi$ limits set to $\pm 180/\sqrt(3)$ [The rectangle 
     
    214251similar reasoning.] This should rotate the sample through the entire 
    215252$\theta$-$\phi$ surface according to the pattern that you see in jitter.py when 
    216 you modify it to use 'rectangle' rather than 'gaussian' for its distribution 
    217 without changing the viewing angle. In order to match the 1-D pattern for 
    218 an arbitrary viewing angle on triaxial shapes, we need to integrate 
     253you use 'rectangle' rather than 'gaussian' for its distribution without 
     254changing the viewing angle. In order to match the 1-D pattern for an arbitrary 
     255viewing angle on triaxial shapes, we need to integrate 
    219256over $\theta$, $\phi$ and $\Psi$. 
    220257 
    221 When computing the dispersity integral, weights are scaled by 
    222 $|\cos(\delta \theta)|$ to account for the points in $\phi$ getting closer 
    223 together as $\delta \theta$ increases. 
    224 [This will probably change so that instead of adjusting the weights, we will 
    225 adjust $\delta\theta$-$\delta\phi$ mesh so that the point density in 
    226 $\delta\phi$ is lower at larger $\delta\theta$. The flag USE_SCALED_PHI in 
    227 *kernel_iq.c* selects an alternative algorithm.] 
    228  
    229 The integrated dispersion is computed at a set of $(qx, qy)$ points $(q 
    230 \cos(\alpha), q \sin(\alpha))$ at some angle $\alpha$ (currently angle=0) for 
    231 each $q$ used in the 1-D integration. The individual $q$ points should be 
    232 equivalent to asymint rect-n when the viewing angle is set to 
    233 $(\theta,\phi,\Psi) = (90, 0, 0)$. Such tests can help to validate that 2d 
    234 models are consistent with 1d models. 
    235  
    236 :mod:`sascomp -sphere=n` uses the same rectangular distribution as check1d to 
    237 compute the pattern of the $q_x$-$q_y$ grid. 
    238  
    239 The :mod:`sascomp` utility can be used for 2d as well as 1d calculations to 
    240 compare results for two sets of parameters or processor types, for example 
    241 these two oriented cylinders here should be equivalent. 
    242  
    243 :mod:`\./sascomp -2d cylinder theta=0 phi=0,90 theta_pd_type=rectangle phi_pd_type=rectangle phi_pd=10,1 theta_pd=1,10 length=500 radius=10` 
    244  
     258*sascomp -sphere=n* uses the same rectangular distribution as check1d to 
     259compute the pattern of the $q_x$-$q_y$ grid.  This ought to produce a 
     260spherically symmetric pattern. 
     261 
     262*explore/precision.py* investigates the accuracy of individual functions 
     263on the different execution platforms.  This includes the basic special 
     264functions as well as more complex expressions used in scattering.  In many 
     265cases the OpenCL function in sasmodels will use a polynomial approximation 
     266over part of the range to improve accuracy, as shown in:: 
     267 
     268    python explore/precision.py 3j1/x 
     269 
     270*explore/realspace.py* allows you to set up a Monte Carlo simulation of your 
     271model by sampling random points within and computing the 1D and 2D scattering 
     272patterns.  This was used to check the core shell parallelepiped example.  This 
     273is similar to the general sas calculator in sasview, though it uses different 
     274code. 
     275 
     276*explore/jitter.py* is for exploring different options for handling 
     277orientation and orientation dispersity.  It uses *explore/guyou.py* to 
     278generate the Guyou projection. 
     279 
     280*explore/asymint.py* is a direct implementation of the 1D integration for 
     281a number of different models.  It calculates the integral for a particular 
     282$q$ using several different integration schemes, including mpmath with 100 
     283bits of precision (33 digits), so you can use it to check the target values 
     284for the 1D model tests.  This is not a general purpose tool; you will need to 
     285edit the file to change the model parameters. It does not currently 
     286apply the $u=cos(\theta)$ substitution that many models are using 
     287internally so the 76-point gaussian quadrature may not match the value 
     288produced by the eqivalent model from sasmodels. 
     289 
     290*explore/symint.py* is for rotationally symmetric models (just cylinder for 
     291now), so it can compute an entire curve rather than a single $q$ point.  It 
     292includes code to compare the long cylinder approximation to cylinder. 
     293 
     294*explore/rpa.py* is for checking different implementations of the RPA model 
     295against calculations for specific blends.  This is a work in (suspended) 
     296progress. 
     297 
     298There are a few modules left over in *explore* that can probably be removed. 
     299*angular_pd.py* was an early version of *jitter.py*.  *sc.py* and *sc.c* 
     300was used to explore different calculations for paracrystal models; it has 
     301been absorbed into *asymint.py*. *transform_angles.py* translates orientation 
     302parameters from the SasView 3.x definition to sasmodels. 
     303 
     304*explore/angles.py* generates code for the view and jitter transformations. 
     305Keep this around since it may be needed if we add new projections. 
    245306 
    246307Testing 
  • explore/precision.py

    r2a602c7 r2a7e20e  
    345345) 
    346346add_function( 
    347     name="(1/2+(1-cos(x))/x^2-sin(x)/x)/x", 
     347    name="(1/2-sin(x)/x+(1-cos(x))/x^2)/x", 
    348348    mp_function=lambda x: (0.5 - mp.sin(x)/x + (1-mp.cos(x))/(x*x))/x, 
    349349    np_function=lambda x: (0.5 - np.sin(x)/x + (1-np.cos(x))/(x*x))/x, 
     
    609609    names = ", ".join(sorted(ALL_FUNCTIONS)) 
    610610    print("""\ 
    611 usage: precision.py [-f/a/r] [-x<range>] name... 
     611usage: precision.py [-f/a/r] [-x<range>] "name" ... 
    612612where 
    613613    -f indicates that the function value should be plotted, 
     
    620620      zoom indicates linear stepping in [1000, 1010] 
    621621      neg indicates linear stepping in [-100.1, 100.1] 
    622 and name is "all [first]" or one of: 
     622and name is "all" or one of: 
    623623    """+names) 
    624624    sys.exit(1) 
  • explore/realspace.py

    ra1c32c2 r8fb2a94  
    4444    """ 
    4545    return Rx(phi)*Ry(theta)*Rz(psi) 
     46 
     47def apply_view(points, view): 
     48    """ 
     49    Apply the view transform (theta, phi, psi) to a set of points. 
     50 
     51    Points are stored in a 3 x n numpy array. 
     52 
     53    View angles are in degrees. 
     54    """ 
     55    theta, phi, psi = view 
     56    return np.asarray((Rz(phi)*Ry(theta)*Rz(psi))*np.matrix(points.T)).T 
     57 
     58 
     59def invert_view(qx, qy, view): 
     60    """ 
     61    Return (qa, qb, qc) for the (theta, phi, psi) view angle at detector 
     62    pixel (qx, qy). 
     63 
     64    View angles are in degrees. 
     65    """ 
     66    theta, phi, psi = view 
     67    q = np.vstack((qx.flatten(), qy.flatten(), 0*qx.flatten())) 
     68    return np.asarray((Rz(-psi)*Ry(-theta)*Rz(-phi))*np.matrix(q)) 
     69 
    4670 
    4771class Shape: 
     
    219243        values = self.value.repeat(points.shape[0]) 
    220244        return values, self._adjust(points) 
     245 
     246NUMBA = False 
     247if NUMBA: 
     248    from numba import njit 
     249    @njit("f8[:](f8[:],f8[:],f8[:],f8[:],f8[:],f8[:],f8[:])") 
     250    def _Iqxy(values, x, y, z, qa, qb, qc): 
     251        Iq = np.zeros_like(qa) 
     252        for j in range(len(Iq)): 
     253            total = 0. + 0j 
     254            for k in range(len(Iq)): 
     255                total += values[k]*np.exp(1j*(qa[j]*x[k] + qb[j]*y[k] + qc[j]*z[k])) 
     256            Iq[j] = abs(total)**2 
     257        return Iq 
     258 
     259def calc_Iqxy(qx, qy, rho, points, volume=1.0, view=(0, 0, 0)): 
     260    qx, qy = np.broadcast_arrays(qx, qy) 
     261    qa, qb, qc = invert_view(qx, qy, view) 
     262    rho, volume = np.broadcast_arrays(rho, volume) 
     263    values = rho*volume 
     264    x, y, z = points.T 
     265 
     266    # I(q) = |sum V(r) rho(r) e^(1j q.r)|^2 / sum V(r) 
     267    if NUMBA: 
     268        Iq = _Iqxy(values, x, y, z, qa.flatten(), qb.flatten(), qc.flatten()) 
     269    else: 
     270        Iq = [abs(np.sum(values*np.exp(1j*(qa_k*x + qb_k*y + qc_k*z))))**2 
     271              for qa_k, qb_k, qc_k in zip(qa.flat, qb.flat, qc.flat)] 
     272    return np.asarray(Iq).reshape(qx.shape) / np.sum(volume) 
    221273 
    222274def _calc_Pr_nonuniform(r, rho, points): 
     
    239291            print("processing %d of %d"%(k, len(rho)-1)) 
    240292    Pr = extended_Pr[1:-1] 
    241     return Pr / Pr.max() 
    242  
    243 def calc_Pr(r, rho, points): 
    244     # P(r) with uniform steps in r is 3x faster; check if we are uniform 
    245     # before continuing 
    246     if np.max(np.abs(np.diff(r) - r[0])) > r[0]*0.01: 
    247         return _calc_Pr_nonuniform(r, rho, points) 
    248  
     293    return Pr 
     294 
     295def _calc_Pr_uniform(r, rho, points): 
    249296    # Make Pr a little be bigger than necessary so that only distances 
    250297    # min < d < max end up in Pr 
    251     n_max = len(r) 
     298    dr, n_max = r[0], len(r) 
    252299    extended_Pr = np.zeros(n_max+1, 'd') 
    253300    t0 = time.clock() 
    254301    t_next = t0 + 3 
    255     row_zero = np.zeros(len(rho), 'i') 
    256302    for k, rho_k in enumerate(rho[:-1]): 
    257303        distances = sqrt(np.sum((points[k] - points[k+1:])**2, axis=1)) 
    258304        weights = rho_k * rho[k+1:] 
    259         index = np.minimum(np.asarray(distances/r[0], 'i'), n_max) 
     305        index = np.minimum(np.asarray(distances/dr, 'i'), n_max) 
    260306        # Note: indices may be duplicated, so "Pr[index] += w" will not work!! 
    261307        extended_Pr += np.bincount(index, weights, n_max+1) 
     
    269315    # we are only accumulating the upper triangular distances. 
    270316    #Pr = Pr * 2 / len(rho)**2 
    271     return Pr / Pr.max() 
     317    return Pr 
    272318 
    273319    # Can get an additional 2x by going to C.  Cuda/OpenCL will allow even 
     
    291337} 
    292338""" 
     339 
     340if NUMBA: 
     341    @njit("f8[:](f8[:], f8[:], f8[:,:])") 
     342    def _calc_Pr_uniform_jit(r, rho, points): 
     343        dr = r[0] 
     344        n_max = len(r) 
     345        Pr = np.zeros_like(r) 
     346        for j in range(len(rho) - 1): 
     347            x, y, z = points[j, 0], points[j, 1], points[j, 2] 
     348            for k in range(j+1, len(rho)): 
     349                distance = sqrt((x - points[k, 0])**2 
     350                                + (y - points[k, 1])**2 
     351                                + (z - points[k, 2])**2) 
     352                index = int(distance/dr) 
     353                if index < n_max: 
     354                    Pr[index] += rho[j] * rho[k] 
     355        # Make Pr independent of sampling density.  The factor of 2 comes because 
     356        # we are only accumulating the upper triangular distances. 
     357        #Pr = Pr * 2 / len(rho)**2 
     358        return Pr 
     359 
     360 
     361def calc_Pr(r, rho, points): 
     362    # P(r) with uniform steps in r is 3x faster; check if we are uniform 
     363    # before continuing 
     364    if np.max(np.abs(np.diff(r) - r[0])) > r[0]*0.01: 
     365        Pr = _calc_Pr_nonuniform(r, rho, points) 
     366    else: 
     367        if NUMBA: 
     368            Pr = _calc_Pr_uniform_jit(r, rho, points) 
     369        else: 
     370            Pr = _calc_Pr_uniform(r, rho, points) 
     371    return Pr / Pr.max() 
     372 
    293373 
    294374def j0(x): 
     
    333413        plt.legend() 
    334414 
     415def plot_calc_2d(qx, qy, Iqxy, theory=None): 
     416    import matplotlib.pyplot as plt 
     417    qx, qy = bin_edges(qx), bin_edges(qy) 
     418    #qx, qy = np.meshgrid(qx, qy) 
     419    if theory is not None: 
     420        plt.subplot(121) 
     421    plt.pcolormesh(qx, qy, np.log10(Iqxy)) 
     422    plt.xlabel('qx (1/A)') 
     423    plt.ylabel('qy (1/A)') 
     424    if theory is not None: 
     425        plt.subplot(122) 
     426        plt.pcolormesh(qx, qy, np.log10(theory)) 
     427        plt.xlabel('qx (1/A)') 
     428 
    335429def plot_points(rho, points): 
    336430    import mpl_toolkits.mplot3d 
     
    343437        pass 
    344438    n = len(points) 
    345     index = np.random.choice(n, size=1000) if n > 1000 else slice(None, None) 
     439    #print("len points", n) 
     440    index = np.random.choice(n, size=500) if n > 500 else slice(None, None) 
    346441    ax.scatter(points[index, 0], points[index, 1], points[index, 2], c=rho[index]) 
    347442    #low, high = points.min(axis=0), points.max(axis=0) 
    348443    #ax.axis([low[0], high[0], low[1], high[1], low[2], high[2]]) 
    349444    ax.autoscale(True) 
     445 
     446def check_shape(shape, fn=None): 
     447    rho_solvent = 0 
     448    q = np.logspace(-3, 0, 200) 
     449    r = shape.r_bins(q, r_step=0.01) 
     450    sampling_density = 6*5000 / shape.volume() 
     451    rho, points = shape.sample(sampling_density) 
     452    t0 = time.time() 
     453    Pr = calc_Pr(r, rho-rho_solvent, points) 
     454    print("calc Pr time", time.time() - t0) 
     455    Iq = calc_Iq(q, r, Pr) 
     456    theory = (q, fn(q)) if fn is not None else None 
     457 
     458    import pylab 
     459    #plot_points(rho, points); pylab.figure() 
     460    plot_calc(r, Pr, q, Iq, theory=theory) 
     461    pylab.show() 
     462 
     463def check_shape_2d(shape, fn=None, view=(0, 0, 0)): 
     464    rho_solvent = 0 
     465    nq, qmax = 100, 1.0 
     466    qx = np.linspace(0.0, qmax, nq) 
     467    qy = np.linspace(0.0, qmax, nq) 
     468    Qx, Qy = np.meshgrid(qx, qy) 
     469    sampling_density = 50000 / shape.volume() 
     470    #t0 = time.time() 
     471    rho, points = shape.sample(sampling_density) 
     472    #print("sample time", time.time() - t0) 
     473    t0 = time.time() 
     474    Iqxy = calc_Iqxy(Qx, Qy, rho, points, view=view) 
     475    print("calc time", time.time() - t0) 
     476    theory = fn(Qx, Qy) if fn is not None else None 
     477    Iqxy += 0.001 * Iqxy.max() 
     478    if theory is not None: 
     479        theory += 0.001 * theory.max() 
     480 
     481    import pylab 
     482    #plot_points(rho, points); pylab.figure() 
     483    plot_calc_2d(qx, qy, Iqxy, theory=theory) 
     484    pylab.show() 
     485 
     486def sas_sinx_x(x): 
     487    with np.errstate(all='ignore'): 
     488        retvalue = sin(x)/x 
     489    retvalue[x == 0.] = 1. 
     490    return retvalue 
    350491 
    351492def sas_2J1x_x(x): 
     
    373514    Iq = Iq/Iq[0] 
    374515    return Iq 
     516 
     517def cylinder_Iqxy(qx, qy, radius, length, view=(0, 0, 0)): 
     518    qa, qb, qc = invert_view(qx, qy, view) 
     519    qab = np.sqrt(qa**2 + qb**2) 
     520    Fq = sas_2J1x_x(qab*radius) * j0(qc*length/2) 
     521    Iq = Fq**2 
     522    return Iq.reshape(qx.shape) 
    375523 
    376524def sphere_Iq(q, radius): 
     
    415563    return Iq/Iq[0] 
    416564 
    417 def check_shape(shape, fn=None): 
    418     rho_solvent = 0 
    419     q = np.logspace(-3, 0, 200) 
    420     r = shape.r_bins(q, r_step=0.01) 
    421     sampling_density = 15000 / shape.volume() 
    422     rho, points = shape.sample(sampling_density) 
    423     Pr = calc_Pr(r, rho-rho_solvent, points) 
    424     Iq = calc_Iq(q, r, Pr) 
    425     theory = (q, fn(q)) if fn is not None else None 
    426  
    427     import pylab 
    428     #plot_points(rho, points); pylab.figure() 
    429     plot_calc(r, Pr, q, Iq, theory=theory) 
    430     pylab.show() 
     565def csbox_Iqxy(qx, qy, a, b, c, da, db, dc, slda, sldb, sldc, sld_core, view=(0,0,0)): 
     566    qa, qb, qc = invert_view(qx, qy, view) 
     567 
     568    sld_solvent = 0 
     569    overlapping = False 
     570    dr0 = sld_core - sld_solvent 
     571    drA, drB, drC = slda-sld_solvent, sldb-sld_solvent, sldc-sld_solvent 
     572    tA, tB, tC = a + 2*da, b + 2*db, c + 2*dc 
     573    siA = a*sas_sinx_x(a*qa/2) 
     574    siB = b*sas_sinx_x(b*qb/2) 
     575    siC = c*sas_sinx_x(c*qc/2) 
     576    siAt = tA*sas_sinx_x(tA*qa/2) 
     577    siBt = tB*sas_sinx_x(tB*qb/2) 
     578    siCt = tC*sas_sinx_x(tC*qc/2) 
     579    Fq = (dr0*siA*siB*siC 
     580          + drA*(siAt-siA)*siB*siC 
     581          + drB*siA*(siBt-siB)*siC 
     582          + drC*siA*siB*(siCt-siC)) 
     583    Iq = Fq**2 
     584    return Iq.reshape(qx.shape) 
    431585 
    432586def check_cylinder(radius=25, length=125, rho=2.): 
     
    434588    fn = lambda q: cylinder_Iq(q, radius, length) 
    435589    check_shape(shape, fn) 
     590 
     591def check_cylinder_2d(radius=25, length=125, rho=2., view=(0, 0, 0)): 
     592    shape = EllipticalCylinder(radius, radius, length, rho) 
     593    fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view) 
     594    check_shape_2d(shape, fn, view=view) 
     595 
     596def check_cylinder_2d_lattice(radius=25, length=125, rho=2., 
     597                              view=(0, 0, 0)): 
     598    nx, dx = 1, 2*radius 
     599    ny, dy = 30, 2*radius 
     600    nz, dz = 30, length 
     601    dx, dy, dz = 2*dx, 2*dy, 2*dz 
     602    def center(*args): 
     603        sigma = 0.333 
     604        space = 2 
     605        return [(space*n+np.random.randn()*sigma)*x for n, x in args] 
     606    shapes = [EllipticalCylinder(radius, radius, length, rho, 
     607                                 #center=(ix*dx, iy*dy, iz*dz) 
     608                                 orientation=np.random.randn(3)*0, 
     609                                 center=center((ix, dx), (iy, dy), (iz, dz)) 
     610                                ) 
     611              for ix in range(nx) 
     612              for iy in range(ny) 
     613              for iz in range(nz)] 
     614    shape = Composite(shapes) 
     615    fn = lambda qx, qy, view=view: cylinder_Iqxy(qx, qy, radius, length, view=view) 
     616    check_shape_2d(shape, fn, view=view) 
    436617 
    437618def check_sphere(radius=125, rho=2): 
     
    449630    side_c2 = copy(side_c).shift(0, 0, -c-dc) 
    450631    shape = Composite((core, side_a, side_b, side_c, side_a2, side_b2, side_c2)) 
    451     fn = lambda q: csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 
    452     check_shape(shape, fn) 
     632    def fn(q): 
     633        return csbox_Iq(q, a, b, c, da, db, dc, slda, sldb, sldc, sld_core) 
     634    #check_shape(shape, fn) 
     635 
     636    view = (20, 30, 40) 
     637    def fn_xy(qx, qy): 
     638        return csbox_Iqxy(qx, qy, a, b, c, da, db, dc, 
     639                          slda, sldb, sldc, sld_core, view=view) 
     640    check_shape_2d(shape, fn_xy, view=view) 
    453641 
    454642if __name__ == "__main__": 
    455643    check_cylinder(radius=10, length=40) 
     644    #check_cylinder_2d(radius=10, length=40, view=(90,30,0)) 
     645    #check_cylinder_2d_lattice(radius=10, length=50, view=(90,30,0)) 
    456646    #check_sphere() 
    457647    #check_csbox() 
  • sasmodels/compare.py

    r2d81cfe r2a7e20e  
    9292    -accuracy=Low accuracy of the resolution calculation Low, Mid, High, Xhigh 
    9393    -neval=1 sets the number of evals for more accurate timing 
     94    -ngauss=0 overrides the number of points in the 1-D gaussian quadrature 
    9495 
    9596    === precision options === 
  • sasmodels/kernel_iq.c

    r108e70e rec8d4ac  
    174174static double 
    175175qac_apply( 
    176     QACRotation rotation, 
     176    QACRotation *rotation, 
    177177    double qx, double qy, 
    178178    double *qa_out, double *qc_out) 
    179179{ 
    180     const double dqc = rotation.R31*qx + rotation.R32*qy; 
     180    const double dqc = rotation->R31*qx + rotation->R32*qy; 
    181181    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
    182182    const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy); 
     
    247247static double 
    248248qabc_apply( 
    249     QABCRotation rotation, 
     249    QABCRotation *rotation, 
    250250    double qx, double qy, 
    251251    double *qa_out, double *qb_out, double *qc_out) 
    252252{ 
    253     *qa_out = rotation.R11*qx + rotation.R12*qy; 
    254     *qb_out = rotation.R21*qx + rotation.R22*qy; 
    255     *qc_out = rotation.R31*qx + rotation.R32*qy; 
     253    *qa_out = rotation->R11*qx + rotation->R12*qy; 
     254    *qb_out = rotation->R21*qx + rotation->R22*qy; 
     255    *qc_out = rotation->R31*qx + rotation->R32*qy; 
    256256} 
    257257 
     
    454454  // theta, phi, dtheta, dphi are defined below in projection to avoid repeated code. 
    455455  #define BUILD_ROTATION() qac_rotation(&rotation, theta, phi, dtheta, dphi); 
    456   #define APPLY_ROTATION() qac_apply(rotation, qx, qy, &qa, &qc) 
     456  #define APPLY_ROTATION() qac_apply(&rotation, qx, qy, &qa, &qc) 
    457457  #define CALL_KERNEL() CALL_IQ_AC(qa, qc, local_values.table) 
    458458 
     
    468468  local_values.table.psi = 0.; 
    469469  #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi) 
    470   #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc) 
     470  #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc) 
    471471  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
    472472#elif defined(CALL_IQ_XY) 
     
    479479#endif 
    480480 
    481 // Doing jitter projection code outside the previous if block so that we don't 
    482 // need to repeat the identical logic in the IQ_AC and IQ_ABC branches.  This 
    483 // will become more important if we implement more projections, or more 
    484 // complicated projections. 
    485 #if defined(CALL_IQ) || defined(CALL_IQ_A) 
     481// Define APPLY_PROJECTION depending on model symmetries. We do this outside 
     482// the previous if block so that we don't need to repeat the identical 
     483// logic in the IQ_AC and IQ_ABC branches.  This will become more important 
     484// if we implement more projections, or more complicated projections. 
     485#if defined(CALL_IQ) || defined(CALL_IQ_A)  // no orientation 
    486486  #define APPLY_PROJECTION() const double weight=weight0 
    487 #elif defined(CALL_IQ_XY) 
     487#elif defined(CALL_IQ_XY) // pass orientation to the model 
    488488  // CRUFT: support oriented model which define Iqxy rather than Iqac or Iqabc 
    489489  // Need to plug the values for the orientation angles back into parameter 
     
    515515    #define APPLY_PROJECTION() const double weight=weight0 
    516516  #endif 
    517 #else // !spherosymmetric projection 
     517#else // apply jitter and view before calling the model 
    518518  // Grab the "view" angles (theta, phi, psi) from the initial parameter table. 
    519519  const double theta = values[details->theta_par+2]; 
     
    526526  // we go through the mesh. 
    527527  double dtheta, dphi, weight; 
    528   #if PROJECTION == 1 
     528  #if PROJECTION == 1 // equirectangular 
    529529    #define APPLY_PROJECTION() do { \ 
    530530      dtheta = local_values.table.theta; \ 
     
    532532      weight = fabs(cos(dtheta*M_PI_180)) * weight0; \ 
    533533    } while (0) 
    534   #elif PROJECTION == 2 
     534  #elif PROJECTION == 2 // sinusoidal 
    535535    #define APPLY_PROJECTION() do { \ 
    536536      dtheta = local_values.table.theta; \ 
     
    542542    } while (0) 
    543543  #endif 
    544 #endif // !spherosymmetric projection 
     544#endif // done defining APPLY_PROJECTION 
    545545 
    546546// ** define looping macros ** 
  • sasmodels/models/core_shell_parallelepiped.c

    r108e70e re077231  
    4343    // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c 
    4444    // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) 
    45     // Code is rewritten,the code is compliant with Diva Singhs thesis now (Dirk Honecker) 
    46     // Code rewritten (PAK) 
     45    // Code is rewritten, the code is compliant with Diva Singh's thesis now (Dirk Honecker) 
     46    // Code rewritten; cross checked against hollow rectangular prism and realspace (PAK) 
    4747 
    4848    const double half_q = 0.5*q; 
     
    121121    const double drC = crim_sld-solvent_sld; 
    122122 
    123     // The definitions of ta, tb, tc are not the same as in the 1D case because there is no 
    124     // the scaling by B. 
    125123    const double tA = length_a + 2.0*thick_rim_a; 
    126124    const double tB = length_b + 2.0*thick_rim_b; 
  • sasmodels/models/core_shell_parallelepiped.py

    r10ee838 r97be877  
    136136* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    137137* **Converted to sasmodels by:** Miguel Gonzales **Date:** February 26, 2016 
    138 * **Last Modified by:** Wojciech Potrzebowski **Date:** January 11, 2017 
    139 * **Currently Under review by:** Paul Butler 
     138* **Last Modified by:** Paul Kienzle **Date:** October 17, 2017 
     139* Cross-checked against hollow rectangular prism and rectangular prism for 
     140  equal thickness overlapping sides, and by Monte Carlo sampling of points 
     141  within the shape for non-uniform, non-overlapping sides. 
    140142""" 
    141143 
  • sasmodels/models/lib/gauss150.c

    r74768cb r99b84ec  
    1 /* 
    2  *  GaussWeights.c 
    3  *  SANSAnalysis 
    4  * 
    5  *  Created by Andrew Jackson on 4/23/07. 
    6  *  Copyright 2007 __MyCompanyName__. All rights reserved. 
    7  * 
    8  */ 
     1// Created by Andrew Jackson on 4/23/07 
     2 
    93 #ifdef GAUSS_N 
    104 # undef GAUSS_N 
     
    1610 #define GAUSS_W Gauss150Wt 
    1711 
    18 // Note: using array size 152 so that it is a multiple of 4 
    1912 
    20 // Gaussians 
     13// Note: using array size 152 rather than 150 so that it is a multiple of 4. 
     14// Some OpenCL devices prefer that vectors start and end on nice boundaries. 
    2115constant double Gauss150Z[152]={ 
    2216        -0.9998723404457334, 
     
    170164        0.9993274305065947, 
    171165        0.9998723404457334, 
    172         0., 
    173         0. 
     166        0., // zero padding is ignored 
     167        0.  // zero padding is ignored 
    174168}; 
    175169 
     
    325319        0.0007624720924706, 
    326320        0.0003276086705538, 
    327         0., 
    328         0. 
     321        0., // zero padding is ignored 
     322        0.  // zero padding is ignored 
    329323}; 
  • sasmodels/models/lib/gauss20.c

    r74768cb r99b84ec  
    1 /* 
    2  *  GaussWeights.c 
    3  *  SANSAnalysis 
    4  * 
    5  *  Created by Andrew Jackson on 4/23/07. 
    6  *  Copyright 2007 __MyCompanyName__. All rights reserved. 
    7  * 
    8  */ 
     1// Created by Andrew Jackson on 4/23/07 
     2 
    93 #ifdef GAUSS_N 
    104 # undef GAUSS_N 
  • sasmodels/models/lib/gauss76.c

    r74768cb r99b84ec  
    1 /* 
    2  *  GaussWeights.c 
    3  *  SANSAnalysis 
    4  * 
    5  *  Created by Andrew Jackson on 4/23/07. 
    6  *  Copyright 2007 __MyCompanyName__. All rights reserved. 
    7  * 
    8  */ 
     1// Created by Andrew Jackson on 4/23/07 
     2 
    93 #ifdef GAUSS_N 
    104 # undef GAUSS_N 
Note: See TracChangeset for help on using the changeset viewer.