Changes in / [d5cc54c:15ec718] in sasmodels


Ignore:
Files:
4 added
2 deleted
15 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/direct_model.py

    rb397165 r0444c02  
    3131from . import resolution2d 
    3232from .details import make_kernel_args, dispersion_mesh 
     33from sas.sasgui.perspectives.fitting.fitpage import FitPage 
    3334 
    3435try: 
    … …  
    193194        # interpret data 
    194195        if hasattr(data, 'lam'): 
     196        #if not FitPage.no_transform.GetValue(): #if the no_transform radio button is not active DOES NOT WORK! not active before fitting 
    195197            self.data_type = 'sesans' 
    196198        elif hasattr(data, 'qx_data'): 

  • sasmodels/kernel_header.c

    rb7e8b94 r218cdbc  
    148148inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; } 
    149149 
    150 #if 1 
    151 //think cos(theta) should be sin(theta) in new coords, RKH 11Jan2017 
     150#if 0 
    152151#define ORIENT_SYMMETRIC(qx, qy, theta, phi, q, sn, cn) do { \ 
    153152    SINCOS(phi*M_PI_180, sn, cn); \ 
    154153    q = sqrt(qx*qx + qy*qy); \ 
    155     cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * sin(theta*M_PI_180));  \ 
     154    cn  = (q==0. ? 1.0 : (cn*qx + sn*qy)/q * cos(theta*M_PI_180));  \ 
    156155    sn = sqrt(1 - cn*cn); \ 
    157156    } while (0) 

  • sasmodels/kernelpy.py

    rb397165 r46d9f48  
    1515from .generate import F64 
    1616from .kernel import KernelModel, Kernel 
     17import gc 
    1718 
    1819try: 
    … …  
    7980        """ 
    8081        self.q = None 
     82        gc.collect() 
    8183 
    8284class PyKernel(Kernel): 

  • sasmodels/models/barbell.py

    rfcb33e4 r0d6e865  
    2020.. math:: 
    2121 
    22     I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q,\alpha).sin(\alpha)\right> 
     22    I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q)\right> 
    2323 
    24 where the amplitude $A(q,\alpha)$ with the rod axis at angle $\alpha$ to $q$ is given as 
     24where the amplitude $A(q)$ is given as 
    2525 
    2626.. math:: 
    2727 
    2828    A(q) =&\ \pi r^2L 
    29         \frac{\sin\left(\tfrac12 qL\cos\alpha\right)} 
    30              {\tfrac12 qL\cos\alpha} 
    31         \frac{2 J_1(qr\sin\alpha)}{qr\sin\alpha} \\ 
     29        \frac{\sin\left(\tfrac12 qL\cos\theta\right)} 
     30             {\tfrac12 qL\cos\theta} 
     31        \frac{2 J_1(qr\sin\theta)}{qr\sin\theta} \\ 
    3232        &\ + 4 \pi R^3 \int_{-h/R}^1 dt 
    33         \cos\left[ q\cos\alpha 
     33        \cos\left[ q\cos\theta 
    3434            \left(Rt + h + {\tfrac12} L\right)\right] 
    3535        \times (1-t^2) 
    36         \frac{J_1\left[qR\sin\alpha \left(1-t^2\right)^{1/2}\right]} 
    37              {qR\sin\alpha \left(1-t^2\right)^{1/2}} 
     36        \frac{J_1\left[qR\sin\theta \left(1-t^2\right)^{1/2}\right]} 
     37             {qR\sin\theta \left(1-t^2\right)^{1/2}} 
    3838 
    3939The $\left<\ldots\right>$ brackets denote an average of the structure over 
    40 all orientations. $\left<A^2(q,\alpha)\right>$ is then the form factor, $P(q)$. 
     40all orientations. $\left<A^2(q)\right>$ is then the form factor, $P(q)$. 
    4141The scale factor is equivalent to the volume fraction of cylinders, each of 
    4242volume, $V$. Contrast $\Delta\rho$ is the difference of scattering length 
    … …  
    8585* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    8686* **Last Modified by:** Paul Butler **Date:** March 20, 2016 
    87 * **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 
     87* **Last Reviewed by:** Paul Butler **Date:** March 20, 2016 
    8888""" 
    8989from numpy import inf 

  • sasmodels/models/capped_cylinder.py

    rfcb33e4 r0d6e865  
    2121.. math:: 
    2222 
    23     I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q,\alpha).sin(\alpha)\right> 
     23    I(q) = \frac{\Delta \rho^2}{V} \left<A^2(q)\right> 
    2424 
    25 where the amplitude $A(q,\alpha)$ with the rod axis at angle $\alpha$ to $q$ is given as 
     25where the amplitude $A(q)$ is given as 
    2626 
    2727.. math:: 
    2828 
    2929    A(q) =&\ \pi r^2L 
    30         \frac{\sin\left(\tfrac12 qL\cos\alpha\right)} 
    31             {\tfrac12 qL\cos\alpha} 
    32         \frac{2 J_1(qr\sin\alpha)}{qr\sin\alpha} \\ 
     30        \frac{\sin\left(\tfrac12 qL\cos\theta\right)} 
     31            {\tfrac12 qL\cos\theta} 
     32        \frac{2 J_1(qr\sin\theta)}{qr\sin\theta} \\ 
    3333        &\ + 4 \pi R^3 \int_{-h/R}^1 dt 
    34         \cos\left[ q\cos\alpha 
     34        \cos\left[ q\cos\theta 
    3535            \left(Rt + h + {\tfrac12} L\right)\right] 
    3636        \times (1-t^2) 
    37         \frac{J_1\left[qR\sin\alpha \left(1-t^2\right)^{1/2}\right]} 
    38              {qR\sin\alpha \left(1-t^2\right)^{1/2}} 
     37        \frac{J_1\left[qR\sin\theta \left(1-t^2\right)^{1/2}\right]} 
     38             {qR\sin\theta \left(1-t^2\right)^{1/2}} 
    3939 
    4040The $\left<\ldots\right>$ brackets denote an average of the structure over 
    … …  
    8888* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    8989* **Last Modified by:** Paul Butler **Date:** September 30, 2016 
    90 * **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 
    91  
     90* **Last Reviewed by:** Richard Heenan **Date:** March 19, 2016 
    9291""" 
    9392from numpy import inf 

  • sasmodels/models/core_shell_bicelle.py

    rfcb33e4 rb829b16  
    4141 
    4242    I(Q,\alpha) = \frac{\text{scale}}{V_t} \cdot 
    43         F(Q,\alpha)^2.sin(\alpha) + \text{background} 
     43        F(Q,\alpha)^2 + \text{background} 
    4444 
    4545where 
    … …  
    8585* **Author:** NIST IGOR/DANSE **Date:** pre 2010 
    8686* **Last Modified by:** Paul Butler **Date:** September 30, 2016 
    87 * **Last Reviewed by:** Richard Heenan **Date:** January 4, 2017 
     87* **Last Reviewed by:** Richard Heenan **Date:** October 5, 2016 
    8888""" 
    8989 

  • sasmodels/models/core_shell_cylinder.py

    rfcb33e4 r755ecc2  
    99.. math:: 
    1010 
    11     I(q,\alpha) = \frac{\text{scale}}{V_s} F^2(q,\alpha).sin(\alpha) + \text{background} 
     11    I(q,\alpha) = \frac{\text{scale}}{V_s} F^2(q) + \text{background} 
    1212 
    1313where 
    … …  
    1515.. math:: 
    1616 
    17     F(q,\alpha) = &\ (\rho_c - \rho_s) V_c 
     17    F(q) = &\ (\rho_c - \rho_s) V_c 
    1818           \frac{\sin \left( q \tfrac12 L\cos\alpha \right)} 
    1919                {q \tfrac12 L\cos\alpha} 

  • sasmodels/models/core_shell_ellipsoid.py

    r7c2935c r73e08ae  
    162162 
    163163q = 0.1 
    164 # tests had in old coords theta=0, phi=0; new coords theta=90, phi=0 
    165 qx = q*cos(pi/6.0) 
    166 qy = q*sin(pi/6.0) 
    167 # 11Jan2017 RKH sorted tests after redefinition of angles 
    168 tests = [ 
    169      # Accuracy tests based on content in test/utest_coreshellellipsoidXTmodel.py 
    170     [{'radius_equat_core': 200.0, 
    171       'x_core': 0.1, 
    172       'thick_shell': 50.0, 
    173       'x_polar_shell': 0.2, 
    174       'sld_core': 2.0, 
    175       'sld_shell': 1.0, 
    176       'sld_solvent': 6.3, 
    177       'background': 0.001, 
    178       'scale': 1.0, 
    179      }, 1.0, 0.00189402], 
     164phi = pi/6 
     165qx = q*cos(phi) 
     166qy = q*sin(phi) 
     167# After redefinition of angles find new reasonable values for unit test 
     168#tests = [ 
     169#    # Accuracy tests based on content in test/utest_coreshellellipsoidXTmodel.py 
     170#    [{'radius_equat_core': 200.0, 
     171#      'x_core': 0.1, 
     172#      'thick_shell': 50.0, 
     173#      'x_polar_shell': 0.2, 
     174#      'sld_core': 2.0, 
     175#      'sld_shell': 1.0, 
     176#      'sld_solvent': 6.3, 
     177#      'background': 0.001, 
     178#      'scale': 1.0, 
     179#     }, 1.0, 0.00189402], 
    180180 
    181181    # Additional tests with larger range of parameters 
    182     [{'background': 0.01}, 0.1, 11.6915], 
    183  
    184     [{'radius_equat_core': 20.0, 
    185       'x_core': 200.0, 
    186       'thick_shell': 54.0, 
    187       'x_polar_shell': 3.0, 
    188       'sld_core': 20.0, 
    189       'sld_shell': 10.0, 
    190       'sld_solvent': 6.0, 
    191       'background': 0.0, 
    192       'scale': 1.0, 
    193      }, 0.01, 8688.53], 
    194  
    195    # 2D tests 
    196    [{'background': 0.001, 
    197      'theta': 90.0, 
    198      'phi': 0.0, 
    199      }, (0.4, 0.5), 0.00690673], 
    200  
    201    [{'radius_equat_core': 20.0, 
    202       'x_core': 200.0, 
    203       'thick_shell': 54.0, 
    204       'x_polar_shell': 3.0, 
    205       'sld_core': 20.0, 
    206       'sld_shell': 10.0, 
    207       'sld_solvent': 6.0, 
    208       'background': 0.01, 
    209       'scale': 0.01, 
    210       'theta': 90.0, 
    211       'phi': 0.0, 
    212      }, (qx, qy), 0.01000025], 
    213     ] 
     182#    [{'background': 0.01}, 0.1, 11.6915], 
     183 
     184#    [{'radius_equat_core': 20.0, 
     185#      'x_core': 200.0, 
     186#      'thick_shell': 54.0, 
     187#      'x_polar_shell': 3.0, 
     188#      'sld_core': 20.0, 
     189#      'sld_shell': 10.0, 
     190#      'sld_solvent': 6.0, 
     191#      'background': 0.0, 
     192#      'scale': 1.0, 
     193#     }, 0.01, 8688.53], 
     194 
     195#   [{'background': 0.001}, (0.4, 0.5), 0.00690673], 
     196 
     197#   [{'radius_equat_core': 20.0, 
     198#      'x_core': 200.0, 
     199#      'thick_shell': 54.0, 
     200#      'x_polar_shell': 3.0, 
     201#      'sld_core': 20.0, 
     202#      'sld_shell': 10.0, 
     203#      'sld_solvent': 6.0, 
     204#      'background': 0.01, 
     205#      'scale': 0.01, 
     206#     }, (qx, qy), 0.0100002], 
     207#    ] 

  • sasmodels/models/cylinder.py

    rb7e8b94 r4cdd0cc  
    22# Note: model title and parameter table are inserted automatically 
    33r""" 
    4  
     4The form factor is normalized by the particle volume V = \piR^2L. 
    55For information about polarised and magnetic scattering, see 
    66the :ref:`magnetism` documentation. 
    … …  
    1414.. math:: 
    1515 
    16     P(q,\alpha) = \frac{\text{scale}}{V} F^2(q,\alpha).sin(\alpha) + \text{background} 
     16    P(q,\alpha) = \frac{\text{scale}}{V} F^2(q,\alpha) + \text{background} 
    1717 
    1818where 
    … …  
    2525           \frac{J_1 \left(q R \sin \alpha\right)}{q R \sin \alpha} 
    2626 
    27 and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, $V =\pi R^2L$ 
     27and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, $V$ 
    2828is the volume of the cylinder, $L$ is the length of the cylinder, $R$ is the 
    2929radius of the cylinder, and $\Delta\rho$ (contrast) is the scattering length 
    … …  
    3535.. math:: 
    3636 
    37     F^2(q)=\int_{0}^{\pi/2}{F^2(q,\alpha)\sin(\alpha)d\alpha}=\int_{0}^{1}{F^2(q,u)du} 
     37    F^2(q)=\int_{0}^{\pi/2}{F^2(q,\theta)\sin(\theta)d\theta} 
    3838 
    3939 
    40 Numerical integration is simplified by a change of variable to $u = cos(\alpha)$ with  
    41 $sin(\alpha)=\sqrt{1-u^2}$.  
    42  
    43 The output of the 1D scattering intensity function for randomly oriented 
    44 cylinders is thus given by 
    45  
    46 .. math:: 
    47  
    48     P(q) = \frac{\text{scale}}{V} 
    49         \int_0^{\pi/2} F^2(q,\alpha) \sin \alpha\ d\alpha + \text{background} 
    50  
    51  
    52 NB: The 2nd virial coefficient of the cylinder is calculated based on the 
    53 radius and length values, and used as the effective radius for $S(q)$ 
    54 when $P(q) \cdot S(q)$ is applied. 
    55  
    56 For oriented cylinders, we define the direction of the 
    57 axis of the cylinder using two angles $\theta$ (note this is not the 
    58 same as the scattering angle used in q) and $\phi$. Those angles 
     40To provide easy access to the orientation of the cylinder, we define the 
     41axis of the cylinder using two angles $\theta$ and $\phi$. Those angles 
    5942are defined in :numref:`cylinder-angle-definition` . 
    6043 
    … …  
    6548    Definition of the angles for oriented cylinders. 
    6649 
    67 The $\theta$ and $\phi$ parameters only appear in the model when fitting 2d data. 
     50 
     51NB: The 2nd virial coefficient of the cylinder is calculated based on the 
     52radius and length values, and used as the effective radius for $S(q)$ 
     53when $P(q) \cdot S(q)$ is applied. 
     54 
     55The output of the 1D scattering intensity function for randomly oriented 
     56cylinders is then given by 
     57 
     58.. math:: 
     59 
     60    P(q) = \frac{\text{scale}}{V} 
     61        \int_0^{\pi/2} F^2(q,\alpha) \sin \alpha\ d\alpha + \text{background} 
     62 
     63The $\theta$ and $\phi$ parameters are not used for the 1D output. 
    6864 
    6965Validation 
    … …  
    7874 
    7975    P(q) = \int_0^{\pi/2} d\phi 
    80         \int_0^\pi p(\theta) P_0(q,\theta) \sin \theta\ d\theta 
     76        \int_0^\pi p(\alpha) P_0(q,\alpha) \sin \alpha\ d\alpha 
    8177 
    8278 
    8379where $p(\theta,\phi) = 1$ is the probability distribution for the orientation 
    84 and $P_0(q,\theta)$ is the scattering intensity for the fully oriented 
     80and $P_0(q,\alpha)$ is the scattering intensity for the fully oriented 
    8581system, and then comparing to the 1D result. 
    8682 
    … …  
    149145 
    150146qx, qy = 0.2 * np.cos(2.5), 0.2 * np.sin(2.5) 
    151 # After redefinition of angles, find new tests values.  Was 10 10 in old coords 
    152 tests = [[{}, 0.2, 0.042761386790780453], 
    153         [{}, [0.2], [0.042761386790780453]], 
    154 #  new coords     
    155         [{'theta':80.1534480601659, 'phi':10.1510817110481}, (qx, qy), 0.03514647218513852], 
    156         [{'theta':80.1534480601659, 'phi':10.1510817110481}, [(qx, qy)], [0.03514647218513852]], 
    157 # old coords   [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852], 
    158 #              [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]], 
    159         ] 
     147# After redefinition of angles, find new tests values  
     148#tests = [[{}, 0.2, 0.042761386790780453], 
     149#         [{}, [0.2], [0.042761386790780453]], 
     150#         [{'theta':10.0, 'phi':10.0}, (qx, qy), 0.03514647218513852], 
     151#         [{'theta':10.0, 'phi':10.0}, [(qx, qy)], [0.03514647218513852]], 
     152#        ] 
    160153del qx, qy  # not necessary to delete, but cleaner 
    161154# ADDED by:  RKH  ON: 18Mar2016 renamed sld's etc 

  • sasmodels/models/elliptical_cylinder.py

    rfcb33e4 ra807206  
    2020 
    2121    I(\vec q)=\frac{1}{V_\text{cyl}}\int{d\psi}\int{d\phi}\int{ 
    22         p(\theta,\phi,\psi)F^2(\vec q,\alpha,\psi)\sin(\alpha)d\alpha} 
     22        p(\theta,\phi,\psi)F^2(\vec q,\alpha,\psi)\sin(\theta)d\theta} 
    2323 
    2424with the functions 
    … …  
    2626.. math:: 
    2727 
    28     F(q,\alpha,\psi) = 2\frac{J_1(a)\sin(b)}{ab} 
     28    F(\vec q,\alpha,\psi) = 2\frac{J_1(a)\sin(b)}{ab} 
    2929 
    3030where 
    … …  
    3232.. math:: 
    3333 
    34     a = qr'\sin(\alpha) 
    35      
    36     b = q\frac{L}{2}\cos(\alpha) 
    37      
    38     r'=\frac{r_{minor}}{\sqrt{2}}\sqrt{(1+\nu^{2}) + (1-\nu^{2})cos(\psi)} 
     34    a &= \vec q\sin(\alpha)\left[ 
     35        r^2_\text{major}\sin^2(\psi)+r^2_\text{minor}\cos(\psi) \right]^{1/2} 
    3936 
     37    b &= \vec q\frac{L}{2}\cos(\alpha) 
    4038 
    41 and the angle $\psi$ is defined as the orientation of the major axis of the 
     39and the angle $\Psi$ is defined as the orientation of the major axis of the 
    4240ellipse with respect to the vector $\vec q$. The angle $\alpha$ is the angle 
    4341between the axis of the cylinder and $\vec q$. 
    … …  
    9795 
    9896L A Feigin and D I Svergun, *Structure Analysis by Small-Angle X-Ray and 
    99 Neutron Scattering*, Plenum, New York, (1987) [see table 3.4] 
    100  
    101 Authorship and Verification 
    102 ---------------------------- 
    103  
    104 * **Author:** 
    105 * **Last Modified by:**  
    106 * **Last Reviewed by:**  Richard Heenan - corrected equation in docs **Date:** December 21, 2016 
    107  
     97Neutron Scattering*, Plenum, New York, (1987) 
    10898""" 
    10999 

  • sasmodels/models/stacked_disks.py

    rb7e8b94 r07300ea  
    148148            sld_layer=0.0, 
    149149            sld_solvent=5.0, 
    150             theta=90, 
     150            theta=0, 
    151151            phi=0) 
    152 # After redefinition of spherical coordinates - 
    153 # tests had in old coords theta=0, phi=0; new coords theta=90, phi=0 
    154 # but should not matter here as so far all the tests are 1D not 2D 
    155 tests = [ 
    156 # Accuracy tests based on content in test/utest_extra_models.py. 
    157 # Added 2 tests with n_stacked = 5 using SasView 3.1.2 - PDB; for which alas q=0.001 values seem closer to n_stacked=1 not 5, changed assuming my 4.1 code OK, RKH 
    158     [{'thick_core': 10.0, 
    159       'thick_layer': 15.0, 
    160       'radius': 3000.0, 
    161       'n_stacking': 1.0, 
    162       'sigma_d': 0.0, 
    163       'sld_core': 4.0, 
    164       'sld_layer': -0.4, 
    165       'solvent_sd': 5.0, 
    166       'theta': 90.0, 
    167       'phi': 0.0, 
    168       'scale': 0.01, 
    169       'background': 0.001, 
    170      }, 0.001, 5075.12], 
    171     [{'thick_core': 10.0, 
    172       'thick_layer': 15.0, 
    173       'radius': 3000.0, 
    174       'n_stacking': 5, 
    175       'sigma_d': 0.0, 
    176       'sld_core': 4.0, 
    177       'sld_layer': -0.4, 
    178       'solvent_sd': 5.0, 
    179       'theta': 90.0, 
    180       'phi': 0.0, 
    181       'scale': 0.01, 
    182       'background': 0.001, 
    183 #     }, 0.001, 5065.12793824],    n_stacking=1 not 5 ? slight change in value here 11jan2017, check other cpu types 
    184 #     }, 0.001, 5075.11570493], 
    185      }, 0.001, 25325.635693], 
    186     [{'thick_core': 10.0, 
    187       'thick_layer': 15.0, 
    188       'radius': 3000.0, 
    189       'n_stacking': 5, 
    190       'sigma_d': 0.0, 
    191       'sld_core': 4.0, 
    192       'sld_layer': -0.4, 
    193       'solvent_sd': 5.0, 
    194       'theta': 90.0, 
    195       'phi': 0.0, 
    196       'scale': 0.01, 
    197       'background': 0.001, 
    198 #     }, 0.164, 0.0127673597265],    n_stacking=1 not 5 ?  slight change in value here 11jan2017, check other cpu types 
    199 #     }, 0.164, 0.01480812968], 
    200      }, 0.164, 0.0598367986509], 
    201  
    202     [{'thick_core': 10.0, 
    203       'thick_layer': 15.0, 
    204       'radius': 3000.0, 
    205       'n_stacking': 1.0, 
    206       'sigma_d': 0.0, 
    207       'sld_core': 4.0, 
    208       'sld_layer': -0.4, 
    209       'solvent_sd': 5.0, 
    210       'theta': 90.0, 
    211       'phi': 0.0, 
    212       'scale': 0.01, 
    213       'background': 0.001, 
    214 # second test here was at q=90, changed it to q=5, note I(q) is then just value of flat background 
    215      }, [0.001, 5.0], [5075.12, 0.001]], 
    216  
    217     [{'thick_core': 10.0, 
    218       'thick_layer': 15.0, 
    219       'radius': 3000.0, 
    220       'n_stacking': 1.0, 
    221       'sigma_d': 0.0, 
    222       'sld_core': 4.0, 
    223       'sld_layer': -0.4, 
    224       'solvent_sd': 5.0, 
    225       'theta': 90.0, 
    226       'phi': 0.0, 
    227       'scale': 0.01, 
    228       'background': 0.001, 
    229      }, ([0.4, 0.5]), [0.00105074, 0.00121761]], 
    230  
    231     [{'thick_core': 10.0, 
    232       'thick_layer': 15.0, 
    233       'radius': 3000.0, 
    234       'n_stacking': 1.0, 
    235       'sigma_d': 0.0, 
    236       'sld_core': 4.0, 
    237       'sld_layer': -0.4, 
    238      'solvent_sd': 5.0, 
    239       'theta': 90.0, 
    240       'phi': 0.0, 
    241       'scale': 0.01, 
    242       'background': 0.001, 
    243      }, ([1.3, 1.57]), [0.0010039, 0.0010038]], 
    244     ] 
    245 # 11Jan2017   RKH checking unit test agai, note they are all 1D, no 2D 
    246  
     152#After redefinition to spherical coordinates find new reasonable test values 
     153#tests = [ 
     154#    # Accuracy tests based on content in test/utest_extra_models.py. 
     155#    # Added 2 tests with n_stacked = 5 using SasView 3.1.2 - PDB 
     156#    [{'thick_core': 10.0, 
     157#      'thick_layer': 15.0, 
     158#      'radius': 3000.0, 
     159#      'n_stacking': 1.0, 
     160#      'sigma_d': 0.0, 
     161#      'sld_core': 4.0, 
     162#      'sld_layer': -0.4, 
     163#      'solvent_sd': 5.0, 
     164#      'theta': 0.0, 
     165#      'phi': 0.0, 
     166#      'scale': 0.01, 
     167#      'background': 0.001, 
     168#     }, 0.001, 5075.12], 
     169 
     170#    [{'thick_core': 10.0, 
     171#      'thick_layer': 15.0, 
     172#      'radius': 3000.0, 
     173#      'n_stacking': 5.0, 
     174#      'sigma_d': 0.0, 
     175#      'sld_core': 4.0, 
     176#      'sld_layer': -0.4, 
     177#      'solvent_sd': 5.0, 
     178#      'theta': 0.0, 
     179#      'phi': 0.0, 
     180#      'scale': 0.01, 
     181#      'background': 0.001, 
     182#     }, 0.001, 5065.12793824], 
     183 
     184#    [{'thick_core': 10.0, 
     185#      'thick_layer': 15.0, 
     186#      'radius': 3000.0, 
     187#      'n_stacking': 5.0, 
     188#      'sigma_d': 0.0, 
     189#      'sld_core': 4.0, 
     190#      'sld_layer': -0.4, 
     191#      'solvent_sd': 5.0, 
     192#      'theta': 0.0, 
     193#      'phi': 0.0, 
     194#      'scale': 0.01, 
     195#      'background': 0.001, 
     196#     }, 0.164, 0.0127673597265], 
     197 
     198#    [{'thick_core': 10.0, 
     199#      'thick_layer': 15.0, 
     200#      'radius': 3000.0, 
     201#      'n_stacking': 1.0, 
     202#      'sigma_d': 0.0, 
     203#      'sld_core': 4.0, 
     204#      'sld_layer': -0.4, 
     205#      'solvent_sd': 5.0, 
     206#      'theta': 0.0, 
     207#      'phi': 0.0, 
     208#      'scale': 0.01, 
     209#      'background': 0.001, 
     210#     }, [0.001, 90.0], [5075.12, 0.001]], 
     211 
     212#    [{'thick_core': 10.0, 
     213#      'thick_layer': 15.0, 
     214#      'radius': 3000.0, 
     215#      'n_stacking': 1.0, 
     216#      'sigma_d': 0.0, 
     217#      'sld_core': 4.0, 
     218#      'sld_layer': -0.4, 
     219#      'solvent_sd': 5.0, 
     220#      'theta': 0.0, 
     221#      'phi': 0.0, 
     222#      'scale': 0.01, 
     223#      'background': 0.001, 
     224#     }, ([0.4, 0.5]), [0.00105074, 0.00121761]], 
     225 
     226#    [{'thick_core': 10.0, 
     227#      'thick_layer': 15.0, 
     228#      'radius': 3000.0, 
     229#      'n_stacking': 1.0, 
     230#      'sigma_d': 0.0, 
     231#      'sld_core': 4.0, 
     232#      'sld_layer': -0.4, 
     233#     'solvent_sd': 5.0, 
     234#      'theta': 0.0, 
     235#      'phi': 0.0, 
     236#      'scale': 0.01, 
     237#      'background': 0.001, 
     238#     }, ([1.3, 1.57]), [0.0010039, 0.0010038]], 
     239#    ] 
     240# 21Mar2016   RKH notes that unit tests all have n_stacking=1, ought to test other values 
     241 

  • sasmodels/resolution.py

    rb397165 r26b848d  
    99from numpy import sqrt, log, log10, exp, pi  # type: ignore 
    1010import numpy as np  # type: ignore 
     11 
     12from sasmodels import sesans 
     13from sasmodels.sesans import SesansTransform as SesansTransform 
    1114 
    1215__all__ = ["Resolution", "Perfect1D", "Pinhole1D", "Slit1D", 
    … …  
    4346        raise NotImplementedError("Subclass does not define the apply function") 
    4447 
    45  
    4648class Perfect1D(Resolution): 
    4749    """ 
    … …  
    5557    def apply(self, theory): 
    5658        return theory 
    57  
    5859 
    5960class Pinhole1D(Resolution): 

  • sasmodels/resolution2d.py

    r7e94989 r40a87fa  
    6464        # just need q_calc and weights 
    6565        self.data = data 
    66         self.index = index if index is not None else slice(None) 
    67  
    68         self.qx_data = data.qx_data[self.index] 
    69         self.qy_data = data.qy_data[self.index] 
    70         self.q_data = data.q_data[self.index] 
     66        self.index = index 
     67 
     68        self.qx_data = data.qx_data[index] 
     69        self.qy_data = data.qy_data[index] 
     70        self.q_data = data.q_data[index] 
    7171 
    7272        dqx = getattr(data, 'dqx_data', None) 
    … …  
    7474        if dqx is not None and dqy is not None: 
    7575            # Here dqx and dqy mean dq_parr and dq_perp 
    76             self.dqx_data = dqx[self.index] 
    77             self.dqy_data = dqy[self.index] 
     76            self.dqx_data = dqx[index] 
     77            self.dqy_data = dqy[index] 
    7878            ## Remove singular points if exists 
    7979            self.dqx_data[self.dqx_data < SIGMA_ZERO] = SIGMA_ZERO 
    … …  
    126126        # The angle (phi) of the original q point 
    127127        q_phi = np.arctan(q_phi).repeat(nbins)\ 
    128             .reshape([nq, nbins]).transpose().flatten() 
     128            .reshape(nq, nbins).transpose().flatten() 
    129129        ## Find Gaussian weight for each dq bins: The weight depends only 
    130130        #  on r-direction (The integration may not need) 

  • sasmodels/sasview_model.py

    r64614ad r8977226  
    583583                            % type(qdist)) 
    584584 
    585     def get_composition_models(self): 
    586         """ 
    587             Returns usable models that compose this model 
    588         """ 
    589         s_model = None 
    590         p_model = None 
    591         if hasattr(self._model_info, "composition") \ 
    592            and self._model_info.composition is not None: 
    593             p_model = _make_model_from_info(self._model_info.composition[1][0])() 
    594             s_model = _make_model_from_info(self._model_info.composition[1][1])() 
    595         return p_model, s_model 
    596  
    597585    def calculate_Iq(self, qx, qy=None): 
    598586        # type: (Sequence[float], Optional[Sequence[float]]) -> np.ndarray 

  • sasmodels/sesans.py

    rb397165 r15ec718  
    1414import numpy as np  # type: ignore 
    1515from numpy import pi, exp  # type: ignore 
    16 from scipy.special import jv as besselj 
    17 #import direct_model.DataMixin as model 
     16from scipy.special import j0 
     17#from mpmath import j0 as j0 
    1818         
    19 def make_q(q_max, Rmax): 
    20     r""" 
    21     Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$ 
    22     to $q_max$. This is the integration range of the Hankel transform; bigger range and  
    23     more points makes a better numerical integration. 
    24     Smaller q_min will increase reliable spin echo length range.  
    25     Rmax is the "radius" of the largest expected object and can be set elsewhere. 
    26     q_max is determined by the acceptance angle of the SESANS instrument. 
    27     """ 
    28     from sas.sascalc.data_util.nxsunit import Converter 
     19class SesansTransform(object): 
     20    #: Set of spin-echo lengths in the measured data 
     21    SE = None  # type: np.ndarray 
     22    #: Maximum acceptance of scattering vector in the spin echo encoding dimension (for ToF: Q of min(R) and max(lam)) 
     23    zaccept = None # type: float 
     24    #: Maximum size sensitivity; larger radius requires more computation 
     25    Rmax = None  # type: float 
     26    #: q values to calculate when computing transform 
     27    q = None  # type: np.ndarray 
    2928 
    30     q_min = dq = 0.1 * 2*pi / Rmax 
    31     return np.arange(q_min, 
    32                      Converter(q_max[1])(q_max[0], 
    33                                          units="1/A"), 
    34                      dq) 
    35      
    36 def make_all_q(data): 
    37     """ 
    38     Return a $q$ vector suitable for calculating the total scattering cross section for 
    39     calculating the effect of finite acceptance angles on Time of Flight SESANS instruments. 
    40     If no acceptance is given, or unwanted (set "unwanted" flag in paramfile), no all_q vector is needed. 
    41     If the instrument has a rectangular acceptance, 2 all_q vectors are needed. 
    42     If the instrument has a circular acceptance, 1 all_q vector is needed 
    43      
    44     """ 
    45     if not data.has_no_finite_acceptance: 
    46         return [] 
    47     elif data.has_yz_acceptance(data): 
    48         # compute qx, qy 
    49         Qx, Qy = np.meshgrid(qx, qy) 
    50         return [Qx, Qy] 
    51     else: 
    52         # else only need q 
    53         # data.has_z_acceptance 
    54         return [q] 
     29    # transform arrays 
     30    _H = None  # type: np.ndarray 
     31    _H0 = None # type: np.ndarray 
    5532 
    56 def transform(data, q_calc, Iq_calc, qmono, Iq_mono): 
    57     """ 
    58     Decides which transform type is to be used, based on the experiment data file contents (header) 
    59     (2016-03-19: currently controlled from parameters script) 
    60     nqmono is the number of q vectors to be used for the detector integration 
    61     """ 
    62     nqmono = len(qmono) 
    63     if nqmono == 0: 
    64         result = call_hankel(data, q_calc, Iq_calc) 
    65     elif nqmono == 1: 
    66         q = qmono[0] 
    67         result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono) 
    68     else: 
    69         Qx, Qy = [qmono[0], qmono[1]] 
    70         Qx = np.reshape(Qx, nqx, nqy) 
    71         Qy = np.reshape(Qy, nqx, nqy) 
    72         Iq_mono = np.reshape(Iq_mono, nqx, nqy) 
    73         qx = Qx[0, :] 
    74         qy = Qy[:, 0] 
    75         result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono) 
     33    def set_transform(self, SE, zaccept, Rmax): 
     34        if self.SE is None or len(SE) != len(self.SE) or np.any(SE != self.SE) or zaccept != self.zaccept or Rmax != self.Rmax: 
     35            self.SE, self.zaccept, self.Rmax = SE, zaccept, Rmax 
     36            self._set_q() 
     37            self._set_hankel() 
    7638 
    77     return result 
     39    def apply(self, Iq): 
     40        G0 = np.dot(self._H0, Iq) 
     41        G = np.dot(self._H.T, Iq) 
     42        P = G - G0 
     43        return P 
    7844 
    79 def call_hankel(data, q_calc, Iq_calc): 
    80     return hankel((data.x, data.x_unit), 
    81                   (data.lam, data.lam_unit), 
    82                   (data.sample.thickness, 
    83                    data.sample.thickness_unit), 
    84                   q_calc, Iq_calc) 
    85    
    86 def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono): 
    87     return hankel(data.x, data.lam * 1e-9, 
    88                   data.sample.thickness / 10, 
    89                   q_calc, Iq_calc) 
    90                    
    91 def call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono): 
    92     return hankel(data.x, data.y, data.lam * 1e-9, 
    93                   data.sample.thickness / 10, 
    94                   q_calc, Iq_calc) 
    95                          
    96 def TotalScatter(model, parameters):  #Work in progress!! 
    97 #    Calls a model with existing model parameters already in place, then integrate the product of q and I(q) from 0 to (4*pi/lambda) 
    98     allq = np.linspace(0,4*pi/wavelength,1000) 
    99     allIq = 1 
    100     integral = allq*allIq 
    101      
     45    def _set_q(self): 
     46        #q_min = dq = 0.1 * 2*pi / self.Rmax 
    10247 
     48        q_max = 2*pi / (self.SE[1]-self.SE[0]) 
     49        q_min = dq = 0.1 *2*pi / (np.size(self.SE) * self.SE[-1]) 
    10350 
    104 def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still 
    105 #============================================================================== 
    106 #     2D Cosine Transform if "wavelength" is a vector 
    107 #============================================================================== 
    108 #allq is the q-space needed to create the total scattering cross-section 
     51        #q_min = dq = q_max / 100000 
     52        q=np.arange(q_min, q_max, q_min) 
     53        self.q = q 
     54        self.dq = dq 
    10955 
    110     Gprime = np.zeros_like(wavelength, 'd') 
    111     s = np.zeros_like(wavelength, 'd') 
    112     sd = np.zeros_like(wavelength, 'd') 
    113     Gprime = np.zeros_like(wavelength, 'd') 
    114     f = np.zeros_like(wavelength, 'd') 
    115     for i, wavelength_i in enumerate(wavelength): 
    116         z = magfield*wavelength_i 
    117         allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
    118         allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
    119         alldq = (allq[1]-allq[0])*1e10 
    120         sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
    121         s[i]=1-exp(-sigma) 
    122         for j, Iqy_j, qy_j in enumerate(qy): 
    123             for k, Iqz_k, qz_k in enumerate(qz): 
    124                 Iq = np.sqrt(Iqy_j^2+Iqz_k^2) 
    125                 q = np.sqrt(qy_j^2 + qz_k^2) 
    126                 Gintegral = Iq*cos(z*Qz_k) 
    127                 Gprime[i] += Gintegral 
    128 #                sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i] 
    129 #                s[i] += 1-exp(Totalscatter(modelname)*thickness) 
    130 #                For now, work with standard 2-phase scatter 
     56    def _set_hankel(self): 
     57        #Rmax = #value in text box somewhere in FitPage? 
     58        q = self.q 
     59        dq = self.dq 
     60        SElength = self.SE 
    13161 
     62        H0 = dq / (2 * pi) * q 
     63        q=np.array(q,dtype='float32') 
     64        SElength=np.array(SElength,dtype='float32') 
    13265 
    133                 sd[i] += Iq 
    134         f[i] = 1-s[i]+sd[i] 
    135         P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i] 
     66        # Using numpy tile, dtype is conserved 
     67        repq=np.tile(q,(SElength.size,1)) 
     68        repSE=np.tile(SElength,(q.size,1)) 
     69        H = dq / (2 * pi) * j0(repSE*repq.T)*repq.T 
    13670 
     71        # Using numpy meshgrid - meshgrid produces float64 from float32 inputs! Problem for 32-bit OS: Memerrors! 
     72        #H0 = dq / (2 * pi) * q 
     73        #repSE, repq = np.meshgrid(SElength, q) 
     74        #repq=np.array(repq,dtype='float32') 
     75        #repSE=np.array(repSE,dtype='float32') 
     76        #H = dq / (2 * pi) * j0(repSE*repq)*repq 
    13777 
     78        self._H, self._H0 = H, H0 
    13879 
    139  
    140 def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname): 
    141 #============================================================================== 
    142 #     HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS 
    143 #============================================================================== 
    144 #acceptq is the q-space needed to create limited acceptance effect 
    145     SElength= wavelength*magfield 
    146     G = np.zeros_like(SElength, 'd') 
    147     threshold=2*pi*theta/wavelength 
    148     for i, SElength_i in enumerate(SElength): 
    149         allq=np.linspace() #for calculating the Q-range of the  scattering power integral 
    150         allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow 
    151         alldq = (allq[1]-allq[0])*1e10 
    152         sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) 
    153         s[i]=1-exp(-sigma) 
    154  
    155         dq = (q[1]-q[0])*1e10 
    156         a = (x<threshold) 
    157         acceptq = a*q 
    158         acceptIq = a*Iq 
    159  
    160         G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq) 
    161  
    162 #        G[i]=np.sum(integral) 
    163  
    164     G *= dq*1e10*2*pi 
    165  
    166     P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0])) 
    167      
    168 def hankel(SElength, wavelength, thickness, q, Iq): 
    169     r""" 
    170     Compute the expected SESANS polarization for a given SANS pattern. 
    171  
    172     Uses the hankel transform followed by the exponential.  The values for *zz* 
    173     (or spin echo length, or delta), wavelength and sample thickness should 
    174     come from the dataset.  $q$ should be chosen such that the oscillations 
    175     in $I(q)$ are well sampled (e.g., $5 \cdot 2 \pi/d_{\max}$). 
    176  
    177     *SElength* [A] is the set of $z$ points at which to compute the 
    178     Hankel transform 
    179  
    180     *wavelength* [m]  is the wavelength of each individual point *zz* 
    181  
    182     *thickness* [cm] is the sample thickness. 
    183  
    184     *q* [A$^{-1}$] is the set of $q$ points at which the model has been 
    185     computed. These should be equally spaced. 
    186  
    187     *I* [cm$^{-1}$] is the value of the SANS model at *q* 
    188     """ 
    189  
    190     from sas.sascalc.data_util.nxsunit import Converter 
    191     wavelength = Converter(wavelength[1])(wavelength[0],"A") 
    192     thickness = Converter(thickness[1])(thickness[0],"A") 
    193     Iq = Converter("1/cm")(Iq,"1/A") # All models default to inverse centimeters 
    194     SElength = Converter(SElength[1])(SElength[0],"A") 
    195  
    196     G = np.zeros_like(SElength, 'd') 
    197 #============================================================================== 
    198 #     Hankel Transform method if "wavelength" is a scalar; mono-chromatic SESANS 
    199 #============================================================================== 
    200     for i, SElength_i in enumerate(SElength): 
    201         integral = besselj(0, q*SElength_i)*Iq*q 
    202         G[i] = np.sum(integral) 
    203     G0 = np.sum(Iq*q) 
    204  
    205     # [m^-1] step size in q, needed for integration 
    206     dq = (q[1]-q[0]) 
    207  
    208     # integration step, convert q into [m**-1] and 2 pi circle integration 
    209     G *= dq*2*pi 
    210     G0 = np.sum(Iq*q)*dq*2*np.pi 
    211  
    212     P = exp(thickness*wavelength**2/(4*pi**2)*(G-G0)) 
    213  
    214     return P 
     80class SESANS1D(SesansTransform): 
     81    def __init__(self, data, _H0, _H, q_calc): 
     82        # x values of the data (Sasview requires these to be named "q") 
     83        self.q = data.x 
     84        self._H0 = _H0 
     85        self._H = _H 
     86        # Pysmear does some checks on the smearer object, these checks require the "data" object... 
     87        self.data=data 
     88        # q values of the SAS model 
     89        self.q_calc = q_calc # These are the MODEL's q values used by the smearer (in this case: the Hankel transform) 
     90    def apply(self, theory): 
     91        return SesansTransform.apply(self,theory) 
Note: See TracChangeset for help on using the changeset viewer.